/********************************************************************/
/**                                                                **/
/**                         LINEAR ALGEBRA                         **/
/**                         (second part)                          **/
/**                                                                **/
/********************************************************************/
/* $Id: alglin2.c,v 1.3 1999/09/23 17:50:55 karim Exp $ */
#include "pari.h"

/*******************************************************************/
/*                                                                 */
/*                   CHARACTERISTIC POLYNOMIAL                     */
/*                                                                 */
/*******************************************************************/

GEN
charpoly0(GEN x, int v, long flag)
{
  if (v<0) v = 0;
  switch(flag)
  {
    case 0: return caradj(x,v,0);
    case 1: return caract(x,v);
    case 2: return carhess(x,v);
  }
  err(flagerr,"charpoly"); return NULL; /* not reached */
}


static GEN
caract2_i(GEN p, GEN x, int v, GEN (subres_f)(GEN,GEN,GEN*))
{
  long av = avma, d;
  GEN p1, p2 = leading_term(p);

  if (!signe(x)) return gpowgs(polx[v], lgef(p)-3);
  x = gneg_i(x); x[2] = ladd((GEN)x[2], polx[MAXVARN]);
  p1=subres_f(p, x, NULL);
  if (varn(p1)==MAXVARN) setvarn(p1,v); else p1=gsubst(p1,MAXVARN,polx[v]);

  if (!gcmp1(p2) && (d=lgef(x)-3) > 0) p1 = gdiv(p1, gpuigs(p2,d));
  return gerepileupto(av,p1);
}

/* return caract(Mod(x,p)) in variable v */
GEN
caract2(GEN p, GEN x, int v)
{
  return caract2_i(p,x,v, subresall);
}
GEN
caractducos(GEN p, GEN x, int v)
{
  return caract2_i(p,x,v, (GEN (*)(GEN,GEN,GEN*))resultantducos);
}

/* characteristic pol. Easy cases. Return NULL in case it's not so easy. */
static GEN
easychar(GEN x, int v, GEN *py)
{
  long l,tetpil,lx;
  GEN p1,p2;

  switch(typ(x))
  {
    case t_INT: case t_REAL: case t_INTMOD:
    case t_FRAC: case t_FRACN: case t_PADIC:
      p1=cgetg(4,t_POL);
      p1[1]=evalsigne(1) | evallgef(4) | evalvarn(v);
      p1[2]=lneg(x); p1[3]=un;
      if (py)
      {
	p2=cgetg(2,t_MAT);
	p2[1]=lgetg(2,t_COL);
	coeff(p2,1,1)=lcopy(x);
	*py=p2;
      }
      return p1;

    case t_COMPLEX: case t_QUAD:
      if (py) err(typeer,"easychar");
      p1=cgetg(5,t_POL);
      p1[1]=evalsigne(1) | evallgef(5) | evalvarn(v);
      p1[2]=lnorm(x); l=avma; p2=gtrace(x); tetpil=avma;
      p1[3]=lpile(l,tetpil,gneg(p2));
      p1[4]=un; return p1;

    case t_POLMOD:
      if (py) err(typeer,"easychar");
      return caract2((GEN)x[1], (GEN)x[2], v);

    case t_MAT:
      lx=lg(x);
      if (lx==1)
      {
	if (py) *py = cgetg(1,t_MAT);
	return polun[v];
      }
      if (lg(x[1]) != lx) break;
      return NULL;
  }
  err(mattype1,"easychar");
  return NULL; /* not reached */
}

GEN
caract(GEN x, int v)
{
  long n,k,l=avma,tetpil;
  GEN  p1,p2,p3,p4,p5,p6;

  if ((p1 = easychar(x,v,NULL))) return p1;

  p1=gzero; p2=gun;
  n=lg(x)-1; if (n&1) p2=gneg_i(p2);
  p4=cgetg(3,t_RFRACN); p5=dummycopy(polx[v]); p4[2]=(long)p5;
  p6=cgeti(3); p6[1] = evalsigne(-1) | evallgefint(3);
  for (k=0; k<=n; k++)
  {
    p3=det(gsub(gscalmat(stoi(k),n), x));
    p4[1]=lmul(p3,p2); p6[2]=k;
    p1=gadd(p4,p1); p5[2]=(long)p6;
    if (k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
  }
  p2=mpfact(n); tetpil=avma;
  return gerepile(l,tetpil,gdiv((GEN) p1[1],p2));
}

GEN
caradj0(GEN x, long v)
{
  return caradj(x,v,NULL);
}

/* Using traces: return the characteristic polynomial of x (in variable v).
 * If py != NULL, the adjoint matrix is put there.
 */
GEN
caradj(GEN x, long v, GEN *py)
{
  long i,j,k,l,av,tetpil;
  GEN p,y,t,*gptr[2];

  if ((p = easychar(x,v,py))) return p;

  l=lg(x);
  if (l==1) { p=polun[v]; if (py) *py=gcopy(x); return p; }
  if (l==2) { p=gsub(polx[v],gtrace(x)); if (py) *py=idmat(1); return p; }

  p=cgetg(l+2,t_POL); p[1]=evalsigne(1) | evallgef(l+2) | evalvarn(v);
  av=avma; t=gtrace(x); tetpil=avma;
  t=gerepile(av,tetpil,gneg(t));
  p[l]=(long)t; p[l+1]=un;

  av=avma; y=cgetg(l,t_MAT);
  for (j=1; j<l; j++)
  {
    y[j]=lgetg(l,t_COL);
    for (i=1; i<l; i++)
      coeff(y,i,j) = (i==j) ? ladd(gcoeff(x,i,j),t) : coeff(x,i,j);
  }

  for (k=2; k<l-1; k++)
  {
    GEN z=gmul(x,y);

    t=gtrace(z); tetpil=avma;
    t=gdivgs(t,-k); y=cgetg(l,t_MAT);
    for (j=1; j<l; j++)
    {
      y[j]=lgetg(l,t_COL);
      for (i=1;i<l;i++)
	coeff(y,i,j) = (i==j)? ladd(gcoeff(z,i,i),t): lcopy(gcoeff(z,i,j));
    }
    gptr[0]=&y; gptr[1]=&t;
    gerepilemanysp(av,tetpil,gptr,2);
    p[l-k+1]=(long)t; av=avma;
  }

  t=gzero;
  for (i=1; i<l; i++)
    t=gadd(t,gmul(gcoeff(x,1,i),gcoeff(y,i,1)));
  tetpil=avma; t=gneg(t);

  if (py)
  {
    *py = (l&1) ? gneg(y) : gcopy(y);
    gptr[0] = &t; gptr[1] = py;
    gerepilemanysp(av,tetpil,gptr,2);
    p[2]=(long)t;
  }
  else p[2]=lpile(av,tetpil,t);
  i = gvar2(p);
  if (i == v) err(talker,"incorrect variable in caradj");
  if (i < v) p = poleval(p, polx[v]); 
  return p;
}

GEN
adj(GEN x)
{
  GEN y;
  caradj(x,MAXVARN,&y); return y;
}

/*******************************************************************/
/*                                                                 */
/*                       HESSENBERG FORM                           */
/*                                                                 */
/*******************************************************************/

GEN
hess(GEN x)
{
  long lx=lg(x),av=avma,tetpil,m,i,j;
  GEN p,p1,p2;

  if (typ(x) != t_MAT) err(mattype1,"hess");
  if (lx==1) return cgetg(1,t_MAT);
  if (lg(x[1])!=lx) err(mattype1,"hess");
  x = dummycopy(x);

  for (m=2; m<lx-1; m++)
    for (i=m+1; i<lx; i++)
    {
      p = gcoeff(x,i,m-1);
      if (!gcmp0(p))
      {
	for (j=m-1; j<lx; j++)
	{
	  p1 = gcoeff(x,i,j);
	  coeff(x,i,j) = coeff(x,m,j);
	  coeff(x,m,j) = (long)p1;
	}
	p1=(GEN)x[i]; x[i]=x[m]; x[m]=(long)p1;
	for (i=m+1; i<lx; i++)
        {
          p1 = gcoeff(x,i,m-1);
	  if (!gcmp0(p1))
	  {
	    p1 = gdiv(p1,p); p2 = gneg_i(p1);
            coeff(x,i,m-1) = zero;
	    for (j=m; j<lx; j++)
	      coeff(x,i,j) = ladd(gcoeff(x,i,j), gmul(p2,gcoeff(x,m,j)));
	    for (j=1; j<lx; j++)
	      coeff(x,j,m) = ladd(gcoeff(x,j,m), gmul(p1,gcoeff(x,j,i)));
	  }
        }
	break;
      }
    }
  tetpil=avma; return gerepile(av,tetpil,gcopy(x));
}

GEN
carhess(GEN x, long v)
{
  long av,tetpil,lx,r,i;
  GEN *y,p1,p2,p3,p4;

  if ((p1 = easychar(x,v,NULL))) return p1;

  lx=lg(x); av=avma; y = (GEN*) new_chunk(lx);
  y[0]=polun[v]; p1=hess(x); p2=polx[v];
  tetpil=avma;
  for (r=1; r<lx; r++)
  {
    y[r]=gmul(y[r-1], gsub(p2,gcoeff(p1,r,r)));
    p3=gun; p4=gzero;
    for (i=r-1; i; i--)
    {
      p3=gmul(p3,gcoeff(p1,i+1,i));
      p4=gadd(p4,gmul(gmul(p3,gcoeff(p1,i,r)),y[i-1]));
    }
    tetpil=avma; y[r]=gsub(y[r],p4);
  }
  return gerepile(av,tetpil,y[lx-1]);
}

/*******************************************************************/
/*                                                                 */
/*                            NORM                                 */
/*                                                                 */
/*******************************************************************/

GEN
gnorm(GEN x)
{
  long l,lx,i,tetpil, tx=typ(x);
  GEN p1,p2,y;

  switch(tx)
  {
    case t_INT:
      return sqri(x);

    case t_REAL:
      return mulrr(x,x);

    case t_FRAC: case t_FRACN:
      return gsqr(x);

    case t_COMPLEX:
      l=avma; p1=gsqr((GEN) x[1]); p2=gsqr((GEN) x[2]);
      tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));

    case t_QUAD:
      l=avma; p1=(GEN)x[1];
      p2=gmul((GEN) p1[2], gsqr((GEN) x[3]));
      p1 = gcmp0((GEN) p1[3])? gsqr((GEN) x[2])
                             : gmul((GEN) x[2], gadd((GEN) x[2],(GEN) x[3]));
      tetpil=avma; return gerepile(l,tetpil,gadd(p1,p2));

    case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
      l=avma; p1=gmul(gconj(x),x); tetpil=avma;
      return gerepile(l,tetpil,greal(p1));

    case t_POLMOD:
      p1=(GEN)x[1]; p2=leading_term(p1);
      if (gcmp1(p2) || gcmp0((GEN) x[2])) return subres(p1,(GEN) x[2]);
      l=avma; y=subres(p1,(GEN)x[2]); p1=gpuigs(p2,lgef(x[2])-3);
      tetpil=avma; return gerepile(l,tetpil,gdiv(y,p1));

    case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=lnorm((GEN) x[i]);
      return y;
  }
  err(typeer,"gnorm");
  return NULL; /* not reached */
}

GEN
gnorml2(GEN x)
{
  GEN y;
  long i,tx=typ(x),lx,av;

  if (! is_matvec_t(tx)) return gnorm(x);
  lx=lg(x); if (lx==1) return gzero;

  av=avma; y = gnorml2((GEN) x[1]);
  for (i=2; i<lx; i++)
     y = gadd(gnorml2((GEN) x[i]), y);
  return gerepileupto(av,y);
}

/*******************************************************************/
/*                                                                 */
/*                          CONJUGATION                            */
/*                                                                 */
/*******************************************************************/

GEN
gconj(GEN x)
{
  long lx,i,tx=typ(x);
  GEN z;

  switch(tx)
  {
    case t_INT: case t_REAL: case t_INTMOD:
    case t_FRAC: case t_FRACN: case t_PADIC:
      return gcopy(x);

    case t_COMPLEX:
      z=cgetg(3,t_COMPLEX);
      z[1]=lcopy((GEN) x[1]);
      z[2]=lneg((GEN) x[2]);
      break;

    case t_QUAD:
      z=cgetg(4,t_QUAD);
      copyifstack(x[1],z[1]);
      z[2]=gcmp0(gmael(x,1,3))? lcopy((GEN) x[2])
                              : ladd((GEN) x[2],(GEN) x[3]);
      z[3]=lneg((GEN) x[3]);
      break;

    case t_POL:
      lx=lgef(x); z=cgetg(lx,tx); z[1]=x[1];
      for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
      break;

    case t_SER:
      lx=lg(x); z=cgetg(lx,tx); z[1]=x[1];
      for (i=2; i<lx; i++) z[i]=lconj((GEN) x[i]);
      break;

    case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
      lx=lg(x); z=cgetg(lx,tx);
      for (i=1; i<lx; i++) z[i]=lconj((GEN) x[i]);
      break;

    default:
      err(typeer,"gconj");
      return NULL; /* not reached */
  }
  return z;
}

GEN
conjvec(GEN x,long prec)
{
  long lx,s,av,tetpil,i,tx=typ(x);
  GEN z,y,p1,p2,p;

  switch(tx)
  {
    case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN:
      z=cgetg(2,t_COL); z[1]=lcopy(x); break;

    case t_COMPLEX: case t_QUAD:
      z=cgetg(3,t_COL); z[1]=lcopy(x); z[2]=lconj(x); break;

    case t_VEC: case t_COL:
      lx=lg(x); z=cgetg(lx,t_MAT);
      for (i=1; i<lx; i++) z[i]=(long)conjvec((GEN)x[i],prec);
      s=lg(z[1]);
      for (i=2; i<lx; i++)
	if (lg(z[i])!=s) err(talker,"incompatible field degrees in conjvec");
      break;

    case t_POLMOD:
      y=(GEN)x[1]; lx=lgef(y);
      if (lx<=3) return cgetg(1,t_COL);
      av=avma; p=NULL;
      for (i=2; i<lx; i++)
      {
	tx=typ(y[i]);
	if (tx==t_INTMOD) p=gmael(y,i,1);
	else
	  if (tx != t_INT && ! is_frac_t(tx))
	    err(polrationer,"conjvec");
      }
      if (!p)
      {
	p1=roots(y,prec); x = (GEN)x[2]; tetpil=avma;
	z=cgetg(lx-2,t_COL);
	for (i=1; i<=lx-3; i++)
	{
	  p2=(GEN)p1[i]; if (gcmp0((GEN)p2[2])) p2 = (GEN)p2[1];
	  z[i] = (long)poleval(x, p2);
	 }
	return gerepile(av,tetpil,z);
      }
      z=cgetg(lx-2,t_COL); z[1]=lcopy(x);
      for (i=2; i<=lx-3; i++) z[i] = lpui((GEN) z[i-1],p,prec);
      break;

    default:
      err(typeer,"conjvec");
      return NULL; /* not reached */
  }
  return z;
}

/*******************************************************************/
/*                                                                 */
/*                            TRACES                               */
/*                                                                 */
/*******************************************************************/

GEN
assmat(GEN x)
{
  long lx,i,j;
  GEN y,p1,p2;

  if (typ(x)!=t_POL) err(notpoler,"assmat");
  if (gcmp0(x)) err(zeropoler,"assmat");

  lx=lgef(x)-2; y=cgetg(lx,t_MAT);
  for (i=1; i<lx-1; i++)
  {
    p1=cgetg(lx,t_COL); y[i]=(long)p1;
    for (j=1; j<lx; j++)
      p1[j] = (j==(i+1))? un: zero;
  }
  p1=cgetg(lx,t_COL); y[i]=(long)p1;
  if (gcmp1((GEN) x[lx+1]))
    for (j=1; j<lx; j++)
      p1[j] = lneg((GEN)x[j+1]);
  else
  {
    p2 = (GEN)x[lx+1]; gnegz(p2,p2);
    for (j=1; j<lx; j++)
      p1[j] = ldiv((GEN)x[j+1],p2);
    gnegz(p2,p2);
  }
  return y;
}

GEN
gtrace(GEN x)
{
  long i,l,n,tx=typ(x),lx,tetpil;
  GEN y,p1,p2;

  switch(tx)
  {
    case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
      return gmul2n(x,1);

    case t_COMPLEX:
      return gmul2n((GEN)x[1],1);

    case t_QUAD:
      p1=(GEN)x[1];
      if (!gcmp0((GEN) p1[3]))
      {
	l=avma; p2=gmul2n((GEN)x[2],1);
	tetpil=avma; return gerepile(l,tetpil,gadd((GEN)x[3],p2));
      }
      return gmul2n((GEN)x[2],1);

    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
      return y;

    case t_SER:
      lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=ltrace((GEN)x[i]);
      return y;

    case t_POLMOD:
      l=avma; n=(lgef(x[1])-4);
      p1=polsym((GEN)x[1],n); p2=gzero;
      for (i=0; i<=n; i++)
        p2=gadd(p2,gmul(truecoeff((GEN)x[2],i),(GEN)p1[i+1]));
      return gerepileupto(l,p2);

    case t_RFRAC: case t_RFRACN:
      return gadd(x,gconj(x));

    case t_VEC: case t_COL:
      lx=lg(x); y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=ltrace((GEN)x[i]);
      return y;

    case t_MAT:
      lx=lg(x); if (lx!=lg(x[1])) err(mattype1,"gtrace");
      l=avma; p1=gcoeff(x,1,1); if (lx==2) return gcopy(p1);
      for (i=2; i<lx-1; i++)
	p1=gadd(p1,gcoeff(x,i,i));
      tetpil=avma; return gerepile(l,tetpil,gadd(p1,gcoeff(x,i,i)));

  }
  err(typeer,"gtrace");
  return NULL; /* not reached */
}

/* Gauss reduction for positive definite matrix a
 * If a is not positive definite:
 *   if flag is zero: raise an error
 *   else: return NULL.
 */
GEN
sqred1intern(GEN a,long flag)
{
  GEN b,p;
  long av = avma,tetpil,i,j,k, lim=stack_lim(av,1), n=lg(a);

  if (typ(a)!=t_MAT) err(typeer,"sqred1");
  if (lg(a[1])!=n) err(mattype1,"sqred1");
  b=cgetg(n,t_MAT);
  for (j=1; j<n; j++)
  {
    GEN p1=cgetg(n,t_COL), p2=(GEN)a[j];

    b[j]=(long)p1;
    for (i=1; i<=j; i++) p1[i] = p2[i];
    for (   ; i<n ; i++) p1[i] = zero;
  }
  for (k=1; k<n; k++)
  {
    p = gcoeff(b,k,k);
    if (gsigne(p)<=0) /* not positive definite */
    {
      if (flag) { avma=av; return NULL; }
      err(talker,"not a positive definite matrix in sqred1");
    }
    p = ginv(p);
    for (i=k+1; i<n; i++)
      for (j=i; j<n; j++)
	coeff(b,i,j)=lsub(gcoeff(b,i,j),
	                  gmul(gmul(gcoeff(b,k,i),gcoeff(b,k,j)), p));
    for (j=k+1; j<n; j++)
      coeff(b,k,j)=lmul(gcoeff(b,k,j), p);
    if (low_stack(lim, stack_lim(av,1)))
    {
      if (DEBUGMEM>1) err(warnmem,"sqred1");
      tetpil=avma; b=gerepile(av,tetpil,gcopy(b));
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(b));
}

GEN
sqred1(GEN a)
{
  return sqred1intern(a,0);
}

GEN
sqred3(GEN a)
{
  long av=avma,tetpil,i,j,k,l, lim=stack_lim(av,1), n=lg(a);
  GEN p1,b;

  if (typ(a)!=t_MAT) err(typeer,"sqred3");
  if (lg(a[1])!=n) err(mattype1,"sqred3");
  av=avma; b=cgetg(n,t_MAT);
  for (j=1; j<n; j++)
  {
    p1=cgetg(n,t_COL); b[j]=(long)p1;
    for (i=1; i<n; i++) p1[i]=zero;
  }
  for (i=1; i<n; i++)
  {
    for (k=1; k<i; k++)
    {
      p1=gzero;
      for (l=1; l<k; l++)
	p1=gadd(p1, gmul(gmul(gcoeff(b,l,l),gcoeff(b,k,l)), gcoeff(b,i,l)));
      coeff(b,i,k)=ldiv(gsub(gcoeff(a,i,k),p1),gcoeff(b,k,k));
    }
    p1=gzero;
    for (l=1; l<i; l++)
      p1=gadd(p1, gmul(gcoeff(b,l,l), gsqr(gcoeff(b,i,l))));
    coeff(b,i,k)=lsub(gcoeff(a,i,i),p1);
    if (low_stack(lim, stack_lim(av,1)))
    {
      if (DEBUGMEM>1) err(warnmem,"sqred3");
      tetpil=avma; b=gerepile(av,tetpil,gcopy(b));
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(b));
}

/* Gauss reduction (arbitrary symetric matrix, only the part above the
 * diagonal is considered). If no_signature is zero, return only the
 * signature
 */
static GEN
sqred2(GEN a, long no_signature)
{
  GEN r,p,mun;
  long av,tetpil,av1,lim,n,i,j,k,l,sp,sn,t;

  if (typ(a)!=t_MAT) err(typeer,"sqred2");
  n=lg(a); if (lg(a[1]) != n) err(mattype1,"sqred2");

  av=avma; mun=negi(gun); r=new_chunk(n);
  for (i=1; i<n; i++) r[i]=1;
  av1=avma; lim=stack_lim(av1,1); a=dummycopy(a);
  n--; t=n; sp=sn=0;
  while (t)
  {
    k=1; while (k<=n && (gcmp0(gcoeff(a,k,k)) || !r[k])) k++;
    if (k<=n)
    {
      p=gcoeff(a,k,k); if (gsigne(p)>0) sp++; else sn++;
      r[k]=0; t--;
      for (j=1; j<=n; j++)
	coeff(a,k,j)=r[j] ? ldiv(gcoeff(a,k,j),p) : zero;
	
      for (i=1; i<=n; i++) if (r[i])
	for (j=1; j<=n; j++)
	  coeff(a,i,j) = r[j] ? lsub(gcoeff(a,i,j),
	                             gmul(gmul(gcoeff(a,k,i),gcoeff(a,k,j)),p))
			      : zero;
      coeff(a,k,k)=(long)p;
    }
    else
    {
      for (k=1; k<=n; k++) if (r[k])
      {
	l=k+1; while (l<=n && (gcmp0(gcoeff(a,k,l)) || !r[l])) l++;
	if (l<=n)
	{
	  p=gcoeff(a,k,l); r[k]=r[l]=0; sp++; sn++; t-=2;
	  for (i=1; i<=n; i++) if (r[i])
	  {
	    for (j=1; j<=n; j++)
	      coeff(a,i,j) =
		r[j]? lsub(gcoeff(a,i,j),
			   gdiv(gadd(gmul(gcoeff(a,k,i),gcoeff(a,l,j)),
				     gmul(gcoeff(a,k,j),gcoeff(a,l,i))),
				p))
		    : zero;
	    coeff(a,k,i)=ldiv(gadd(gcoeff(a,k,i),gcoeff(a,l,i)),p);
	    coeff(a,l,i)=ldiv(gsub(gcoeff(a,k,i),gcoeff(a,l,i)),p);
	  }
	  coeff(a,k,l)=un; coeff(a,l,k)=(long)mun;
	  coeff(a,k,k)=lmul2n(p,-1); coeff(a,l,l)=lneg(gcoeff(a,k,k));
	  if (low_stack(lim, stack_lim(av1,1)))
	  {
	    if(DEBUGMEM>1) err(warnmem,"sqred2");
	    tetpil=avma; a=gerepile(av1,tetpil,gcopy(a));
	  }
	  break;
	}
      }
      if (k>n) break;
    }
  }
  if (no_signature) { tetpil=avma; return gerepile(av,tetpil,gcopy(a)); }
  avma=av;
  a=cgetg(3,t_VEC); a[1]=lstoi(sp); a[2]=lstoi(sn); return a;
}

GEN
sqred(GEN a) { return sqred2(a,1); }

GEN
signat(GEN a) { return sqred2(a,0); }

/* Diagonalization of a REAL symetric matrix. Return a 2-component vector:
 * 1st comp = vector of eigenvalues
 * 2nd comp = matrix of eigenvectors
 */
GEN
jacobi(GEN a, long prec)
{
  long de,e,e1,e2,l,n,i,j,p,q,av1,av2;
  GEN c,s,t,u,ja,lambda,r,unr,x,y,x1,y1;

  if (typ(a)!=t_MAT) err(mattype1,"jacobi");
  ja=cgetg(3,t_VEC); l=lg(a); n=l-1;
  e1=HIGHEXPOBIT-1;
  lambda=cgetg(l,t_COL); ja[1]=(long)lambda;
  for (j=1; j<=n; j++)
  {
    lambda[j]=lgetr(prec);
    gaffect(gcoeff(a,j,j), (GEN)lambda[j]);
    e=expo(lambda[j]); if (e<e1) e1=e;
  }
  r=cgetg(l,t_MAT); ja[2]=(long)r;
  for (j=1; j<=n; j++)
  {
    r[j]=lgetg(l,t_COL);
    for (i=1; i<=n; i++)
      affsr(i==j, (GEN)(coeff(r,i,j)=lgetr(prec)));
  }
  av1=avma;

  e2=-HIGHEXPOBIT;p=q=1;
  c=cgetg(l,t_MAT);
  for (j=1; j<=n; j++)
  {
    c[j]=lgetg(j,t_COL);
    for (i=1; i<j; i++)
    {
      gaffect(gcoeff(a,i,j), (GEN)(coeff(c,i,j)=lgetr(prec)));
      e=expo(gcoeff(c,i,j)); if (e>e2) { e2=e; p=i; q=j; }
    }
  }
  a=c; unr = realun(prec);
  de=bit_accuracy(prec);

 /* e1 = min des expo des coeff diagonaux
  * e2 = max des expo des coeff extra-diagonaux
  * Test d'arret: e2 < e1-precision
  */
  while (e1-e2<de)
  {
    /*calcul de la rotation associee dans le plan
	  des p et q-iemes vecteurs de base   */
    av2=avma;
    x=divrr(subrr((GEN)lambda[q],(GEN)lambda[p]),shiftr(gcoeff(a,p,q),1));
    y=mpsqrt(addrr(unr,mulrr(x,x)));
    t=(gsigne(x)>0)? divrr(unr,addrr(x,y)) : divrr(unr,subrr(x,y));
    c=divrr(unr,mpsqrt(addrr(unr,mulrr(t,t))));
    s=mulrr(t,c); u=divrr(s,addrr(unr,c));

    /* Recalcul des transformees successives de a et de la matrice
	   cumulee (r) des rotations :  */

    for (i=1; i<p; i++)
    {
      x=gcoeff(a,i,p); y=gcoeff(a,i,q);
      x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
      y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
      affrr(x1,gcoeff(a,i,p)); affrr(y1,gcoeff(a,i,q));
    }
    for (i=p+1; i<q; i++)
    {
      x=gcoeff(a,p,i); y=gcoeff(a,i,q);
      x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
      y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
      affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,i,q));
    }
    for (i=q+1; i<=n; i++)
    {
      x=gcoeff(a,p,i); y=gcoeff(a,q,i);
      x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
      y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
      affrr(x1,gcoeff(a,p,i)); affrr(y1,gcoeff(a,q,i));
    }
    x=(GEN)lambda[p]; y=gcoeff(a,p,q); subrrz(x,mulrr(t,y),(GEN)lambda[p]);
    x=y; y=(GEN)lambda[q]; addrrz(y,mulrr(t,x),y);
    setexpo(x,expo(x)-de-1);

    for (i=1; i<=n; i++)
    {
      x=gcoeff(r,i,p); y=gcoeff(r,i,q);
      x1=subrr(x,mulrr(s,addrr(y,mulrr(u,x))));
      y1=addrr(y,mulrr(s,subrr(x,mulrr(u,y))));
      affrr(x1,gcoeff(r,i,p)); affrr(y1,gcoeff(r,i,q));
    }

    e2=expo(gcoeff(a,1,2)); p=1; q=2;
    for (j=1; j<=n; j++)
    {
      for (i=1; i<j; i++)
	if ((e=expo(gcoeff(a,i,j))) > e2) { e2=e; p=i; q=j; }
      for (i=j+1; i<=n; i++)
	if ((e=expo(gcoeff(a,j,i))) > e2) { e2=e; p=j; q=i; }
    }
    avma=av2;
  }
  avma=av1; return ja;
}

/*************************************************************************/
/**									**/
/**		    MATRICE RATIONNELLE --> ENTIERE			**/
/**									**/
/*************************************************************************/

GEN
matrixqz0(GEN x,GEN p)
{
  if (typ(p)!=t_INT) err(typeer,"matrixqz0");
  if (signe(p)>=0)  return matrixqz(x,p);
  if (!cmpsi(-1,p)) return matrixqz2(x);
  if (!cmpsi(-2,p)) return matrixqz3(x);
  err(flagerr,"matrixqz"); return NULL; /* not reached */
}

GEN
matrixqz(GEN x, GEN p)
{
  long av=avma,av1,tetpil,i,j,j1,m,n,t,lim,nfact;
  GEN p1,p2,p3,unmodp;

  if (typ(x)!=t_MAT) err(typeer,"matrixqz");
  n=lg(x)-1; if (!n) return gcopy(x);
  m=lg(x[1])-1;
  if (n > m) err(talker,"more rows than columns in matrixqz");
  if (n==m)
  {
    p1=det(x);
    if (gcmp0(p1)) err(talker,"matrix of non-maximal rank in matrixqz");
    avma=av; return idmat(n);
  }
  p1=cgetg(n+1,t_MAT);
  for (j=1; j<=n; j++)
  {
    p2=gun; p3=(GEN)x[j];
    for (i=1; i<=m; i++)
    {
      t=typ(p3[i]);
      if (t != t_INT && !is_frac_t(t))
        err(talker,"not a rational or integral matrix in matrixqz");
      p2=ggcd(p2,(GEN)p3[i]);
    }
    p1[j]=ldiv(p3,p2);
  }
  x=p1; unmodp=cgetg(3,t_INTMOD); unmodp[2]=un;

  if (gcmp0(p))
  {
    p1=cgetg(n+1,t_MAT); p2=gtrans(x);
    for (j=1; j<=n; j++) p1[j]=p2[j];
    p3=det(p1); p1[n]=p2[n+1]; p3=ggcd(p3,det(p1));
    if (!signe(p3))
      err(impl,"matrixqz when the first 2 dets are zero");
    if (gcmp1(p3))
      { tetpil=avma; return gerepile(av,tetpil,gcopy(x)); }

    p3=factor(p3); p1=(GEN)p3[1]; nfact=lg(p1)-1;
  }
  else
  {
    p1=cgetg(2,t_VEC); p1[1]=(long)p; nfact=1;
  }
  av1=avma; lim=stack_lim(av1,1);
  for (i=1; i<=nfact; i++)
  {
    p=(GEN)p1[i]; unmodp[1]=(long)p;
    for(;;)
    {
      p2=ker(gmul(unmodp,x));
      if (lg(p2)==1) break;

      p2=centerlift(p2); p3=gdiv(gmul(x,p2),p);
      for (j=1; j<lg(p2); j++)
      {
	j1=n; while (gcmp0(gcoeff(p2,j1,j))) j1--;
	x[j1] = p3[j];
      }
      if (low_stack(lim, stack_lim(av1,1)))
      {
	if (DEBUGMEM>1) err(warnmem,"matrixqz");
	tetpil=avma; x=gerepile(av1,tetpil,gcopy(x));
      }
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(x));
}

static GEN
matrixqz_aux(GEN x, long m, long n)
{
  long av = avma, lim = stack_lim(av,1), tetpil,i,j,k,in[2];
  GEN p1;

  for (i=1; i<=m; i++)
  {
    for(;;)
    {
      long fl=0;

      for (j=1; j<=n; j++)
	if (!gcmp0(gcoeff(x,i,j)))
	  { in[fl]=j; fl++; if (fl==2) break; }
      if (j>n) break;

      j=(gcmp(gabs(gcoeff(x,i,in[0]),DEFAULTPREC),
	      gabs(gcoeff(x,i,in[1]),DEFAULTPREC)) > 0)? in[1]: in[0];
      p1=gcoeff(x,i,j);
      for (k=1; k<=n; k++)
	if (k!=j)
	  x[k]=lsub((GEN)x[k],gmul(ground(gdiv(gcoeff(x,i,k),p1)),(GEN)x[j]));
    }

    j=1; while (j<=n && gcmp0(gcoeff(x,i,j))) j++;
    if (j<=n)
    {
      p1=denom(gcoeff(x,i,j));
      if (!gcmp1(p1)) x[j]=lmul(p1,(GEN)x[j]);
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"matrixqz_aux");
      tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
    }
  }
  return hnf(x);
}

GEN
matrixqz2(GEN x)
{
  long av = avma,m,n;

  if (typ(x)!=t_MAT) err(typeer,"matrixqz2");
  n=lg(x)-1; if (!n) return gcopy(x);
  m=lg(x[1])-1; x=dummycopy(x);
  return gerepileupto(av, matrixqz_aux(x,m,n));
}

GEN
matrixqz3(GEN x)
{
  long av=avma,av1,j,j1,k,m,n,lim;
  GEN c;

  if (typ(x)!=t_MAT) err(typeer,"matrixqz3");
  n=lg(x)-1; if (!n) return gcopy(x);
  m=lg(x[1])-1; x=dummycopy(x); c=new_chunk(n+1);
  for (j=1; j<=n; j++) c[j]=0;
  av1=avma; lim=stack_lim(av1,1);
  for (k=1; k<=m; k++)
  {
    j=1; while (j<=n && (c[j] || gcmp0(gcoeff(x,k,j)))) j++;
    if (j<=n)
    {
      c[j]=k; x[j]=ldiv((GEN)x[j],gcoeff(x,k,j));
      for (j1=1; j1<=n; j1++)
	if (j1!=j) x[j1]=lsub((GEN)x[j1],gmul(gcoeff(x,k,j1),(GEN)x[j]));
      if (low_stack(lim, stack_lim(av1,1)))
      {
	long tetpil = avma;
	if(DEBUGMEM>1) err(warnmem,"matrixqz3");
	x=gerepile(av1,tetpil,gcopy(x));
      }
    }
  }
  return gerepileupto(av, matrixqz_aux(x,m,n));
}

GEN
intersect(GEN x, GEN y)
{
  long av,tetpil,j, lx = lg(x);
  GEN z;

  if (typ(x)!=t_MAT || typ(y)!=t_MAT) err(typeer,"intersect");
  if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);

  av=avma; z=ker(concatsp(x,y));
  for (j=lg(z)-1; j; j--) setlg(z[j],lx);
  tetpil=avma; return gerepile(av,tetpil,gmul(x,z));
}

/**************************************************************/
/**                                                          **/
/**		   HERMITE NORMAL FORM REDUCTION	     **/
/**							     **/
/**************************************************************/

GEN
mathnf0(GEN x, long flag)
{
  switch(flag)
  {
    case 0: return hnf(x);
    case 1: return hnfall(x);
    case 2: return hnfhavas(x);
    case 3: return hnfperm(x);
    case 4: return hnflll(x);
    default: err(flagerr,"mathnf");
  }
  return NULL; /* not reached */
}

static GEN
init_hnf(GEN x, GEN *denx, long *co, long *li, long *av)
{
  if (typ(x) != t_MAT) err(typeer,"mathnf");
  *co=lg(x); if (*co==1) return NULL; /* empty matrix */
  *li=lg(x[1]); *denx=denom(x); *av=avma;

  if (gcmp1(*denx)) /* no denominator */
    { *denx = NULL; return dummycopy(x); }
  return gmul(x,*denx);
}

/* return c * X */
#ifdef INLINE
INLINE
#endif
GEN
int_col_mul(GEN c, GEN X)
{
  long i,m = lg(X);
  GEN A = new_chunk(m);
  for (i=1; i<m; i++) A[i] = lmulii(c,(GEN)X[i]);
  A[0] = X[0]; return A;
}

/* X,Y columns; u,v scalars; everybody is integral. return A = u*X + v*Y */
GEN
lincomb_integral(GEN u, GEN v, GEN X, GEN Y)
{
  long av,i,lx,m;
  GEN p1,p2,A;

  if (!signe(u)) return int_col_mul(v,Y);
  if (!signe(v)) return int_col_mul(u,X);
  lx = lg(X); A = cgetg(lx,t_COL); m = lgefint(u)+lgefint(v)+4;
  if (gcmp1(u))
  {
    for (i=1; i<lx; i++)
    {
      p1=(GEN)X[i]; p2=(GEN)Y[i];
      if      (!signe(p1)) A[i] = lmulii(v,p2);
      else if (!signe(p2)) A[i] = licopy(p1);
      else
      {
        av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2));
        p2 = mulii(v,p2);
        avma = av; A[i] = laddii(p1,p2);
      }
    }
  }
  else
    for (i=1; i<lx; i++)
    {
      p1=(GEN)X[i]; p2=(GEN)Y[i];
      if      (!signe(p1)) A[i] = lmulii(v,p2);
      else if (!signe(p2)) A[i] = lmulii(u,p1);
      else
      {
        av = avma; (void)new_chunk(m+lgefint(p1)+lgefint(p2));
        p1 = mulii(u,p1);
        p2 = mulii(v,p2);
        avma = av; A[i] = laddii(p1,p2);
      }
    }
  return A;
}

GEN
hnf_special(GEN x, long remove)
{
  long av0, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
  GEN p1,u,v,d,denx,a,b, x2,res;

  if (typ(x) != t_VEC || lg(x) !=3) err(typeer,"hnf_special");
  res = cgetg(3,t_VEC);
  av0 = avma;
  x2 = (GEN)x[2];
  x  = (GEN)x[1];
  x = init_hnf(x,&denx,&co,&li,&av);
  if (!x) return cgetg(1,t_MAT);

  lim = stack_lim(av,1);
  def=co-1; ldef=(li>co)? li-co: 0;
  if (lg(x2) != co) err(talker,"incompatible matrices in hnf_special");
  x2 = dummycopy(x2);
  for (i=li-1; i>ldef; i--)
  {
    for (j=def-1; j; j--)
    {
      while (j && !signe(gcoeff(x,i,j))) j--;
      if (!j) break;
      k = (j==1)? def: j-1;
      a = gcoeff(x,i,j);
      b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
      if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
      if (DEBUGLEVEL>5) { outerr(u); outerr(v); }
      p1 = (GEN)x[j]; b = negi(b);
      x[j] = (long)lincomb_integral(a,b, (GEN)x[k], p1);
      x[k] = (long)lincomb_integral(u,v, p1, (GEN)x[k]);
      p1 = (GEN)x2[j];
      x2[j]=  ladd(gmul(a, (GEN)x2[k]), gmul(b,p1));
      x2[k] = ladd(gmul(u,p1), gmul(v, (GEN)x2[k]));
      if (low_stack(lim, stack_lim(av,1)))
      {
        GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
        if (DEBUGMEM>1) err(warnmem,"hnf_special[1]. i=%ld",i);
        gerepilemany(av,gptr,2);
      }
    }
    p1 = gcoeff(x,i,def); s = signe(p1);
    if (s)
    {
      if (s < 0)
      { 
        x[def] = lneg((GEN)x[def]); p1 = gcoeff(x,i,def);
        x2[def]= lneg((GEN)x2[def]);
      }
      for (j=def+1; j<co; j++)
      {
	b = negi(gdivent(gcoeff(x,i,j),p1));
	x[j] = (long)lincomb_integral(gun,b, (GEN)x[j], (GEN)x[def]);
        x2[j]= ladd((GEN)x2[j], gmul(b, (GEN)x2[def]));
      }
      def--;
    }
    else
      if (ldef && i==ldef+1) ldef--;
    if (low_stack(lim, stack_lim(av,1)))
    {
      GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
      if (DEBUGMEM>1) err(warnmem,"hnf_special[2]. i=%ld",i);
      gerepilemany(av,gptr,2);
    }
  }
  if (remove)
  {                            /* remove null columns */
    for (i=1,j=1; j<co; j++)
      if (!gcmp0((GEN)x[j])) 
      {
        x[i] = x[j];
        x2[i] = x2[j]; i++;
      }
    setlg(x,i);
    setlg(x2,i);
  }
  tetpil=avma;
  x = denx? gdiv(x,denx): gcopy(x);
  x2 = gcopy(x2);
  {
    GEN *gptr[2]; gptr[0]=&x; gptr[1]=&x2;
    gerepilemanysp(av0,tetpil,gptr,2);
  }
  res[1] = (long)x;
  res[2] = (long)x2;
  return res;
}

GEN
hnf0(GEN x, long remove)       /* remove: throw away lin.dep.columns, GN */
{
  long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
  GEN p1,u,v,d,denx,a,b;

  if (typ(x) == t_VEC) return hnf_special(x,remove);
  x = init_hnf(x,&denx,&co,&li,&av);
  if (!x) return cgetg(1,t_MAT);

  lim = stack_lim(av,1);
  def=co-1; ldef=(li>co)? li-co: 0;
  for (i=li-1; i>ldef; i--)
  {
    for (j=def-1; j; j--)
    {
      while (j && !signe(gcoeff(x,i,j))) j--;
      if (!j) break;
      k = (j==1)? def: j-1;
      a = gcoeff(x,i,j);
      b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
      if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
      if (DEBUGLEVEL>5) { outerr(u); outerr(v); }
      p1 = (GEN)x[j];
      x[j] = (long)lincomb_integral(a,negi(b), (GEN)x[k],p1);
      x[k] = (long)lincomb_integral(u,v, p1,(GEN)x[k]);
      if (low_stack(lim, stack_lim(av,1)))
      {
        if (DEBUGMEM>1) err(warnmem,"hnf[1]. i=%ld",i);
        tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
      }
    }
    p1 = gcoeff(x,i,def); s = signe(p1);
    if (s)
    {
      if (s < 0) { x[def] = lneg((GEN)x[def]); p1 = gcoeff(x,i,def); }
      for (j=def+1; j<co; j++)
      {
	b = negi(gdivent(gcoeff(x,i,j),p1));
	x[j] = (long)lincomb_integral(gun,b, (GEN)x[j], (GEN)x[def]);
      }
      def--;
    }
    else
      if (ldef && i==ldef+1) ldef--;
    if (low_stack(lim, stack_lim(av,1)))
    {
      if (DEBUGMEM>1) err(warnmem,"hnf[2]. i=%ld",i);
      tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
    }
  }
  if (remove)
  {                            /* remove null columns */
    for (i=1,j=1; j<co; j++)
      if (!gcmp0((GEN)x[j])) x[i++] = x[j];
    setlg(x,i);
  }
  tetpil=avma;
  x = denx? gdiv(x,denx): gcopy(x);
  return gerepile(av0,tetpil,x);
}

GEN
hnf(GEN x) { return hnf0(x,1); }

#define cmod(x,u,us2) \
  {GEN a=modii((GEN)x,u); if (cmpii(a,us2)>0) a=subii(a,u); x=(long)a;}

/* dm = multiple of diag element (usually detint(x)) */
/* flag: don't/do append dm*matid. */
static GEN
allhnfmod(GEN x,GEN dm,long flag)
{
  long li,co,av0,av,tetpil,i,j,k,def,ldef,lim,ldm;
  GEN a,b,w,p1,p2,d,u,v,dms2;

  if (typ(dm)!=t_INT) err(typeer,"allhnfmod");
  if (!signe(dm)) return hnf(x);
  if (typ(x)!=t_MAT) err(typeer,"allhnfmod");
  if (DEBUGLEVEL>6) fprintferr("Enter hnfmod");

  co=lg(x); if (co==1) return cgetg(1,t_MAT);
  av0=avma; lim=stack_lim(av0,1);
  dms2=shifti(dm,-1); li=lg(x[1]);
  av=avma;

  if (flag)
  {
    p1 = cgetg(co,t_MAT); for (i=1; i<co; i++) p1[i]=x[i];
    if (li > co) err(talker,"nb lines > nb columns in hnfmod");
    x = p1;
  }
  else
  {
    /* concatenate dm * Id to x */
    x = concatsp(x, idmat_intern(li-1,dm,gzero));
    for (i=1; i<co; i++) x[i] = lmod((GEN)x[i], dm);
    co += li-1;
  }
  def=co-1; ldef=0;
  for (i=li-1; i>ldef; i--,def--)
  {
    if (DEBUGLEVEL>6) { fprintferr(" %ld",i); flusherr(); }
    for (j=def-1; j; j--)
    {
      while (j && !signe(gcoeff(x,i,j))) j--;
      if (!j) break;
      if (DEBUGLEVEL>8) { fprintferr(" %ld",j); flusherr(); }
      k = (j==1)? def: j-1;
      a = gcoeff(x,i,j);
      b = gcoeff(x,i,k); d = bezout(a,b,&u,&v);
      if (!is_pm1(d)) { a = divii(a,d); b = divii(b,d); }
      p1 = lincomb_integral(u,v, (GEN)x[j], (GEN)x[k]);
      p2 = lincomb_integral(a, negi(b), (GEN)x[k], (GEN)x[j]);
      x[k] = (long)p1;
      x[j] = (long)p2;
      for (k=1; k<=i; k++)
      {
        cmod(p1[k], dm, dms2);
        cmod(p2[k], dm, dms2);
      }
      if (low_stack(lim, stack_lim(av0,1)))
      {
        if (DEBUGMEM>1) err(warnmem,"allhnfmod[1]. i=%ld",i);
	tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
      }
    }
  }
  w=cgetg(li,t_MAT); b=dm;
  for (i=li-1; i>=1; i--)
  {
    d = bezout(gcoeff(x,i,i+def),b,&u,&v);
    w[i] = lmod(gmul(u,(GEN)x[i+def]), b);
    if (!signe(gcoeff(w,i,i))) coeff(w,i,i)=(long)d;
    if (i > 1 && flag) b=divii(b,d);
  }
  ldm = lgefint(dm);
  for (i=li-2; i>=1; i--)
  {
    GEN diag = gcoeff(w,i,i);
    for (j=i+1; j<li; j++)
    {
      b = negi(gdivent(gcoeff(w,i,j), diag));
      p1 = lincomb_integral(gun,b, (GEN)w[j], (GEN)w[i]);
      w[j] = (long)p1;
      for (k=1; k<i; k++)
      {
        p2 = (GEN)p1[k];
        if (lgefint(p2) > ldm) p1[k] = lmodii(p2, dm);
      }
      if (low_stack(lim, stack_lim(av0,1)))
      {
        if (DEBUGMEM>1) err(warnmem,"allhnfmod[2]. i=%ld", i);
        tetpil=avma; w=gerepile(av,tetpil,gcopy(w));
        diag = gcoeff(w,i,i);
      }
    }
  }
  if (DEBUGLEVEL>6) { fprintferr("\nEnd hnfmod\n"); flusherr(); }
  tetpil=avma; return gerepile(av0,tetpil,gcopy(w));
}
#undef cmod

GEN
hnfmod(GEN x, GEN detmat) { return allhnfmod(x,detmat,1); }

GEN
hnfmodid(GEN x, GEN p) { return allhnfmod(x,p,0); }

/* Return [y,U,V] such that y=V.x.U, V permutation vector, U unimodular
 * matrix, and y in HNF form
 */
GEN
hnfhavas(GEN x)
{
  long av0=avma, av,av1,tetpil,li,co,i,j,k,def,ldef,lim,imin,jmin,vpk;
  long jpro,com,vi;
  GEN p1,p2,z,u,denx,vperm,mat1,col2,lil2,s,pmin,apro,bpro,cpro;

  if (DEBUGLEVEL>6)
    { fprintferr("Entering hnfhavas: AVMA = %ld\n",avma); flusherr(); }
  if (typ(x)!=t_MAT) err(typeer,"hnfhavas");
  co=lg(x);
  if (co==1)
  {
    z=cgetg(4,t_VEC); z[1]=lgetg(1,t_MAT);
    z[2]=lgetg(1,t_MAT); z[3]=lgetg(1,t_VEC);
    return z;
  }
  li=lg(x[1]); denx=denom(x);
  vperm=new_chunk(li); for (i=1; i<li; i++) vperm[i]=i;

  av=avma; lim=stack_lim(av,1); u=idmat(co-1);
  x = gcmp1(denx)? dummycopy(x): gmul(denx,x);
  def=co; ldef=(li>co)?li-co+1:1;
  for (i=li-1; i>=ldef; i--)
  {
    def--; av1=avma; mat1=cgetg(def+1,t_MAT); col2=cgetg(def+1,t_COL);
    for (j=1; j<=def; j++)
    {
      p1=cgetg(i+1,t_COL); mat1[j]=(long)p1; s=gzero;
      for (k=1; k<=i; k++)
      {
	p2=gsqr(gcoeff(x,vperm[k],j));
	p1[k]=(long)p2; s=gadd(s,p2);
      }
      col2[j]=(long)s;
    }
    lil2=cgetg(i+1,t_COL);
    for (k=1; k<=i; k++)
    {
      s=gzero;
      for (j=1; j<=def; j++) s=gadd(s,gcoeff(mat1,k,j));
      lil2[k]=(long)s;
    }

    pmin = NULL;
    for (k=i; k>=1; k--)
    {
      while (k>=1 && !signe(lil2[k])) k--;
      if (!k) goto comterm;
      vpk=vperm[k];
      if (!pmin || cmpii((GEN)lil2[k],pmin) <= 0)
      {
	j=1; while (!signe(gcoeff(x,vpk,j))) j++;
	if (!pmin)
	{
	  imin=k; jmin=j; pmin=mulii((GEN)lil2[k],(GEN)col2[j]);
	  cpro=absi(gcoeff(x,vpk,j));
	}
	jpro=j; apro=absi(gcoeff(x,vpk,j)); j++;
	for (; j<=def; j++)
	{
	  com=cmpii((GEN)col2[j],(GEN)col2[jpro]);
	  if (signe(gcoeff(x,vpk,j)) && com <=0)
	  {
	    if (com<0) { jpro=j; apro=absi(gcoeff(x,vpk,j)); }
	    else
	    {
	      bpro=absi(gcoeff(x,vpk,j));
	      if (cmpii(bpro,apro)<0) { jpro=j; apro=bpro; }
	    }
	  }
	}
	p1=mulii((GEN)lil2[k],(GEN)col2[jpro]);
	com=cmpii(p1,pmin);
	if (com<0 || (com==0 && cmpii(apro,cpro)<0))
	{
	  imin=k; jmin=jpro; pmin=p1; cpro=apro;
	}
      }
    }
    avma=av1;
    if (jmin!=def)
    {
      p1=(GEN)x[def]; x[def]=x[jmin]; x[jmin]=(long)p1;
      p1=(GEN)u[def]; u[def]=u[jmin]; u[jmin]=(long)p1;
    }
    if (imin!=i) { vpk=vperm[i]; vperm[i]=vperm[imin]; vperm[imin]=vpk; }
    vi=vperm[i];
    for(;;)
    {
      GEN p3,p12,p13;

      if (signe(gcoeff(x,vi,def))<0)
      {
	x[def]=lneg((GEN)x[def]); u[def]=lneg((GEN)u[def]);
      }
      p1=gcoeff(x,vi,def); p12=shifti(p1,-1); p13=negi(p12);
      for (j=1; j<def; j++)
      {
	p2=dvmdii(gcoeff(x,vi,j),p1,&p3);
	if (cmpii(p3,p13)<0) p2=addis(p2,-1);
	else { if (cmpii(p3,p12)>0) p2=addis(p2,1); }
	if (DEBUGLEVEL>5) outerr(p2);
        setsigne(p2,-signe(p2));
	x[j]=ladd((GEN)x[j],gmul(p2,(GEN)x[def]));
	u[j]=ladd((GEN)u[j],gmul(p2,(GEN)u[def]));
      }
      j=1; while (!signe(gcoeff(x,vi,j))) j++;
      if (j<def)
      {
	pmin=gnorml2((GEN)x[j]); jmin=j; apro=absi(gcoeff(x,vi,j));
	j++;
	for (; j<def; j++)
	{
	  if (signe(gcoeff(x,vi,j)))
	  {
	    p1=gnorml2((GEN)x[j]);
            com=cmpii(p1,pmin);
	    if (com<0)
	    {
	      pmin=p1; jmin=j;
	    }
	    else if (com==0)
	    {
	      bpro=absi(gcoeff(x,vi,j));
              if (cmpii(bpro,apro)<0)
	      {
	        pmin=p1; jmin=j; apro=bpro;
	      }
	    }
	  }
	}
	p1=(GEN)x[def]; x[def]=x[jmin]; x[jmin]=(long)p1;
	p1=(GEN)u[def]; u[def]=u[jmin]; u[jmin]=(long)p1;
      }
      else break;
    }
    vi=vperm[i]; p1=gcoeff(x,vi,def);
    for (j=def+1; j<co; j++)
    {
      p2=gdivent(gcoeff(x,vi,j),p1); setsigne(p2,-signe(p2));
      if (DEBUGLEVEL>5) outerr(p2);
      x[j]=ladd((GEN)x[j],gmul(p2,(GEN)x[def]));
      u[j]=ladd((GEN)u[j],gmul(p2,(GEN)u[def]));
    }

    if (low_stack(lim, stack_lim(av,1)))
    {
      GEN *gptr[2];
      if (DEBUGMEM>1) err(warnmem,"hnfhavas");
      gptr[0]=&x; gptr[1]=&u;
      gerepilemany(av,gptr,2);
    }
  }

comterm:
  tetpil=avma; z=cgetg(4,t_VEC); p1=cgetg(co,t_MAT);
  if (gcmp1(denx))
  {
    for (j=1; j<co; j++)
    {
      p2=cgetg(li,t_COL); p1[j]=(long)p2;
      for (i=1; i<li; i++)
	p2[i] = lcopy(gcoeff(x,vperm[i],j));
    }
  }
  else
  {
    for (j=1; j<co; j++)
    {
      p2=cgetg(li,t_COL); p1[j]=(long)p2;
      for (i=1; i<li; i++)
	p2[i] = ldiv(gcoeff(x,vperm[i],j),denx);
    }
  }
  z[1]=(long)p1; z[2]=lcopy(u);
  p1=cgetg(li,t_VEC); for (i=1; i<li; i++) p1[i]=lstoi(vperm[i]);
  z[3]=(long)p1; return gerepile(av0,tetpil,z);
}

/* HNF by Bo Majewski and Allan Steele */

/* premier indice non nul de la j-eme ligne de la matrice b */
static long
depthvector(GEN b,long j)
{
  long lv = lg(b), i = 1;

  while (i<lv && gcmp0(gcoeff(b,j,i))) i++;
  return (i==lv)?-1:i;
}

static GEN
incompleteprod(GEN b,long k,long l,long deb,long fin)
{
  GEN p1 = gzero;
  long j;

  for (j=deb; j<=fin; j++)
    p1 = addii(p1,mulii(gcoeff(b,k,j),gcoeff(b,l,j)));
  return p1;
}

static void
redlll(GEN b,GEN mu,long l,long c,long k)
{
  long i,j, lb;
  GEN q, p1=gcoeff(b,l,c);

  if (signe(p1)) q=gdivround(gcoeff(b,k,c),p1); else q=ground(gcoeff(mu,k,l));
  if (signe(q))
  {
    q=negi(q); lb=lg(b);
    for (j=1; j<lb; j++)
      coeff(b,k,j) = laddii(gcoeff(b,k,j),mulii(q,gcoeff(b,l,j)));
    coeff(mu,k,l)=ladd(gcoeff(mu,k,l),q);
    for (i=1; i<l; i++)
    {
      p1=gcoeff(mu,l,i);
      if (gsigne(p1))
	coeff(mu,k,i) = ladd(gcoeff(mu,k,i),gmul(q,p1));
    }
  }
}

GEN
hnflll(GEN x)
{
  long n,m,i,j,k,ii,jj,p,c,s,av=avma,tetpil,kmax,ok;
  GEN q,qneg,p1,bnew,B,mu,mmu,cst,E,U,y,t,BB;

  if (typ(x)!=t_MAT) err(typeer,"hnflll");
  n=lg(x)-1;
  if (!n)
  {
    y=cgetg(3,t_VEC);
    y[1]=lgetg(1,t_MAT);
    y[2]=lgetg(1,t_MAT); return y;
  }
  cst=gdivgs(stoi(9),10);
  x=gcopy(x); n=lg(x)-1; m=lg(x[1])-1; p=n+m;
  bnew=cgetg(p+1,t_MAT);
  for (j=1; j<=n; j++) bnew[j]=x[j];
  for (j=n+1; j<=p; j++)
  {
    p1=cgetg(m+1,t_COL); bnew[j]=(long)p1;
    for (i=1; i<=m; i++) p1[i]=(i==(j-n))?un:zero;
  }
  x=bnew; c=n+1;
  for (i=1; i<=m; i++) c=min(c,depthvector(x,i));
  s=0; mu=cgetg(m+2,t_MAT);
  for (j=1; j<=m; j++) mu[j]=lgetg(m+2,t_COL); B=cgetg(m+2,t_COL);
  while (c<=n)
  {
    k=2; kmax=1;
    B[1]=(long)incompleteprod(x,1,1,c+1,p);
    while (k<=m-c+1)
    {
      if (k>kmax)
      {
	kmax=k;
	for (j=1; j<k; j++)
	{
	  mmu=incompleteprod(x,k,j,c+1,p);
	  for (i=1; i<j; i++) mmu=gsub(mmu,gmul(gcoeff(mu,j,i),gcoeff(mu,k,i)));
	  coeff(mu,k,j)=(long)mmu;
	}
	for (j=1; j<k; j++) coeff(mu,k,j)=ldiv(gcoeff(mu,k,j),(GEN)B[j]);
	B[k]=(long)incompleteprod(x,k,k,c+1,p);
	for (j=1; j<k; j++)
	  B[k]=lsub((GEN)B[k],gmul((GEN)B[j],gsqr(gcoeff(mu,k,j))));
      }
      redlll(x,mu,k-1,c,k);
      ok = (absi(gcoeff(x,k-1,c))>absi(gcoeff(x,k,c))) ||
	  (gegal(gcoeff(x,k-1,c),gcoeff(x,k,c)) &&
	    (gcmp((GEN) B[k],
	          gmul(gsub(cst,gsqr(gcoeff(mu,k,k-1))), (GEN) B[k-1])) < 0) );
      while (ok)
      {
	for (j=1; j<=p; j++)
	{
	  t=gcoeff(x,k,j); coeff(x,k,j)=coeff(x,k-1,j);
	  coeff(x,k-1,j)=(long)t;
	}
	for (j=1; j<=k-2; j++)
	{
	  t=gcoeff(mu,k,j); coeff(mu,k,j)=coeff(mu,k-1,j);
	  coeff(mu,k-1,j)=(long)t;
	}
	mmu=gcoeff(mu,k,k-1);
	BB=gadd((GEN)B[k],gmul(gmul(mmu,mmu),(GEN)B[k-1]));
	q=gdiv((GEN)B[k-1],BB);
	coeff(mu,k,k-1)=lmul(mmu,q);
	B[k]=lmul((GEN)B[k],q); B[k-1]=(long)BB;
	for (i=k+1; i<=kmax; i++)
	{
	  t=gcoeff(mu,i,k);
	  coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mmu,t));
	  coeff(mu,i,k-1)=ladd(t,gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
	}
	if (k>2) k--;
	redlll(x,mu,k-1,c,k);
        ok=(absi(gcoeff(x,k-1,c))>absi(gcoeff(x,k,c))) ||
	   (gegal(gcoeff(x,k-1,c),gcoeff(x,k,c)) &&
	     (gcmp((GEN) B[k],
	           gmul(gsub(cst,gsqr(gcoeff(mu,k,k-1))), (GEN) B[k-1])) < 0));
      }
      for (i=k-2; i; i--) redlll(x,mu,i,c,k);
      k++;
    }
    s++; c=n+1;
    for (i=1; i<=m-s; i++) c=min(c,depthvector(x,i));
  }
  s=m-s+1;
  if (signe(gcoeff(x,s,depthvector(x,s)))<0)
    for (j=1; j<=p; j++)
      coeff(x,s,j)=lnegi(gcoeff(x,s,j));
  for (i=s+1; i<=m; i++)
  {
    if (signe(gcoeff(x,i,depthvector(x,i)))<0)
      for (j=1; j<=p; j++)
	coeff(x,i,j)=lnegi(gcoeff(x,i,j));
    for (j=i-1; j>=s; j--)
    {
      k=depthvector(x,j);
      qneg=negi(gdivent(gcoeff(x,i,k),gcoeff(x,j,k)));
      for (jj=1; jj<=p; jj++)
	coeff(x,i,jj)=laddii(gcoeff(x,i,jj),mulii(qneg,gcoeff(x,j,jj)));
    }
  }
  for (k=s; k<=m; k++)
  {
    for (j=1; j<s; j++)
    {
      mmu=incompleteprod(x,k,j,n+1,p);
      for (i=1; i<j; i++) mmu=gsub(mmu,gmul(gcoeff(mu,j,i),gcoeff(mu,k,i)));
      coeff(mu,k,j)=(long)mmu;
    }
    for (j=1; j<s; j++) coeff(mu,k,j)=ldiv(gcoeff(mu,k,j),(GEN)B[j]);
    B[k]=(long)incompleteprod(x,k,k,n+1,p);
    for (j=1; j<s; j++)
      B[k]=lsub((GEN)B[k],gmul(gsqr(gcoeff(mu,k,j)),(GEN)B[j]));
    for (j=s-1; j; j--)
    {
      qneg=negi(ground(gcoeff(mu,k,j)));
      if (signe(qneg))
      {
	for (jj=1; jj<=p; jj++)
	  coeff(x,k,jj)=laddii(gcoeff(x,k,jj),mulii(qneg,gcoeff(x,j,jj)));
	for (i=1; i<j; i++)
	  if (gsigne(gcoeff(mu,j,i)))
	    coeff(mu,k,i)=ladd(gcoeff(mu,k,i),gmul(qneg,gcoeff(mu,j,i)));
      }
    }
  }
  tetpil=avma; y=cgetg(3,t_VEC);
  E=cgetg(n+1,t_MAT);
  for (i=1; i<=n; i++)
  {
    p1=cgetg(m-s+2,t_COL); E[i]=(long)p1;
    for (ii=1; ii<=m-s+1; ii++)
      p1[ii]=lcopy(gcoeff(x,m-ii+1,i));
  }
  y[1]=(long)E; U=cgetg(m+1,t_MAT);
  for (i=1; i<=m; i++)
  {
    p1=cgetg(m+1,t_COL); U[i]=(long)p1;
    for (ii=m; ii>=1; ii--)
      p1[m-ii+1]=lcopy(gcoeff(x,ii,i+n));
  }
  y[2]=(long)U; return gerepile(av,tetpil,y);
}

/* HNF with permutation */
GEN
hnfperm(GEN A)
{
  GEN U,c,l,perm,d,u,v,p,q,y,a,b,p1;
  long r,t,i,j,j1,k,m,n,av=avma,av1,tetpil,lim;

  if (typ(A)!=t_MAT) err(typeer,"hnfperm");
  n=lg(A)-1;
  if (!n)
  {
    y=cgetg(4,t_VEC);
    y[1]=lgetg(1,t_MAT);
    y[2]=lgetg(1,t_MAT);
    y[3]=lgetg(1,t_VEC); return y;
  }
  m=lg(A[1])-1;
  c=new_chunk(m+1); for (i=1; i<=m; i++) c[i]=0;
  l=new_chunk(n+1); for (j=1; j<=n; j++) l[j]=0;
  perm=new_chunk(m+1);
  av1=avma; lim=stack_lim(av1,1);
  U=idmat(n); A=dummycopy(A);
/* U base change matrix : A0*U=A all along */

  for (r=0,k=1; k<=n; k++)
  {
    for (j=1; j<k; j++) if (l[j])
    {
      t=l[j]; b=gcoeff(A,t,k);
      if (signe(b) == 0) continue;

      a=gcoeff(A,t,j); d=bezout(a,b,&u,&v);
      if (!is_pm1(d)) { a=divii(a,d); b=divii(b,d); }
      p1 = (GEN)A[j]; b = negi(b);
      A[j] = (long)lincomb_integral(u,v, p1,(GEN)A[k]);
      A[k] = (long)lincomb_integral(a,b, (GEN)A[k],p1);
      p1 = (GEN)U[j];
      U[j] = (long)lincomb_integral(u,v, p1,(GEN)U[k]);
      U[k] = (long)lincomb_integral(a,b, (GEN)U[k],p1);
      for (j1=1; j1<j; j1++) if (l[j1])
      {
        q=truedvmdii(gcoeff(A,t,j1),d,NULL);
        if (signe(q))
        {
          q = negi(q);
          A[j1] = (long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[j]);
          U[j1] = (long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[j]);
        }
      }
    }
    t=m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
    if (t)
    {
      p = gcoeff(A,t,k);
      for (i=t-1; i; i--)
      {
        q = gcoeff(A,i,k);
	if (signe(q) && absi_cmp(p,q) > 0) { p=q; t=i; }
      }
      perm[++r]=l[k]=t; c[t]=k;
      if (signe(p) < 0)
      {
	for (i=1; i<=m; i++) coeff(A,i,k)= lnegi(gcoeff(A,i,k));
	for (i=1; i<=n; i++) coeff(U,i,k)= lnegi(gcoeff(U,i,k));
	p=gcoeff(A,t,k);
      }
      for (j=1; j<k; j++) if (l[j])
      {
	q=truedvmdii(gcoeff(A,t,j),p,NULL);
	if (signe(q))
	{
          q = negi(q);
          A[j]=(long)lincomb_integral(gun,q,(GEN)A[j],(GEN)A[k]);
          U[j]=(long)lincomb_integral(gun,q,(GEN)U[j],(GEN)U[k]);
	}
      }
    }
    if (low_stack(lim, stack_lim(av1,1)))
    {
      GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
      if (DEBUGMEM>1) err(warnmem,"hnfperm");
      gerepilemany(av1,gptr,2);
    }
  }
  if (r < m)
  {
    for (i=1,k=r; i<=m; i++)
      if (c[i]==0) perm[++k] = i;
  }

/* We have A0*U=A, U in Gl(n,Z)
 * basis for Im(A):  columns of A s.t l[j]>0 (r   cols)
 * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols)
 */
  tetpil=avma; y=cgetg(4,t_VEC);
  p=cgetg(r+1,t_MAT); u=cgetg(n+1,t_MAT);
  for (t=1,k=r,j=1; j<=n; j++)
    if (l[j])
    {
      q=cgetg(m+1,t_COL); p[k]=(long)q;
      for (i=1; i<=m; i++) q[i]=lcopy(gcoeff(A,perm[m-i+1],j));
      u[k+n-r]=lcopy((GEN)U[j]);
      k--;
    }
    else u[t++]=lcopy((GEN)U[j]);
  y[1]=(long)p; y[2]=(long)u;
  q = cgetg(m+1,t_VEC); y[3]=(long)q;
  for (i=1; i<=m; i++) q[m-i+1]=lstoi(perm[i]);
  return gerepile(av,tetpil,y);
}

GEN
colint_copy(GEN x)
{
  long i, lx = lg(x);
  GEN y = cgetg(lx, t_COL);
  for (i=1; i<lx; i++) y[i] = signe(x[i])? licopy((GEN)x[i]): zero;
  return y;
}

GEN
matint_copy(GEN x)
{
  long i, lx = lg(x);
  GEN y = cgetg(lx, t_MAT);

  for (i=1; i<lx; i++)
    y[i] = (long)colint_copy((GEN)x[i]);
  return y;
}

/*====================================================================
 *	    Forme Normale d'Hermite (Version par colonnes 31/01/94)
 *====================================================================*/
GEN
hnfall0(GEN A, long remove)
{
  GEN U,c,h,x,y,u,v,p,q,d,a,b,p1;
  long m,n,r,i,j,j1,k,li,z,av=avma,av1,tetpil,lim;

  if (typ(A)!=t_MAT) err(typeer,"hnfall");
  n=lg(A)-1;
  if (!n)
  {
    y=cgetg(3,t_VEC);
    y[1]=lgetg(1,t_MAT);
    y[2]=lgetg(1,t_MAT); return y;
  }
  m=lg(A[1])-1;
  c=new_chunk(m+1); for (i=1; i<=m; i++) c[i]=0;
  h=new_chunk(n+1); for (j=1; j<=n; j++) h[j]=m;
  av1=avma; lim=stack_lim(av1,1);
  A=dummycopy(A); U=idmat(n); r=n+1;
  for (li=m; li; li--)
  {
    for (j=1; j<r; j++)
    {
      for (i=h[j]; i>li; i--)
      {
        b = gcoeff(A,i,j);
	if (signe(b))
	{
	  k=c[i]; a=gcoeff(A,i,k); /* annuler bij A l'aide de p=bik */
	  d=bezout(a,b,&u,&v);
	  if (DEBUGLEVEL>5)
            { fprintferr("(u,v) = (%Z, %Z); ",u,v); flusherr(); }
	  if (!is_pm1(d)) { a=divii(a,d); b =divii(b,d); }
          p1 = (GEN)A[k]; b = negi(b);
          A[k] = (long)lincomb_integral(u,v, p1,(GEN)A[j]);
          A[j] = (long)lincomb_integral(a,b, (GEN)A[j],p1);
          p1 = (GEN)U[k];
          U[k] = (long)lincomb_integral(u,v, p1,(GEN)U[j]);
          U[j] = (long)lincomb_integral(a,b, (GEN)U[j],p1);
	  for (j1=k+1; j1<=n; j1++)
	  {
	    q=truedvmdii(gcoeff(A,i,j1),d,NULL);
	    if (signe(q))
	    {
              q = negi(q);
              A[j1]=(long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[k]);
              U[j1]=(long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[k]);
	    }
	  }
          if (low_stack(lim, stack_lim(av1,1)))
          {
            GEN *gptr[2]; tetpil = avma;
            A = matint_copy(A); gptr[0]=&A;
            U = matint_copy(U); gptr[1]=&U;
            if (DEBUGMEM>1) err(warnmem,"hnfall[1], li = %ld", li);
            gerepilemanysp(av1,tetpil,gptr,2);
          }	
	}
      }
      x=gcoeff(A,li,j);
      if (signe(x))
      {
        r--;
        if (j<r)
        {
          z=A[j]; A[j]=A[r]; A[r]=z;
          z=U[j]; U[j]=U[r]; U[r]=z;
          h[j]=h[r]; h[r]=li; c[li]=r;
        }
        if (signe(gcoeff(A,li,r))<0)
        {
          p1=(GEN)A[r]; for (i=1; i<=li; i++) p1[i]=lnegi((GEN)p1[i]);
          p1=(GEN)U[r]; for (i=1; i<=n ; i++) p1[i]=lnegi((GEN)p1[i]);
        }
        p=gcoeff(A,li,r);
        for (j=r+1; j<=n; j++)
        {
          q = truedvmdii(gcoeff(A,li,j),p,NULL);
          if (signe(q))
          {
            q = negi(q);
            A[j]=(long)lincomb_integral(gun,q,(GEN)A[j],(GEN)A[r]);
            U[j]=(long)lincomb_integral(gun,q,(GEN)U[j],(GEN)U[r]);
          }
        }
        if (low_stack(lim, stack_lim(av1,1)))
        {
          GEN *gptr[2]; gptr[0]=&A; gptr[1]=&U;
          if (DEBUGMEM>1) err(warnmem,"hnfall[2], li = %ld", li);
          gerepilemany(av1,gptr,2);
        }	
        break;
      }
      h[j]=li-1;
    }
  }
  if (DEBUGLEVEL>5) fprintferr("\nhnfall, final phase: ");
  r--; /* first r cols are in the image the n-r (independent) last ones */
  for (j=1; j<=r; j++)
    for (i=h[j]; i; i--)
      if (signe(b=gcoeff(A,i,j)))
      {
	k=c[i]; a=gcoeff(A,i,k);
	d=bezout(a,b,&u,&v);
	if (!is_pm1(d)) { a=divii(a,d); b=divii(b,d); }
        if (DEBUGLEVEL>5)
          { fprintferr("(u,v) = (%Z, %Z); ",u,v); flusherr(); }
        p1 = (GEN)A[k]; b = negi(b);
        A[k] = (long)lincomb_integral(u,v, p1,(GEN)A[j]);
        A[j] = (long)lincomb_integral(a,b, (GEN)A[j],p1);
        p1 = (GEN)U[k];
        U[k] = (long)lincomb_integral(u,v, p1,(GEN)U[j]);
        U[j] = (long)lincomb_integral(a,b, (GEN)U[j],p1);
	for (j1=k+1; j1<=n; j1++)
	{
	  q=truedvmdii(gcoeff(A,i,j1),d,NULL);
	  if (signe(q))
	  {
            q = negi(q);
            A[j1] = (long)lincomb_integral(gun,q,(GEN)A[j1],(GEN)A[k]);
            U[j1] = (long)lincomb_integral(gun,q,(GEN)U[j1],(GEN)U[k]);
	  }
	}
        if (low_stack(lim, stack_lim(av1,1)))
        {
          GEN *gptr[2]; tetpil = avma;
          A = matint_copy(A); gptr[0]=&A;
          U = matint_copy(U); gptr[1]=&U;
          if (DEBUGMEM>1) err(warnmem,"hnfall[3], j = %ld", j);
          gerepilemanysp(av1,tetpil,gptr,2);
        }	
      }
  if (DEBUGLEVEL>5) fprintferr("\n");
  tetpil=avma; y=cgetg(3,t_VEC);
  if (remove)
  { /* remove the first r columns */
    A += r; A[0] = evaltyp(t_MAT) | evallg(n-r+1);
  }
  y[1]=lcopy(A); y[2]=lcopy(U);
  return gerepile(av,tetpil,y);
}

GEN
hnfall(GEN x) {return hnfall0(x,1);}

/***************************************************************/
/**							      **/
/**      	    SMITH NORMAL FORM REDUCTION	              **/
/**							      **/
/***************************************************************/

static GEN
col_mul(GEN x, GEN c)
{
  long s = signe(x);
  GEN xc = NULL;
  if (s)
  {
    if (!is_pm1(x)) xc = gmul(x,c);
    else xc = (s>0)? c: gneg_i(c);
  }
  return xc;
}

static void
do_zero(GEN x)
{
  long i, lx = lg(x);
  for (i=1; i<lx; i++) x[i] = zero;
}

/* c1 <-- u.c1 + v.c2; c2 <-- a.c2 - b.c1 */
static void
update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
{
  GEN p1,p2;

  u = col_mul(u,*c1);
  v = col_mul(v,*c2);
  if (u) p1 = v? gadd(u,v): u;
  else   p1 = v? v: (GEN)NULL;

  a = col_mul(a,*c2);
  b = col_mul(gneg_i(b),*c1);
  if (a) p2 = b? gadd(a,b): a;
  else   p2 = b? b: (GEN)NULL;

  if (!p1) do_zero(*c1); else *c1 = p1;
  if (!p2) do_zero(*c2); else *c2 = p2;
}

static GEN
trivsmith(long all)
{
  GEN z;
  if (!all) return cgetg(1,t_VEC);
  z=cgetg(4,t_VEC);
  z[1]=lgetg(1,t_MAT);
  z[2]=lgetg(1,t_MAT);
  z[3]=lgetg(1,t_VEC); return z;
}

/* Return the smith normal form d of matrix x. If all != 0 return [d,u,v],
 * where d = u.x.v
 */
static GEN
smithall(GEN x, long all)
{
  long av,tetpil,i,j,k,l,c,fl,n,s1,s2,lim;
  GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr,mun,mdet,ys;

  if (typ(x)!=t_MAT) err(typeer,"smithall");
  if (DEBUGLEVEL>=9) outerr(x);
  n=lg(x)-1;
  if (!n) return trivsmith(all);
  if (lg(x[1]) != n+1) err(mattype1,"smithall");
  for (i=1; i<=n; i++)
    for (j=1; j<=n; j++)
      if (typ(coeff(x,i,j)) != t_INT)
        err(talker,"non integral matrix in smithall");

  av = avma; lim = stack_lim(av,1);
  x = dummycopy(x); mdet = detint(x);
  if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }
  else
  {
    if (signe(mdet))
    {
      p1=hnfmod(x,mdet);
      if (all) { ml=idmat(n); mr=gauss(x,p1); }
    }
    else
    {
      if (all)
      {
        p1 = hnfall0(x,0);
        ml = idmat(n); mr = (GEN)p1[2]; p1 = (GEN)p1[1];
      }
      else
        p1 = hnf0(x,0);
    }
    x = p1;
  }
  p1=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) p1[i]=lnegi(gcoeff(x,i,i));
  p2=sindexsort(p1); ys=cgetg(n+1,t_MAT);
  for (j=1; j<=n; j++)
  {
    p1=cgetg(n+1,t_COL); ys[j]=(long)p1;
    for (i=1; i<=n; i++) p1[i]=coeff(x,p2[i],p2[j]);
  }
  x = ys;
  if (all)
  {
    p3=cgetg(n+1,t_MAT); p4=cgetg(n+1,t_MAT);
    for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }
    ml=p3; mr=p4;
  }
  if (signe(mdet))
  {
    p1 = hnfmod(x,mdet);
    if (all) mr=gmul(mr,gauss(x,p1));
  }
  else
  {
    if (all)
    {
      p1 = hnfall0(x,0);
      mr = gmul(mr,(GEN)p1[2]); p1 = (GEN)p1[1];
    }
    else
      p1 = hnf0(x,0);
  }
  x=p1; mun = negi(gun);

  if (DEBUGLEVEL>7) {fprintferr("starting SNF loop");flusherr();}
  for (i=n; i>=2; i--)
  {
    if (DEBUGLEVEL>7) {fprintferr("\ni = %ld: ",i);flusherr();}
    for(;;)
    {
      c = 0;
      for (j=i-1; j>=1; j--)
      {
	p1=gcoeff(x,i,j); s1 = signe(p1);
	if (s1)
	{
	  p2=gcoeff(x,i,i);
          if (!absi_cmp(p1,p2))
          {
            s2=signe(p2);
            if (s1 == s2) { d=p1; u=gun; p4=gun; }
            else
	    {
              if (s2>0) { u = gun; p4 = mun; }
              else      { u = mun; p4 = gun; }
	      d=(s1>0)? p1: absi(p1);
	    }
            v = gzero; p3 = u;
          }
          else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
	  for (k=1; k<=i; k++)
	  {
	    b=addii(mulii(u,gcoeff(x,k,i)),mulii(v,gcoeff(x,k,j)));
	    coeff(x,k,j)=lsubii(mulii(p3,gcoeff(x,k,j)),
	                        mulii(p4,gcoeff(x,k,i)));
	    coeff(x,k,i)=(long)b;
	  }
	  if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
          if (low_stack(lim, stack_lim(av,1)))
	  {
	    if (DEBUGMEM>1) err(warnmem,"[1]: smithall");
	    if (all)
	    {
	      GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
	      gerepilemany(av,gptr,3);
	    }
	    else x=gerepileupto(av,gcopy(x));
	  }
	}
      }
      if (DEBUGLEVEL>=8) {fprintferr("; ");flusherr();}
      for (j=i-1; j>=1; j--)
      {
	p1=gcoeff(x,j,i); s1 = signe(p1);
	if (s1)
	{
	  p2=gcoeff(x,i,i);
	  if (!absi_cmp(p1,p2))
          {
            s2 = signe(p2);
            if (s1 == s2) { d=p1; u=gun; p4=gun; }
            else
	    {
              if (s2>0) { u = gun; p4 = mun; }
              else      { u = mun; p4 = gun; }
	      d=(s1>0)? p1: absi(p1);
	    }
            v = gzero; p3 = u;
          }
          else { d=bezout(p2,p1,&u,&v); p3=divii(p2,d); p4=divii(p1,d); }
	  for (k=1; k<=i; k++)
	  {
	    b=addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
	    coeff(x,j,k)=lsubii(mulii(p3,gcoeff(x,j,k)),
	                        mulii(p4,gcoeff(x,i,k)));
	    coeff(x,i,k)=(long)b;
	  }
	  if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
	  c = 1;
	}
      }
      if (!c)
      {
	b=gcoeff(x,i,i); fl=1;
	if (signe(b))
	{
	  for (k=1; k<i && fl; k++)
	    for (l=1; l<i && fl; l++)
	      fl = (int)!signe(resii(gcoeff(x,k,l),b));
          /* cast to (int) necessary for gcc-2.95 on sparcv9-64 (IS) */
	  if (!fl)
	  {
	    k--;
	    for (l=1; l<=i; l++)
	      coeff(x,i,l)=laddii(gcoeff(x,i,l),gcoeff(x,k,l));
	    if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
	  }
	}
        if (fl) break;
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
	if (DEBUGMEM>1) err(warnmem,"[2]: smithall");
	if (all)
	{
	  GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
	  gerepilemany(av,gptr,3);
	}
	else x=gerepileupto(av,gcopy(x));
      }
    }
  }
  if (DEBUGLEVEL>7) {fprintferr("\n");flusherr();}
  if (all)
  {
    for (k=1; k<=n; k++)
      if (signe(gcoeff(x,k,k))<0)
        { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lnegi(gcoeff(x,k,k)); }
    tetpil=avma; z=cgetg(4,t_VEC);
    z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
    return gerepile(av,tetpil,z);
  }
  tetpil=avma; z=cgetg(n+1,t_VEC); j=n;
  for (k=n; k; k--)
    if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));
  for (   ; j; j--) z[j]=zero;
  return gerepile(av,tetpil,z);
}

GEN
smith(GEN x) { return smithall(x,0); }

GEN
smith2(GEN x) { return smithall(x,1); }

/* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
GEN
smithclean(GEN z)
{
  long i,j,l,c;
  GEN u,v,y,d,p1;

  if (typ(z) != t_VEC) err(typeer,"smithclean");
  l = lg(z); if (l == 1) return cgetg(1,t_VEC);
  u=(GEN)z[1];
  if (l != 4 || typ(u) != t_MAT) 
  {
    if (typ(u) != t_INT) err(typeer,"smithclean");
    for (c=1; c<l; c++)
      if (gcmp1((GEN)z[c])) break;
    return gcopy_i(z, c);
  }
  v=(GEN)z[2]; d=(GEN)z[3]; l = lg(d);
  for (c=1; c<l; c++)
    if (gcmp1(gcoeff(d,c,c))) break;

  y=cgetg(4,t_VEC);
  y[1]=(long)(p1 = cgetg(l,t_MAT));
  for (i=1; i<l; i++) p1[i] = (long)gcopy_i((GEN)u[i], c);
  y[2]=(long)gcopy_i(v, c);
  y[3]=(long)(p1 = cgetg(c,t_MAT));
  for (i=1; i<c; i++)
  {
    GEN p2 = cgetg(c,t_COL); p1[i] = (long)p2;
    for (j=1; j<c; j++) p2[j] = i==j? lcopy(gcoeff(d,i,i)): zero;
  }
  return y;
}

static GEN
gsmithall(GEN x,long all)
{
  long av,tetpil,i,j,k,l,c,fl,n, lim;
  GEN p1,p2,p3,p4,z,b,u,v,d,ml,mr;

  if (typ(x)!=t_MAT) err(typeer,"gsmithall");
  n=lg(x)-1;
  if (!n) return trivsmith(all);
  if (lg(x[1]) != n+1) err(mattype1,"gsmithall");
  av = avma; lim = stack_lim(av,1);
  x = dummycopy(x);
  if (all) { ml=idmat(n); mr=idmat(n); }
  for (i=n; i>=2; i--)
  {
    do
    {
      c=0;
      for (j=i-1; j>=1; j--)
      {
	p1=gcoeff(x,i,j);
	if (signe(p1))
	{
	  p2=gcoeff(x,i,i); v=gdiventres(p1,p2);
	  if (gcmp0((GEN)v[2])) { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
	  else { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
	  for (k=1; k<=i; k++)
	  {
	    b=gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
	    coeff(x,k,j)=lsub(gmul(p3,gcoeff(x,k,j)),gmul(p4,gcoeff(x,k,i)));
	    coeff(x,k,i)=(long)b;
	  }
	  if (all) update(u,v,p3,p4,(GEN*)(mr+i),(GEN*)(mr+j));
	}
      }
      for (j=i-1; j>=1; j--)
      {
	p1=gcoeff(x,j,i);
	if (signe(p1))
	{
	  p2=gcoeff(x,i,i); v=gdiventres(p1,p2);
	  if (gcmp0((GEN)v[2])) { d=p2; p4=(GEN)v[1]; v=gzero; p3=gun; u=gun; }
	  else { d=gbezout(p2,p1,&u,&v); p3=gdiv(p2,d); p4=gdiv(p1,d); }
	  for (k=1; k<=i; k++)
	  {
	    b=gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
	    coeff(x,j,k)=lsub(gmul(p3,gcoeff(x,j,k)),gmul(p4,gcoeff(x,i,k)));
	    coeff(x,i,k)=(long)b;
	  }
	  if (all) update(u,v,p3,p4,(GEN*)(ml+i),(GEN*)(ml+j));
          c = 1;
	}
      }
      if (!c)
      {
	b=gcoeff(x,i,i); fl=1;
	if (signe(b))
	{
	  for (k=1; (k<i)&&fl; k++)
	    for (l=1; (l<i)&&fl; l++)
	      fl= !signe(gmod(gcoeff(x,k,l),b));
	  if (!fl)
	  {
	    k--;
	    for (l=1; l<=i; l++)
	      coeff(x,i,l)=ladd(gcoeff(x,i,l),gcoeff(x,k,l));
	    if (all) ml[i]=ladd((GEN)ml[i],(GEN)ml[k]);
	  }
	}
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
	if (DEBUGMEM>1) err(warnmem,"[5]: smithall");
	tetpil=avma;
	if (all)
	{
	  GEN *gptr[3]; gptr[0]=&x; gptr[1]=&ml; gptr[2]=&mr;
	  gerepilemany(av,gptr,3);
	}
	else x=gerepile(av,tetpil,gcopy(x));
      }
    }
    while (c || !fl);
  }
  if (all)
  {
    for (k=1; k<=n; k++)
      if (signe(gcoeff(x,k,k))<0)
      { mr[k]=lneg((GEN)mr[k]); coeff(x,k,k)=lneg(gcoeff(x,k,k)); }
    tetpil=avma; z=cgetg(4,t_VEC);
    z[1]=ltrans(ml); z[2]=lcopy(mr); z[3]=lcopy(x);
  }
  else
  {
    tetpil=avma; z=cgetg(n+1,t_VEC);
    for (j=0,k=1; k<=n; k++) if (!signe(gcoeff(x,k,k))) z[++j]=zero;
    for (k=1; k<=n; k++)
      if (signe(p1=gcoeff(x,k,k))) z[++j]=(long)gabs(p1,0);
  }
  return gerepile(av,tetpil,z);
}

GEN
matsnf0(GEN x,long flag)
{
  long av = avma;
  if (flag > 7) err(flagerr,"matsnf");
  if (typ(x) == t_VEC && flag & 4) return smithclean(x);
  if (flag & 2) x = gsmithall(x,flag & 1);
  else          x = smithall(x, flag & 1);
  if (flag & 4) x = smithclean(x);
  return gerepileupto(av, x);
}

GEN
gsmith(GEN x) { return gsmithall(x,0); }

GEN
gsmith2(GEN x) { return gsmithall(x,1); }