/********************************************************************/
/**                                                                **/
/**                 LLL Algorithm and close friends                **/
/**                                                                **/
/********************************************************************/
/* $Id: bibli1.c,v 1.4 1999/09/24 10:51:24 karim Exp $ */
#include "pari.h"
#include "parinf.h"
GEN lincomb_integral(GEN u, GEN v, GEN X, GEN Y);

/* scalar product of x with himself */
static GEN
sqscal(GEN x)
{
  long i,av=avma,lx=lg(x);
  GEN z=gzero;

  for (i=1; i<lx; i++)
    z = gadd(z, gsqr((GEN)x[i]));
  return gerepileupto(av,z);
}

/* scalar product x . y */
static GEN
gscal(GEN x,GEN y)
{
  long i,av=avma,lx=lg(x);
  GEN z=gzero;

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

static GEN
sqscali(GEN x)
{
  GEN z = gzero;
  long i,lx=lg(x);

  for (i=1; i<lx; i++)
    z = addii(z,sqri((GEN)x[i]));
  return z;
}

static GEN
gscali(GEN x,GEN y)
{
  GEN z = gzero;
  long i,lx=lg(x);

  for (i=1; i<lx; i++)
    z = addii(z, mulii((GEN)x[i],(GEN)y[i]));
  return z;
}

static GEN
lllall_trivial(GEN x, long n, long flag)
{
  GEN y;
  if (!n)
  {
    if (flag != lll_ALL) return cgetg(1,t_MAT);
    y=cgetg(3,t_VEC);
    y[1]=lgetg(1,t_MAT);
    y[2]=lgetg(1,t_MAT); return y;
  }
  /* here n = 1 */
  if (gcmp0((GEN)x[1]))
  {
    switch(flag)
    {
      case lll_KER: return idmat(1);
      case lll_IM : return cgetg(1,t_MAT);
      default: y=cgetg(3,t_VEC);
        y[1]=(long)idmat(1);
        y[2]=lgetg(1,t_MAT); return y;
    }
  }
  if (flag & lll_GRAM) flag ^= lll_GRAM; else x = NULL;
  switch (flag)
  {
    case lll_KER: return cgetg(1,t_MAT);
    case lll_IM : return idmat(1);
    default: y=cgetg(3,t_VEC);
      y[1]=lgetg(1,t_MAT);
      y[2]=x? lcopy(x): (long)idmat(1); return y;
  }
}

static GEN
lllgramall_finish(GEN h,GEN fl,long flag,long n)
{
  long k;
  GEN y;

  k=1; while (k<=n && !fl[k]) k++;
  switch(flag)
  {
    case lll_KER: setlg(h,k);
      y = gcopy(h); break;

    case lll_IM: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
      y = gcopy(h); break;

    default: setlg(h,k); y=cgetg(3,t_VEC);
      y[1] = lcopy(h);
      h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
      y[2] = lcopy(h);
      break;
  }
  return y;
}

/********************************************************************/
/**                                                                **/
/**                          LLL with content                      **/
/**                                                                **/
/********************************************************************/

/* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */
static GEN
lllgramintwithcontent(GEN x, GEN veccon, long flag)
{
  long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
  GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon;
  GEN *gptr[8];

  if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent");
  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent");
  fl = new_chunk(lx);

  av=avma; lim=stack_lim(av,1);
  x=dummycopy(x); veccon=dummycopy(veccon);
  B=cgetg(lx+1,t_COL); B[1]=un;
  /* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */

  for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; }
  lam=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++)
  { p2=cgetg(lx,t_COL); lam[j]=(long)p2; for (i=1; i<lx; i++) p2[i]=zero; }
/* le vrai lam[i,j] est
   lam[i,j]*veccon[i]*veccon[j]*prod(1,l=1,j-1,veccon[l]^2) */
  k=2; h=idmat(n); kmax=1;
  u=gcoeff(x,1,1); if (typ(u)!=t_INT) err(lllger4);
  if (signe(u)) { B[2]=(long)u; coeff(lam,1,1)=un; fl[1]=1; }
  else { B[2]=un; fl[1]=0; }
  if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); }
  for(;;)
  {
    if (k>kmax)
    {
      kmax=k;
      for (j=1; j<=k; j++)
      {
	if (j==k || fl[j])
	{
	  u=gcoeff(x,k,j); if (typ(u)!=t_INT) err(lllger4);
	  for (i=1; i<j; i++)
	    if (fl[i])
	      u=divii(subii(mulii((GEN)B[i+1],u),mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[i]);
	  if (j<k) coeff(lam,k,j)=(long)u;
	  else
	  {
	    if (signe(u)) { B[k+1]=(long)u; coeff(lam,k,k)=un; fl[k]=1; }
	    else { B[k+1]=B[k]; fl[k]=0; }
	  }
	}
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
	if(DEBUGMEM>1) err(warnmem,"[1]: lllgramintwithcontent");
	gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
	gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
      }
    }
    u=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1);
    u2=mulii((GEN)B[k],(GEN)veccon[k-1]);
    if (cmpii(absi(u),u2)>0)
    {
      q=dvmdii(addii(u,u2),shifti(u2,1),&r);
      if (signe(r)<0) q=addsi(-1,q);
      h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[k-1]));
      newcon=mppgcd((GEN)veccon[k],(GEN)veccon[k-1]);
      corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
      if(!gcmp1(corr))
      {
	corr2=sqri(corr);
	for (j=1; j<=n; j++)
	  coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
	coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
	for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
	for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
	for (i=k+1; i<=kmax; i++)
	{
	  coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
	  for (j=k+1; j<i; j++)
	    coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
	}
      }
      r=negi(mulii(q,divii((GEN)veccon[k-1],(GEN)veccon[k])));
      p1=gcoeff(x,k,k-1);
      for (j=1; j<=n; j++)
	coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,k-1)));
      coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,k-1,k)));
      coeff(lam,k,k-1)=laddii(gcoeff(lam,k,k-1),mulii(r,(GEN)B[k]));
      for (i=1; i<k-1; i++)
	coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,k-1,i)));
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"[2]: lllgramintwithcontent");
      gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
      gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
    }
    p3=mulii((GEN)B[k-1],(GEN)B[k+1]);la=gcoeff(lam,k,k-1);p4=mulii(la,la);
    p2=mulsi(100,mulii(mulii((GEN)veccon[k],(GEN)veccon[k]),addii(p3,p4)));
    p1=mulii((GEN)veccon[k-1],(GEN)B[k]);p1=mulsi(99,mulii(p1,p1));
    if (fl[k-1] && (cmpii(p1,p2)>0 || !fl[k]))
    {
      if (DEBUGLEVEL>=4 && k==n)
	{ fprintferr(" (%ld)", expi(p1)-expi(p2)); flusherr(); }
      p1=(GEN)h[k-1]; h[k-1]=h[k]; h[k]=(long)p1;
      p1=(GEN)x[k-1]; x[k-1]=x[k]; x[k]=(long)p1;
      for (j=1; j<=n; j++)
      { p1=gcoeff(x,k-1,j); coeff(x,k-1,j)=coeff(x,k,j); coeff(x,k,j)=(long)p1; }
      p1=(GEN)veccon[k-1]; veccon[k-1]=veccon[k]; veccon[k]=(long)p1;
      for (j=1; j<=k-2; j++)
      { p1=gcoeff(lam,k-1,j); coeff(lam,k-1,j)=coeff(lam,k,j); coeff(lam,k,j)=(long)p1; }
      if (fl[k])
      {
	for (i=k+1; i<=kmax; i++)
	{
	  bb=gcoeff(lam,i,k);
	  coeff(lam,i,k)=ldivii(subii(mulii((GEN)B[k+1],gcoeff(lam,i,k-1)),mulii(la,bb)),(GEN)B[k]);
	  coeff(lam,i,k-1)=ldivii(addii(mulii(la,gcoeff(lam,i,k-1)),mulii((GEN)B[k-1],bb)),(GEN)B[k]);
          if (low_stack(lim, stack_lim(av,1)))
	  {
	    if(DEBUGMEM>1) err(warnmem,"[3]: lllgramintwithcontent");
	    gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
	    gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p3;
	    gptr[7]=&p4; gerepilemany(av,gptr,8);
	  }
	}
	B[k]=ldivii(addii(p3,p4),(GEN)B[k]);
      }
      else
      {
	if (signe(la))
	{
	  p2=(GEN)B[k]; p1=divii(p4,p2);
	  for (i=k+1; i<=kmax; i++)
	    coeff(lam,i,k-1)=ldivii(mulii(la,gcoeff(lam,i,k-1)),p2);
	  for (j=k+1; j<kmax; j++)
	  {
	    for (i=j+1; i<=kmax; i++)
	      coeff(lam,i,j)=ldivii(mulii(p1,gcoeff(lam,i,j)),p2);
            if (low_stack(lim, stack_lim(av,1)))
	    {
	      if(DEBUGMEM>1) err(warnmem,"[4]: lllgramintwithcontent");
	      gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
	      gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p1;
	      gptr[7]=&p2; gerepilemany(av,gptr,8);
	    }
	  }
	  B[k+1]=B[k]=(long)p1;
	  for (i=k+2; i<=lx; i++)
	    B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
	}
	else
	{
	  coeff(lam,k,k-1)=zero;
	  for (i=k+1; i<=kmax; i++)
	  {
	    coeff(lam,i,k)=coeff(lam,i,k-1);
	    coeff(lam,i,k-1)=zero;
	  }
	  B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
	}

        if (low_stack(lim, stack_lim(av,1)))
	{
	  if(DEBUGMEM>1) err(warnmem,"[5]: lllgramintwithcontent");
	  gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
	  gptr[3]=&x; gptr[4]=&veccon;
	  gerepilemany(av,gptr,5);
	}
      }
      if (k>2) k--;
      if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
    }
    else
    {
      for (l=k-2; l>=1; l--)
      {
	u=shifti(mulii(gcoeff(lam,k,l),(GEN)veccon[k]),1);
	u2=mulii((GEN)B[l+1],(GEN)veccon[l]);
	if (cmpii(absi(u),u2)>0)
	{
	  q=dvmdii(addii(u,u2),shifti(u2,1),&r);
	  if (signe(r)<0) q=addsi(-1,q);
	  h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
	  newcon=mppgcd((GEN)veccon[k],(GEN)veccon[l]);
	  corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
	  if(!gcmp1(corr))
	  {
	    corr2=sqri(corr);
	    for (j=1; j<=n; j++)
	      coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
	    coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
	    for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
	    for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
	    for (i=k+1; i<=kmax; i++)
	    {
	      coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
	      for (j=k+1; j<i; j++)
		coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
	    }
	  }
	  r=negi(mulii(q,divii((GEN)veccon[l],(GEN)veccon[k])));
	  p1=gcoeff(x,k,l);
	  for (j=1; j<=n; j++)
	    coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,l)));
	  coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,l,k)));
	  coeff(lam,k,l)=laddii(gcoeff(lam,k,l),mulii(r,(GEN)B[l+1]));
	  for (i=1; i<l; i++)
	    coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,l,i)));
	}
        if (low_stack(lim, stack_lim(av,1)))
	{
	  if(DEBUGMEM>1) err(warnmem,"[6]: lllgramintwithcontent");
	  gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
	  gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
	}
      }
      k++; if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
      if (k>n)
      {
        if (DEBUGLEVEL>5) { fprintferr("\n"); flusherr(); }
	tetpil=avma;
	return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
      }
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"[7]: lllgramintwithcontent");
      gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
      gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
    }
  }
}

static GEN
lllintwithcontent(GEN x)
{
  long lx=lg(x),i,j,av,tetpil;
  GEN veccon,con,xred,g;

  if (typ(x) != t_MAT) err(typeer,"lllintwithcontent");
  if (lx==1) return cgetg(1,t_MAT);
  av=avma; veccon=cgetg(lx,t_VEC); g=cgetg(lx,t_MAT); xred=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++)
  {
    g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]);
    xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con;
  }
  for (i=1; i<lx; i++)
    for (j=1; j<=i; j++)
      coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)xred[i],(GEN)xred[j]);
  tetpil=avma;
  return gerepile(av,tetpil,lllgramintwithcontent(g,veccon,2));
}

/********************************************************************/
/**                                                                **/
/**                          LLL ALGORITHM                         **/
/**                                                                **/
/********************************************************************/
#define swap(x,y) { long _t=x; x=y; y=_t; }
#define gswap(x,y) { GEN _t=x; x=y; y=_t; }

static void
lllupdate_int(GEN x, GEN h, GEN L, GEN B, long K, long k, long l)
{
  long i,lx;
  GEN r = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
  GEN *hk,*hl,*xk,*xl;
  if (!signe(r)) return;
  r = negi(r); lx = lg(x);
  hk = (GEN*)h[k]; hl = (GEN*)h[l];
  xk = (GEN*)x[k]; xl = (GEN*)x[l];
  if (is_pm1(r))
  {
    if (signe(r) > 0)
    {
      for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
      for (i=1;i<lx;i++) xk[i]=addii(xk[i],xl[i]);
      for (i=1;i<lx;i++) coeff(x,k,i)=laddii(gcoeff(x,k,i),gcoeff(x,l,i));
      for (i=1;i<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i));
    } else {
      for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
      for (i=1;i<lx;i++) xk[i]=subii(xk[i],xl[i]);
      for (i=1;i<lx;i++) coeff(x,k,i)=lsubii(gcoeff(x,k,i),gcoeff(x,l,i));
      for (i=1;i<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
    }
  } else {
    for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(r,hl[i]));
    for(i=1;i<lx;i++) xk[i]=addii(xk[i],mulii(r,xl[i]));
    for(i=1;i<lx;i++) coeff(x,k,i)=laddii(gcoeff(x,k,i),mulii(r,gcoeff(x,l,i)));
    for(i=1;i<l;i++)  coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
  }
  coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,B));
}

static void
lllupdate(GEN x, GEN h, GEN L, long K, long k, long l)
{
  long e,i,lx;
  GEN r = grndtoi(gcoeff(L,k,l),&e);
  GEN *hk,*hl,*xk,*xl;
  if (DEBUGLEVEL>8)
    fprintferr("error bits when rounding in lllgram: %ld\n",e);
  if (!signe(r)) return;
  r = negi(r); lx = lg(x);
  hk = (GEN*)h[k]; hl = (GEN*)h[l];
  xk = (GEN*)x[k]; xl = (GEN*)x[l];
  if (is_pm1(r))
  {
    if (signe(r) > 0)
    {
      for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
      for (i=1;i<lx;i++) xk[i]=gadd(xk[i],xl[i]);
      for (i=1;i<lx;i++) coeff(x,k,i)=ladd(gcoeff(x,k,i),gcoeff(x,l,i));
      for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gcoeff(L,l,i));
    } else {
      for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
      for (i=1;i<lx;i++) xk[i]=gsub(xk[i],xl[i]);
      for (i=1;i<lx;i++) coeff(x,k,i)=lsub(gcoeff(x,k,i),gcoeff(x,l,i));
      for (i=1;i<l; i++) coeff(L,k,i)=lsub(gcoeff(L,k,i),gcoeff(L,l,i));
    }
  } else {
    for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(r,hl[i]));
    for (i=1;i<lx;i++) xk[i]=gadd(xk[i],gmul(r,xl[i]));
    for (i=1;i<lx;i++) coeff(x,k,i)=ladd(gcoeff(x,k,i),gmul(r,gcoeff(x,l,i)));
    for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gmul(r,gcoeff(L,l,i)));
  }
  coeff(L,k,l)=ladd(gcoeff(L,k,l),r);
}

/* x integer matrix */
GEN
lllgramall(GEN x, long alpha, long flag)
{
  long av0=avma,av,tetpil,lim,lx=lg(x),i,j,k,l,n,kmax;
  GEN u,B,L,h,la,p1,p2,p3,p4,fl, *gptr[6];

  if (typ(x) != t_MAT) err(typeer,"lllgramall");
  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall");
  fl = new_chunk(lx);

  av=avma; lim=stack_lim(av,1); x=dummycopy(x);
  B=gscalcol(gun, lx);
  L=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++)
  {
    for (i=1; i<lx; i++)
      if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
    fl[j] = 0; L[j] = (long)zerocol(n);
  }
  k=2; h=idmat(n); kmax=1;
  u=gcoeff(x,1,1);
  if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
  if (DEBUGLEVEL>5) fprintferr("k =");
  for(;;)
  {
    if (k > kmax)
    {
      if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
      kmax = k;
      for (j=1; j<=k; j++)
	if (j==k || fl[j])
	{
          long av1 = avma;
	  u = gcoeff(x,k,j);
	  for (i=1; i<j; i++) if (fl[i])
            u = divii(subii(mulii((GEN)B[i+1],u),
                            mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
                      (GEN)B[i]);
          u = gerepileuptoint(av1, u);
	  if (j<k) coeff(L,k,j)=(long)u;
	  else
	  {
	    if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
	    else { B[k+1]=B[k]; fl[k]=0; }
	  }
	}
    } else if (DEBUGLEVEL>5) {fprintferr(" %ld",k); flusherr();}
    lllupdate_int(x,h,L,(GEN)B[k],kmax,k,k-1);
    if (fl[k-1] &&
     (cmpii(mulsi(alpha-1,sqri((GEN)B[k])),
            mulsi(alpha,p3=addii(mulii((GEN)B[k-1],(GEN)B[k+1]),
                               p4=sqri(la=gcoeff(L,k,k-1))))) > 0
      || !fl[k]))
    {
      if (DEBUGLEVEL>3 && k==n)
      {
        fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri((GEN)B[k])))
                           - expi(mulsi(alpha,p3)));
        flusherr();
      }
      swap(h[k-1], h[k]);
      swap(x[k-1], x[k]);
      for (j=1; j<=n; j++)  swap(coeff(x,k-1,j), coeff(x,k,j));
      for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
      if (fl[k])
      {
        GEN Bk = (GEN)B[k];
        long av1 = avma;
	for (i=k+1; i<=kmax; i++)
	{
	  GEN bb=gcoeff(L,i,k);
          p1 = divii(subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)),
                           mulii(la,bb)), Bk);
	  av1=avma=coeff(L,i,k) = (long)icopy_av(p1,(GEN)av1);
          p1 = divii(addii(mulii(la,gcoeff(L,i,k-1)),
                           mulii((GEN)B[k-1],bb)), Bk);
	  av1=avma=coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av1);
	}
	B[k] = ldivii(p3,Bk);
      }
      else
      {
	if (signe(la))
	{
	  p2=(GEN)B[k]; p1=divii(p4,p2);
	  B[k+1]=B[k]=(long)p1;
	  for (i=k+2; i<=lx; i++)
	    B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
	  for (i=k+1; i<=kmax; i++)
	    coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
	  for (j=k+1; j<kmax; j++)
	  {
	    for (i=j+1; i<=kmax; i++)
	      coeff(L,i,j)=ldivii(mulii(p1,gcoeff(L,i,j)),p2);
	    if (low_stack(lim, stack_lim(av,1)))
	    {
	      if(DEBUGMEM>1) err(warnmem,"lllgramall[1]");
	      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
	      gptr[3]=&x; gptr[4]=&p2;
	      gerepilemany(av,gptr,5); p1=(GEN)B[k];
	    }
	  }
	}
	else
	{
	  for (i=k+1; i<=kmax; i++)
	  { coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
	  B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
	}
      }
      if (k>2) k--;
    }
    else
    {
      for (l=k-2; l; l--)
      {
        lllupdate_int(x,h,L,(GEN)B[l+1],kmax,k,l);
        if (low_stack(lim, stack_lim(av,1)))
        {
          if(DEBUGMEM>1) err(warnmem,"lllgramall[2]");
          gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
          gptr[3]=&x; gerepilemany(av,gptr,4);
        }
      }
      if (++k > n) break;
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"lllgramall[3]");
      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
      gptr[3]=&x; gerepilemany(av,gptr,4);
    }
  }
  if (DEBUGLEVEL>3) fprintferr("\n");
  tetpil=avma; return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
}

static GEN
lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); }

static int
pslg(GEN x)
{
  long tx;
  if (gcmp0(x)) return 2;
  tx=typ(x); return is_scalar_t(tx)? 3: lgef(x);
}

static GEN
lllgramallgen(GEN x, long flag)
{
  long av0=avma,av,tetpil,lx=lg(x),tu,i,j,k,l,n,lim;
  GEN u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl;
  int ps1,ps2,flc;

  if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
  n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen");

  fl = new_chunk(lx);

  av=avma; lim=stack_lim(av,1);
  B=cgetg(lx+1,t_COL);
  B[1]=un; lam=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL);
  for (i=1; i<lx; i++)
    for (j=1; j<=i; j++)
    {
      if (j<i && !fl[j]) coeff(lam,i,j)=coeff(lam,j,i)=zero;
      else
      {
	u=gcoeff(x,i,j); tu=typ(u);
	if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
	for (k=1; k<j; k++)
	  if (fl[k])
	    u=gdiv(gsub( gmul((GEN)B[k+1],u),
                         gmul(gcoeff(lam,i,k),gcoeff(lam,j,k)) ),
		   (GEN)B[k]);
	if (j<i) { coeff(lam,i,j)=(long)u; coeff(lam,j,i)=zero; }
	else
	{
	  if (!gcmp0(u)) { B[i+1]=(long)u; coeff(lam,i,i)=un; fl[i]=1; }
	  else { B[i+1]=B[i]; coeff(lam,i,i)=zero; fl[i]=0; }
	}
      }
    }
  k=2; h=idmat(n); flc=0;
  for(;;)
  {
    u=gcoeff(lam,k,k-1);
    if (pslg(u) >= pslg((GEN)B[k]))
    {
      q=gdeuc(u,(GEN)B[k]); cq=gdivsg(1,content(q)); q=gmul(q,cq); flc=1;
      h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1]));
      coeff(lam,k,k-1)=lsub(gmul(cq,gcoeff(lam,k,k-1)),gmul(q,(GEN)B[k]));
      for (i=1; i<k-1; i++)
	coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,k-1,i)));
    }
    ps1 = pslg(gsqr((GEN)B[k]));
    p3 = gmul((GEN)B[k-1],(GEN)B[k+1]);
    la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1));
    ps2=pslg(gadd(p3,p4));
    if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k]))
    {
      flc = (ps1!=ps2);
      swap(h[k-1],h[k]);
      for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j));
      if (fl[k])
      {
	for (i=k+1; i<=n; i++)
	{
	  bb=gcoeff(lam,i,k);
	  coeff(lam,i,k)=ldiv(gsub(gmul((GEN)B[k+1],gcoeff(lam,i,k-1)),gmul(la,bb)),(GEN)B[k]);
	  coeff(lam,i,k-1)=ldiv(gadd(gmul(la,gcoeff(lam,i,k-1)),gmul((GEN)B[k-1],bb)),(GEN)B[k]);
	 }
	B[k]=ldiv(gadd(p3,p4),(GEN)B[k]);
      }
      else
      {
	if (!gcmp0(la))
	{
	  p2=(GEN)B[k]; p1=gdiv(p4,p2);
	  for (i=k+1; i<lx; i++)
	    coeff(lam,i,k-1)=ldiv(gmul(la,gcoeff(lam,i,k-1)),p2);
	  for (j=k+1; j<lx-1; j++)
	    for (i=j+1; i<lx; i++)
	      coeff(lam,i,j)=ldiv(gmul(p1,gcoeff(lam,i,j)),p2);
	  B[k+1]=B[k]=(long)p1;
	  for (i=k+2; i<=lx; i++)
	    B[i]=ldiv(gmul(p1,(GEN)B[i]),p2);
	 }
	else
	{
	  coeff(lam,k,k-1)=zero;
	  for (i=k+1; i<lx; i++)
	  { coeff(lam,i,k)=coeff(lam,i,k-1); coeff(lam,i,k-1)=zero; }
	  B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
	 }
      }
      if (k>2) k--;
    }
    else
    {
      for (l=k-2; l>=1; l--)
      {
	u=gcoeff(lam,k,l);
	if (pslg(u)>=pslg((GEN)B[l+1]))
	{
	  q=gdeuc(u,(GEN)B[l+1]); cq=gdivsg(1,content(q));
          q=gmul(q,cq); flc=1;
	  h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[l]));
	  coeff(lam,k,l)=lsub(gmul(cq,gcoeff(lam,k,l)),gmul(q,(GEN)B[l+1]));
	  for (i=1; i<l; i++)
            coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,l,i)));
	}
      }
      if (++k > n) break;
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      GEN *gptr[4];
      if(DEBUGMEM>1) err(warnmem,"lllgramallgen");
      gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
      gerepilemany(av,gptr,3);
    }
  }
  tetpil=avma;
  return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
}

/* return x[k,j] - (mu.A)[j] */
#ifdef INLINE
INLINE
#endif
GEN
get_Aj(GEN x, GEN mu, GEN A, long j, long k)
{
  long av,i;
  GEN s;

  if (j==1) return gcopy(gcoeff(x,k,1));
  av = avma; s = gmul(gcoeff(mu,j,1),(GEN)A[1]);
  for (i=2; i<j; i++) s=gadd(s,gmul(gcoeff(mu,j,i),(GEN)A[i]));
  s=gneg(s); return gerepileupto(av, gadd(gcoeff(x,k,j),s));
}

/* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
 * Quality ratio = (alpha-1)/alpha. Suggested value: alpha = 100
 */
GEN
lllgramintern(GEN x, long alpha, long flag, long prec)
{
  GEN xinit,L,h,A,B,L1,L2,q,cst;
  long retry = 2, av = avma,tetpil,lim,l,i,j,k,k1,lx=lg(x),n,kmax;

  if (typ(x) != t_MAT) err(typeer,"lllgram");
  n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram");
  if (n<=1) return idmat(n);
  lim = stack_lim(av,1); xinit=x;
  for (k=2,j=1; j<lx; j++)
  {
    GEN p1=(GEN)x[j];
    for (i=1; i<lx; i++) /* FIXME: lg <-> expo */
      if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; }
  }
  if (k == 2)
  {
    if (!prec) return lllgramint(x);
    x = gmul(x, realun(prec+1));
  }
  else if (prec < k) prec = k;
  h = idmat(n); 
  x = gprec_w(x, prec+1);

LABLLLGRAM:
  switch(retry--)
  {
    case 2: /* entry */ break;
    case 1: /* failed already */
      tetpil = avma; h = gcopy(h);
      prec = (prec<<1)-2;
      if (DEBUGLEVEL > 3) fprintferr("\n");
      if (DEBUGLEVEL) err(warnprec,"lllgramintern",prec);
      x = qf_base_change(gprec_w(xinit,prec),h,1);
      {
        GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x;
        gerepilemanysp(av, tetpil, gsav, 2);
      }
      if (DEBUGLEVEL) err(warner,"lllgramintern starting over");
      break;
    case 0: /* give up */
      if (DEBUGLEVEL > 3) fprintferr("\n");
      if (DEBUGLEVEL) err(warner,"lllgramintern giving up");
      if (flag) { avma=av; return NULL; }
      if (DEBUGLEVEL) outerr(xinit);
      err(lllger3);
  }
  cst = cgetr(prec+1); affsr(alpha-1,cst);
  cst = divrs(cst,alpha);

  L=cgetg(lx,t_MAT);
  B=cgetg(lx,t_COL);
  A=cgetg(lx,t_VEC);
  for (j=1; j<lx; j++)
  {
    L[j] = (long)zerocol(n);
    A[j] = B[j] = zero;
  }
  k=2; kmax=1; B[1]=coeff(x,1,1);
  if (gcmp0((GEN)B[1]))
  {
    if (flag) return NULL;
    err(lllger3);
  }
  if (DEBUGLEVEL>5) fprintferr("k =");
  for(;;)
  {
    if (k>kmax)
    {
      if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
      kmax=k;
      for (j=1; j<k; j++)
      {
        A[j] = (long)get_Aj(x,L,A,j,k);
        coeff(L,k,j) = ldiv((GEN)A[j],(GEN)B[j]);
      }
      B[k] = (long)get_Aj(x,L,A,k,k);
      if (gsigne((GEN)B[k]) <= 0) 
      {
        if (kmax == 2) retry = 0;
        goto LABLLLGRAM;
      }
    }
    else if (DEBUGLEVEL>5) fprintferr(" %ld",k);
    L1=gcoeff(L,k,k-1);
    if (DEBUGLEVEL>9)
    {
      fprintferr(" %ld", gexpo(L1) - bit_accuracy(lg(L1)));
      for (i=1; i<lx; i++)
        fprintferr("%ld: %Z\n",i,qfeval(x,(GEN)h[i]));
      flusherr();
    }
    if (2*gexpo(L1) > bit_accuracy(lg(L1)))
    {
      if (DEBUGLEVEL>3)
      {
	fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n",kmax);
        if (DEBUGLEVEL>9) fprintferr("Old B = %Z\n",B);
      }
      for (k1=1; k1<=kmax; k1++)
      {
	for (j=1; j<k1; j++)
	{
	  A[j] = (long)get_Aj(x,L,A,j,k1);
          coeff(L,k1,j) = ldiv((GEN)A[j],(GEN)B[j]);
	}
        B[k1] = (long)get_Aj(x,L,A,k1,k1);
	if (gsigne((GEN)B[k1]) <= 0) goto LABLLLGRAM;
      }
      if (DEBUGLEVEL>9) fprintferr("New B = %Z\n",B);
    }
    lllupdate(x,h,L,kmax,k,k-1);
    L1 = gcoeff(L,k,k-1);
    L2 = gsqr(L1);
    q = gmul((GEN)B[k-1], gsub(cst,L2));
    if (gcmp(q,(GEN)B[k]) > 0)
    {
      GEN BK,BB;
      if (DEBUGLEVEL>3 && k==kmax)
      { fprintferr(" (%ld)",gexpo(q)-gexpo((GEN)B[k])); flusherr(); }
      BB = gadd((GEN)B[k], gmul((GEN)B[k-1],L2));
      if (gcmp0(BB)) goto LABLLLGRAM;
      coeff(L,k,k-1) = ldiv(gmul(L1,(GEN)B[k-1]), BB);
      BK = gdiv((GEN)B[k],BB);
      B[k] = lmul((GEN)B[k-1], BK);
      B[k-1] = (long)BB;
      swap(h[k-1],h[k]);
      swap(x[k-1],x[k]);
      for (j=1; j<=n; j++)  swap(coeff(x,k-1,j), coeff(x,k,j));
      for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
      for (i=k+1; i<=kmax; i++)
      {
	GEN p=gcoeff(L,i,k);
        coeff(L,i,k) = lsub(gcoeff(L,i,k-1),gmul(L1,p));
	coeff(L,i,k-1)=ladd(gmul(BK,p), gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)));
      }
      if (k>2) k--;
    }
    else
    {
      for (l=k-2; l; l--) lllupdate(x,h,L,kmax,k,l);
      if (++k > n) break;
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      GEN *gptr[6];
      if(DEBUGMEM>1) 
      {
        if (DEBUGLEVEL > 3) fprintferr("\n");
        err(warnmem,"lllgram");
      }
      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&A;
      gptr[4]=&x; gptr[5]=&cst; gerepilemany(av,gptr,6);
    }
  }
  if (DEBUGLEVEL>3) fprintferr("\n");
  tetpil=avma; return gerepile(av,tetpil,gcopy(h));
}

static GEN
lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); }

GEN
lllgram(GEN x,long prec) { return lllgramintern(x,100,0,prec); }

GEN
qflll0(GEN x, long flag, long prec)
{
  switch(flag)
  {
    case 0: return lll(x,prec);
    case 1: return lllint(x);
    case 2: return lllintpartial(x);
    case 3: return lllrat(x);
    case 4: return lllkerim(x);
    case 5: return lllkerimgen(x);
    case 7: return lll1(x,prec);
    case 8: return lllgen(x);
    case 9: return lllintwithcontent(x);
    default: err(flagerr,"qflll");
  }
  return NULL; /* not reached */
}

GEN
qflllgram0(GEN x, long flag, long prec)
{
  switch(flag)
  {
    case 0: return lllgram(x,prec);
    case 1: return lllgramint(x);
    case 4: return lllgramkerim(x);
    case 5: return lllgramkerimgen(x);
    case 7: return lllgram1(x,prec);
    case 8: return lllgramgen(x);
    default: err(flagerr,"qflllgram");
  }
  return NULL; /* not reached */
}

/* x est la matrice d'une base b_i; retourne la matrice u (entiere
 * unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base
 * reduite est c=b*u).
 */
static GEN
lll_proto(GEN x, GEN f(GEN, long), long prec)
{
  long lx=lg(x),i,j,av,av1;
  GEN g;

  if (typ(x) != t_MAT) err(typeer,"lll_proto");
  if (lx==1) return cgetg(1,t_MAT);
  av=avma; g=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
  for (i=1; i<lx; i++)
    for (j=1; j<=i; j++)
      coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
  av1=avma; x = f(g,prec);
  if (!x) { avma=av; return NULL; }
  return gerepile(av,av1,x);
}

GEN
lllintern(GEN x,long flag,long prec)
{
  return lll_proto(x,flag? lllgram_noerr: lllgram,prec);
}

GEN
lll(GEN x,long prec) { return lll_proto(x,lllgram,prec); }

GEN
lll1(GEN x,long prec) { return lll_proto(x,lllgram1,prec); }

GEN
lllrat(GEN x) { return lll_proto(x,lllgram,lll_ALL); }

GEN
lllint(GEN x) { return lll_proto(x,lllgramall0,lll_IM); }

GEN
lllgen(GEN x) { return lll_proto(x,lllgramallgen,lll_IM); }

GEN
lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); }

GEN
lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); }

static GEN
lllkerim_proto(GEN x, GEN f(GEN,long))
{
  long lx=lg(x), i,j,av,av1;
  GEN g;

  if (typ(x) != t_MAT) err(typeer,"lllkerim_proto");
  if (lx==1)
  {
    g=cgetg(3,t_VEC);
    g[1]=lgetg(1,t_MAT);
    g[2]=lgetg(1,t_MAT); return g;
  }
  if (lg((GEN)x[1])==1)
  {
    g=cgetg(3,t_VEC);
    g[1]=(long)idmat(lx-1);
    g[2]=lgetg(1,t_MAT); return g;
  }
  av=avma; g=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
  for (i=1; i<lx; i++)
    for (j=1; j<=i; j++)
      coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
  av1=avma; return gerepile(av,av1,f(g,lll_ALL));
}

GEN
lllkerim(GEN x) { return lllkerim_proto(x,lllgramall0); }

GEN
lllkerimgen(GEN x) { return lllkerim_proto(x,lllgramallgen); }

/* x est ici la matrice de GRAM des bi.  */
GEN
lllgram1(GEN x, long prec)
{
  GEN mu,u,B,BB,BK,p,q,r,cst,unreel,sv,mu1,mu2;
  long av,tetpil,lim,l,i,j,k,lx=lg(x),n,e;

  if (typ(x) != t_MAT) err(typeer,"lllgram1");
  if (lg((GEN)x[1])!=lx) err(mattype1,"lllgram1"); n=lx-1;
  if (n<=1) return idmat(n);
  cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
  if (prec)
  {
    unreel = realun(prec+1);
    x = gmul(x,unreel);
    cst = gmul(cst,unreel);
  }
  av=avma; lim=stack_lim(av,1);
  mu=gtrans(sqred(x)); B=cgetg(lx,t_COL);
  for (i=1,l=0; i<=n; i++)
  {
    if (gsigne((GEN)(B[i]=coeff(mu,i,i)))>0) l++;
    coeff(mu,i,i)=un;
  }
  if (l<n) err(lllger3);

  u=idmat(n); k=2;
  do
  {
    if (!gcmp0(r=grndtoi(gcoeff(mu,k,k-1),&e)))
    {
      u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[k-1]));
      for (j=1; j<k-1; j++)
	coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,k-1,j)));
      mu1=(GEN)(coeff(mu,k,k-1)=lsub(gcoeff(mu,k,k-1),r));
    }
    else mu1=gcoeff(mu,k,k-1);
    q=gmul((GEN)B[k-1],gsub(cst,mu2=gsqr(mu1)));
    if (gcmp(q,(GEN)B[k])>0)
    {
      BB=gadd((GEN)B[k],gmul((GEN)B[k-1],mu2));
      coeff(mu,k,k-1)=ldiv(gmul(mu1,(GEN)B[k-1]),BB);
      B[k]=lmul((GEN)B[k-1],BK=gdiv((GEN)B[k],BB));
      B[k-1]=(long)BB;
      swap(u[k-1],u[k]);
      for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j), coeff(mu,k,j));
      for (i=k+1; i<=n; i++)
      {
	p=gcoeff(mu,i,k);
	coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mu1,p));
	coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k-1)));
      }
      if (k>2) k--;
    }
    else
    {
      for (l=k-2; l; l--)
      {
	if (!gcmp0(r=grndtoi(gcoeff(mu,k,l),&e)))
	{
	  u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[l]));
	  for (j=1; j<l; j++)
	    coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,l,j)));
	  coeff(mu,k,l)=lsub(gcoeff(mu,k,l),r);
	 }
      }
      k++;
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"lllgram1");
      tetpil=avma;
      sv=cgetg(4,t_VEC);
      sv[1]=lcopy(B); sv[2]=lcopy(u); sv[3]=lcopy(mu);
      sv=gerepile(av,tetpil,sv);
      B=(GEN)sv[1]; u=(GEN)sv[2]; mu=(GEN)sv[3];
    }
  }
  while (k<=n);
  tetpil=avma; return gerepile(av,tetpil,gcopy(u));
}

GEN
lllgramint(GEN x)
{
  return lllgramall0(x,lll_IM);
}

GEN
lllgramkerim(GEN x)
{
  return lllgramall0(x,lll_ALL);
}

/*  This routine is functionally similar to lllint when all = 0.
 *
 *    The input is an integer matrix mat (not necessarily square) whose
 *  columns are linearly independent.  It outputs another matrix T such that
 *  mat * T is partially reduced.  If all = 0, the size of mat * T is the
 *  same as the size of mat.  If all = 1 the number of columns of mat * T
 *  is at most equal to its number of rows.  A matrix M is said to be
 *  -partially reduced- if | m1 +- m2 | >= |m1| for any two distinct
 *  columns m1, m2, in M.
 *
 *    This routine is designed to quickly reduce lattices in which one row
 *  is huge compared to the other rows.  For example, when searching for a
 *  polynomial of degree 3 with root a mod p, the four input vectors might
 *  be the coefficients of
 *      X^3 - (a^3 mod p), X^2 - (a^2 mod p), X - (a mod p), p.
 *  All four constant coefficients are O(p) and the rest are O(1). By the
 *  pigeon-hole principle, the coefficients of the smallest vector in the
 *  lattice are O(p^(1/4)), Hence significant reduction of vector lengths
 *  can be anticipated.
 *
 *  		Peter Montgomery (July, 1994)
 *
 *  If flag = 1 complete the reduction using lllint, otherwise return
 *  partially reduced basis.
 */
GEN
lllintpartialall(GEN m, long flag)
{
  const long ncol = lg(m)-1;
  const long ltop1 = avma;
  long ncolnz;
  GEN tm1, tm2, mid, *gptr[4];

  if (typ(m) != t_MAT) err(typeer,"lllintpartialall");
  if (ncol <= 1) return idmat(ncol);

  {
    GEN dot11 = sqscali((GEN)m[1]);
    GEN dot22 = sqscali((GEN)m[2]);
    GEN dot12 = gscali((GEN)m[1], (GEN)m[2]);
    GEN tm  = idmat(2); /* For first two columns only */

    int progress = 0;
    long npass2 = 0;

/* Try to row reduce the first two columns of m.
 * Our best result so far is (first two columns of m)*tm.
 *
 * Initially tm = 2 x 2 identity matrix.
 * The inner products of the reduced matrix are in
 * dot11, dot12, dot22.
 */
    while (npass2 < 2 || progress)
    {
      GEN dot12new,q = gdivround(dot12, dot22);

      npass2++; progress = signe(q);
      if (progress)
      {
       /* Conceptually replace (v1, v2) by (v1 - q*v2, v2),
        * where v1 and v2 represent the reduced basis for the
        * first two columns of the matrix.
        *
        * We do this by updating tm and the inner products.
        *
        * An improved algorithm would look only at the leading
        * digits of dot11, dot12, dot22.  It would use
        * single-precision calculations as much as possible.
        */
        q = negi(q);
        dot12new = addii(dot12, mulii(q, dot22));
        dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
        dot12 = dot12new;
        tm[1] = (long)lincomb_integral(gun,q, (GEN)tm[1],(GEN)tm[2]);
      }

      /* Interchange the output vectors v1 and v2.  */
      gswap(dot11,dot22); swap(tm[1],tm[2]);

      /* Occasionally (including final pass) do garbage collection.  */
      if (npass2 % 8 == 0 || !progress)
      {
        gptr[0] = &dot11; gptr[1] = &dot12;
        gptr[2] = &dot22; gptr[3] = &tm;
        gerepilemany(ltop1, gptr, 4);
      }
    } /* while npass2 < 2 || progress */

    {
      long icol;
      GEN det12 = subii(mulii(dot11, dot22), mulii(dot12, dot12));

      tm1 = idmat(ncol);
      mid = cgetg(ncol+1, t_MAT);
      for (icol = 1; icol <= 2; icol++)
      {
        coeff(tm1,1,icol) = coeff(tm,1,icol);
	coeff(tm1,2,icol) = coeff(tm,2,icol);
        mid[icol] = (long)lincomb_integral(
           gcoeff(tm,1,icol),gcoeff(tm,2,icol), (GEN)m[1],(GEN)m[2]);
      }
      for (icol = 3; icol <= ncol; icol++)
      {
        GEN curcol = (GEN)m[icol];
	GEN dot1i = gscali((GEN)mid[1], curcol);
        GEN dot2i = gscali((GEN)mid[2], curcol);
       /* Try to solve
        *
        * ( dot11  dot12 ) (q1)    ( dot1i )
        * ( dot12  dot22 ) (q2)  = ( dot2i )
        *
        * Round -q1 and -q2 to the nearest integer.
        * Then compute curcol - q1*mid[1] - q2*mid[2].
        * This will be approximately orthogonal to the
        * first two vectors in the new basis.
        */
	GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
        GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));

        q1neg = gdivround(q1neg, det12);
        q2neg = gdivround(q2neg, det12);
        coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)),
				   gmul(q2neg, gcoeff(tm,1,2)));
        coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)),
				   gmul(q2neg, gcoeff(tm,2,2)));
        mid[icol] = ladd(curcol,
          lincomb_integral(q1neg,q2neg, (GEN)mid[1],(GEN)mid[2]));
      } /* for icol */
    } /* local block */
  }
  if (DEBUGLEVEL>4)
  {
    fprintferr("tm1 = "); outbeauterr(tm1);
    fprintferr("mid = "); outbeauterr(mid);
  }
  gptr[0] = &tm1; gptr[1] = &mid;
  gerepilemany(ltop1, gptr, 2);
  {
   /* For each pair of column vectors v and w in mid * tm2,
    * try to replace (v, w) by (v, v - q*w) for some q.
    * We compute all inner products and check them repeatedly.
    */
    const long ltop3 = avma; /* Excludes region with tm1 and mid */
    long icol, lim, reductions, npass = 0;
    GEN dotprd = cgetg(ncol+1, t_MAT);

    tm2 = idmat(ncol);
    for (icol=1; icol <= ncol; icol++)
    {
      long jcol;

      dotprd[icol] = lgetg(ncol+1,t_COL);
      for (jcol=1; jcol <= icol; jcol++)
	coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) =
          (long)gscali((GEN)mid[icol], (GEN)mid[jcol]);
    } /* for icol */
    lim = stack_lim(ltop3,1);
    for(;;)
    {
      reductions = 0;
      for (icol=1; icol <= ncol; icol++)
      {
	long ijdif, jcol, k1, k2;
	GEN codi, q;

        for (ijdif=1; ijdif < ncol; ijdif++)
	{
          const long previous_avma = avma;

          jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */
          k1 = (cmpii(gcoeff(dotprd,icol,icol),
		      gcoeff(dotprd,jcol,jcol) ) > 0)
		? icol : jcol; 	/* index of larger column */
	  k2 = icol + jcol - k1; 	/* index of smaller column */
	  codi = gcoeff(dotprd,k2,k2);
	  q = gcmp0(codi)? gzero: gdivround(gcoeff(dotprd,k1,k2), codi);

	  /* Try to subtract a multiple of column k2 from column k1.  */
	  if (gcmp0(q)) avma = previous_avma;
          else
	  {
	    long dcol;

	    reductions++; q = negi(q);
	    tm2[k1]=(long)
              lincomb_integral(gun,q, (GEN)tm2[k1], (GEN)tm2[k2]);
	    dotprd[k1]=(long)
              lincomb_integral(gun,q, (GEN)dotprd[k1], (GEN)dotprd[k2]);
	    coeff(dotprd, k1, k1) = laddii(gcoeff(dotprd,k1,k1),
				           mulii(q, gcoeff(dotprd,k2,k1)));
	    for (dcol = 1; dcol <= ncol; dcol++)
	      coeff(dotprd,k1,dcol) = coeff(dotprd,dcol,k1);
	  } /* if q != 0 */
        } /* for ijdif */
        if (low_stack(lim, stack_lim(ltop3,1)))
	{
          if(DEBUGMEM>1) err(warnmem,"lllintpartialall");
	  gptr[0] = &dotprd; gptr[1] = &tm2;
	  gerepilemany(ltop3, gptr, 2);
        }
      } /* for icol */
      if (!reductions) break;
      if (DEBUGLEVEL>4)
      {
	GEN diag_prod = dbltor(1.0);
	for (icol = 1; icol <= ncol; icol++)
	  diag_prod = gmul(diag_prod, gcoeff(dotprd,icol,icol));
        npass++;
	fprintferr("npass = %ld, red. last time = %ld, diag_prod = %Z\n\n",
	            npass, reductions, diag_prod);
      }
    } /* for(;;)*/

   /* Sort columns so smallest comes first in m * tm1 * tm2.
    * Use insertion sort.
    */
    for (icol = 1; icol < ncol; icol++)
    {
      long jcol, s = icol;

      for (jcol = icol+1; jcol <= ncol; jcol++)
	if (cmpii(gcoeff(dotprd,s,s),gcoeff(dotprd,jcol,jcol)) > 0) s = jcol;
      if (icol != s)
      { /* Exchange with proper column */
        /* Only diagonal of dotprd is updated */
        swap(tm2[icol], tm2[s]);
        swap(coeff(dotprd,icol,icol), coeff(dotprd,s,s));
      }
    } /* for icol */
    icol=1;
    while (icol <= ncol && !signe(gcoeff(dotprd,icol,icol))) icol++;
    ncolnz = ncol - icol + 1;
  } /* local block */

  if (flag)
  {
    if (ncolnz == lg((GEN)m[1])-1)
    {
      tm2 += (ncol-ncolnz);
      tm2[0]=evaltyp(t_MAT)|evallg(ncolnz+1);
    }
    else
    {
      tm1 = gmul(tm1, tm2);
      tm2 = lllint(gmul(m, tm1));
    }
  }
  if (DEBUGLEVEL>4)
    { fprintferr("lllintpartial output = "); outbeauterr(gmul(tm1, tm2)); }
  return gerepileupto(ltop1, gmul(tm1, tm2));
}

GEN
lllintpartial(GEN mat)
{
  return lllintpartialall(mat,1);
}

/********************************************************************/
/**                                                                **/
/**                   LINEAR & ALGEBRAIC DEPENDANCE                **/
/**                                                                **/
/********************************************************************/

GEN
lindep0(GEN x,long bit,long prec)
{
  if (!bit) return lindep(x,prec);
  if (bit>0) return lindep2(x,bit);
  return deplin(x);
}

GEN
lindep2(GEN x, long bit)
{
  long tx=typ(x),lx=lg(x),ly,i,j,flag,e,tetpil, av = avma;
  GEN re,im,p1,p2;

  if (! is_vec_t(tx)) err(typeer,"lindep2");
  if (lx<=2) return cgetg(1,t_VEC);
  re=greal(x); im=gimag(x); flag = !gcmp0(im);
  if (lx == 3)
  { /* independant over R ? */
    if (gexpo(gsub(gmul((GEN)re[1],(GEN)im[2]),
                   gmul((GEN)re[2],(GEN)im[1]))) > - bit)
    {
      avma = av; return cgetg(1, t_VEC);
    }
  }
  ly = flag? lx+2: lx+1;
  p2=cgetg(lx,t_MAT); bit = (long) (bit/L2SL10);
  for (i=1; i<lx; i++)
  {
    p1=cgetg(ly,t_COL); p2[i]=(long)p1;
    for (j=1; j<lx; j++) p1[j]=(i==j) ? un : zero;
    p1[lx]=lcvtoi(gshift((GEN)re[i],bit),&e);
    if (flag) p1[lx+1]=lcvtoi(gshift((GEN)im[i],bit),&e);
  }
  p1=gmul(p2,lllint(p2)); p1=(GEN)p1[1];
  p1[0]=evaltyp(t_VEC) | evallg(lx); tetpil=avma;
  return gerepile(av,tetpil,gcopy(p1));
}

#define quazero(x) (gcmp0(x) || (typ(x)==t_REAL && expo(x) < EXP))
GEN
lindep(GEN x, long prec)
{
  GEN *b,*be,*bn,**m,qzer;
  GEN c1,c2,c3,px,py,pxy,re,im,p3,p4,r,f,em;
  long i,j,fl,i1, lx = lg(x), tx = typ(x), n = lx-1;
  long av = avma, lim = stack_lim(av,1), av0,av1,tetpil;
  const long EXP = - bit_accuracy(prec) + 2*n;

  if (! is_vec_t(tx)) err(typeer,"lindep");
  if (lx<=2) return cgetg(1,t_VEC);
  x = gmul(x, realun(prec));
  re=greal(x); im=gimag(x);
  if (lx == 3)
  { /* independant over R ? */
    if (gexpo(gsub(gmul((GEN)re[1],(GEN)im[2]),
                   gmul((GEN)re[2],(GEN)im[1]))) > - bit_accuracy(prec))
    {
      avma = av; return cgetg(1, t_VEC);
    }
  }
  qzer = new_chunk(lx);
  b = (GEN*)idmat(n);
  be= (GEN*)new_chunk(lx);
  bn= (GEN*)new_chunk(lx);
  m = (GEN**)new_chunk(lx);
  for (i=1; i<=n; i++)
  {
    bn[i]=cgetr(prec+1);
    be[i]=cgetg(lx,t_COL);
    m[i] = (GEN*)new_chunk(lx);
    for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1);
    for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1);
  }
  px=sqscal(re);
  py=sqscal(im); pxy=gscal(re,im);
  p3=mpsub(mpmul(px,py),gsqr(pxy));
  if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p3);
  av0 = av1 = avma;
  for (i=1; i<=n; i++)
  {
    p4 = gscal(b[i],re);
    if (fl) p4=gmul(gdiv(p4,px),re);
    else
    {
      GEN p5,p6,p7;
      p5=gscal(b[i],im);
      p6=gdiv(gsub(gmul(py,p4),gmul(pxy,p5)),p3);
      p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p4)),p3);
      p4=gadd(gmul(p6,re),gmul(p7,im));
    }
    if (tx!=t_COL) p4=gtrans(p4);
    p4=gsub(b[i],p4);
    for (j=1; j<i; j++)
      if (qzer[j]) affrr(bn[j],m[i][j]);
      else
      {
        gdivz(gscal(b[i],be[j]),bn[j],m[i][j]);
        p4=gsub(p4,gmul(m[i][j],be[j]));
      }
    gaffect(p4,be[i]); affrr(sqscal(be[i]),bn[i]);
    qzer[i]=quazero(bn[i]); avma=av1;
  }
  while (qzer[n])
  {
    long e;
    if (DEBUGLEVEL>9)
    {
      fprintferr("qzer[%ld]=%ld\n",n,qzer[n]);
      for (i1=1; i1<=n; i1++)
	for (i=1; i<i1; i++) output(m[i1][i]);
    }
    em=bn[1]; j=1;
    for (i=2; i<n; i++)
    {
      p3=shiftr(bn[i],i);
      if (cmprr(p3,em)>0) { em=p3; j=i; }
    }
    i=j; i1=i+1;
    avma = av1; r = grndtoi(m[i1][i], &e);
    if (e >= 0) err(talker,"precision too low in lindep");
    r  = negi(r);
    p3 = lincomb_integral(gun,r, b[i1],b[i]);
    av1 = avma;
    b[i1]=b[i]; b[i]=p3; f=addir(r,m[i1][i]);
    for (j=1; j<i; j++)
      if (!qzer[j])
      {
        p3=mpadd(m[i1][j],mulir(r,m[i][j]));
        affrr(m[i][j],m[i1][j]); mpaff(p3,m[i][j]);
      }
    c1=addrr(bn[i1],mulrr(gsqr(f),bn[i]));
    if (!quazero(c1))
    {
      c2=divrr(mulrr(bn[i],f),c1); affrr(c2,m[i1][i]);
      c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]);
      affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0;
      for (j=i+2; j<=n; j++)
      {
        p3=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2));
        subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]);
        affrr(p3,m[j][i]);
      }
    }
    else
    {
      qzer[i1]=qzer[i]; qzer[i]=1;
      affrr(bn[i],bn[i1]); affrr(c1,bn[i]);
      for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]);
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"lindep");
      b = (GEN*)gerepileupto(av0, gcopy((GEN)b));
      av1 = avma;
    }
  }
  p3=cgetg(lx,t_COL); p3[n]=un; for (i=1; i<n; i++) p3[i]=zero;
  p4 = (GEN)b; p4[0] = evaltyp(t_MAT) | evallg(lx);
  p4=gauss(gtrans(p4),p3); tetpil=avma;
  return gerepile(av,tetpil,gtrans(p4));
}

GEN
algdep0(GEN x, long n, long bit, long prec)
{
  long tx=typ(x),av,i,k;
  GEN y,p1;

  if (! is_scalar_t(tx)) err(typeer,"algdep0");
  if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
  if (gcmp0(x)) return gzero;
  if (!n) return gun;

  av=avma; p1=cgetg(n+2,t_COL); p1[1]=un;
  for (i=2; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);

  p1 = bit? lindep2(p1,bit): lindep(p1,prec);
  if (lg(p1) < 2)
    err(talker,"higher degree than expected in algdep");

  y=cgetg(n+3,t_POL);
  y[1] = evalsigne(1) | evalvarn(0);
  k=1; while (gcmp0((GEN)p1[k])) k++;
  for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
  normalizepol_i(y, n+4-k);
  y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
  return gerepileupto(av,y);
}

GEN
algdep2(GEN x, long n, long bit)
{
  return algdep0(x,n,bit,0);
}

GEN
algdep(GEN x, long n, long prec)
{
  return algdep0(x,n,0,prec);
}

/********************************************************************/
/**                                                                **/
/**                   INTEGRAL KERNEL (LLL REDUCED)                **/
/**                                                                **/
/********************************************************************/

GEN
matkerint0(GEN x, long flag)
{
  switch(flag)
  {
    case 0: return kerint(x);
    case 1: return kerint1(x);
    case 2: return kerint2(x);
    default: err(flagerr,"matkerint");
  }
  return NULL; /* not reached */
}

GEN
kerint1(GEN x)
{
  long av=avma,tetpil;
  GEN p1,p2;

  p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma;
  return gerepile(av,tetpil,gmul(p1,p2));
}

GEN
kerint2(GEN x)
{
  long lx=lg(x), i,j,av,av1;
  GEN g,p1;

  if (typ(x) != t_MAT) err(typeer,"kerint2");
  av=avma; g=cgetg(lx,t_MAT);
  for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
  for (i=1; i<lx; i++)
    for (j=1; j<=i; j++)
      coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)x[i],(GEN)x[j]);
  g=lllgramall0(g,lll_KER); p1=lllint(g);
  av1=avma; return gerepile(av,av1,gmul(g,p1));
}

static GEN
lllall0(GEN x, long flag)
{
  long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
  GEN u,B,L,q,r,h,la,p1,p2,p4,fl;

  if (typ(x) != t_MAT) err(typeer,"lllall0");
  n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM);
  fl = new_chunk(lx);

  av=avma; lim=stack_lim(av,1); x=dummycopy(x);
  B=gscalcol(gun, lx);
  L=cgetg(lx,t_MAT);
  for (k=lg(x[1]),j=1; j<lx; j++)
  {
    for (i=1; i<k; i++)
      if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
    fl[j] = 0; L[j] = (long)zerocol(n);
  }
  k=2; h=idmat(n); kmax=1;
  u=sqscali((GEN)x[1]);
  if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
  for(;;)
  {
    if (k > kmax)
    {
      kmax = k;
      for (j=1; j<=k; j++)
      {
	if (j==k || fl[j])
	{
          long av1 = avma;
	  u=gscali((GEN)x[k],(GEN)x[j]);
	  for (i=1; i<j; i++)
	    if (fl[i])
              u = divii(subii(mulii((GEN)B[i+1],u),
                              mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
                        (GEN)B[i]);
          u = gerepileuptoint(av1, u);
	  if (j<k) coeff(L,k,j)=(long)u;
	  else
	  {
	    if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
	    else { B[k+1]=B[k]; fl[k]=0; }
	  }
	}
      }
    }
    if (fl[k-1] && !fl[k])
    {
      u = shifti(gcoeff(L,k,k-1),1);
      if (absi_cmp(u, (GEN)B[k]) > 0)
      {
	q = truedvmdii(addii(u,(GEN)B[k]),shifti((GEN)B[k],1), NULL);
        r = negi(q);
        h[k] = (long)lincomb_integral(gun,r, (GEN)h[k],(GEN)h[k-1]);
        x[k] = (long)lincomb_integral(gun,r, (GEN)x[k],(GEN)x[k-1]);
	coeff(L,k,k-1)=laddii(gcoeff(L,k,k-1),mulii(r,(GEN)B[k]));
	for (i=1; i<k-1; i++)
	  coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,k-1,i)));
      }
      la=gcoeff(L,k,k-1); p4=sqri(la);
      swap(h[k-1], h[k]);
      swap(x[k-1], x[k]);
      for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
      if (signe(la))
      {
	p2=(GEN)B[k]; p1=divii(p4,p2);
	B[k+1]=B[k]=(long)p1;
	for (i=k+1; i<=kmax; i++)
	  coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
	for (j=k+1; j<kmax; j++)
	  for (i=j+1; i<=kmax; i++)
	    coeff(L,i,j)=ldivii(mulii((GEN)p1,gcoeff(L,i,j)),p2);
	for (i=k+2; i<=kmax+1; i++)
	  B[i]=ldivii(mulii((GEN)p1,(GEN)B[i]),p2);
      }
      else
      {
	for (i=k+1; i<=kmax; i++)
	{ coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
	B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
      }
      if (k>2) k--;
    }
    else
    {
      for (l=k-1; l>=1; l--)
      {
        u = shifti(gcoeff(L,k,l),1);
	if (absi_cmp(u,(GEN)B[l+1]) > 0)
	{
	  q = truedvmdii(addii(u,(GEN)B[l+1]),shifti((GEN)B[l+1],1), NULL);
          r = negi(q);
	  h[k] = (long)lincomb_integral(gun,r,(GEN)h[k],(GEN)h[l]);
	  x[k] = (long)lincomb_integral(gun,r,(GEN)x[k],(GEN)x[l]);
	  coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,(GEN)B[l+1]));
	  for (i=1; i<l; i++)
	    coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
	 }
      }
      if (++k > n) break;
    }
    if (low_stack(lim, stack_lim(av,1)))
    {
      GEN *gptr[4];
      if(DEBUGMEM>1) err(warnmem,"lllall0");
      gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
      gptr[3]=&x; gerepilemany(av,gptr,4);
    }
  }
  tetpil=avma;
  return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
}

GEN
kerint(GEN x)
{
  long av=avma,av1;
  GEN g,p1;

  g=lllall0(x,lll_KER); if (lg(g)==1) return g;
  p1=lllint(g); av1=avma;
  return gerepile(av,av1,gmul(g,p1));
}

/********************************************************************/
/**                                                                **/
/**                        POLRED & CO.                            **/
/**                                                                **/
/********************************************************************/
/* remove duplicate polynomials in y, updating a (same indexes), in place */
static long
remove_duplicates(GEN y, GEN a)
{
  long k,i, nv = lg(y), av = avma;
  GEN z;

  if (nv < 2) return nv;
  z = new_chunk(3);
  z[1] = (long)y;
  z[2] = (long)a; (void)sort_factor(z, gcmp);
  for  (k=1, i=2; i<nv; i++)
    if (!gegal((GEN)y[i], (GEN)y[k]))
    {
      k++;
      a[k] = a[i];
      y[k] = y[i];
    }
  nv = k+1; setlg(a,nv); setlg(y,nv);
  avma = av; return nv;
}

/* in place; make sure second non-zero coeff is negative (choose x or -x) */
int
canon_pol(GEN z)
{
  long i,s;

  for (i=lgef(z)-2; i>=2; i-=2)
  {
    s = signe(z[i]);
    if (s > 0)
    {
      for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
      return -1;
    }
    if (s) return 1;
  }
  return 0;
}

static GEN
pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta)
{
  long i,j, v = varn(x), n = lg(base);
  GEN p1,p2,p3,y, a = cgetg(n,t_VEC);

  for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]);
  y=cgetg(n,t_VEC);
  for (i=1; i<n; i++)
  {
    if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
    p1=(GEN)a[i]; p3=content(p1);
    if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3);
    p1 = caract2(x,p1,v);
    if (p3)
      for (p2=gun, j=lgef(p1)-2; j>=2; j--)
      {
        p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2);
      }
    p2=modulargcd(derivpol(p1),p1); p3=leading_term(p2);
    if (!gcmp1(p3)) p2=gdiv(p2,p3);
    p1 = gdiv(p1,p2);
    if (canon_pol(p1) < 0 && pta)
      a[i] = (long) gneg_i((GEN)a[i]);
    y[i] = (long)p1; if (DEBUGLEVEL>=4) outerr(p1);
  }
  (void)remove_duplicates(y,a);
  if (pta) *pta = a;
  return y;
}

GEN
nf_get_T2(GEN base, GEN polr)
{
  long i,j, n = lg(base);
  GEN p1,p2=cgetg(n,t_MAT);

  for (i=1; i<n; i++)
  {
    p1=cgetg(n,t_COL); p2[i]=(long)p1;
    for (j=1; j<n; j++)
      p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]);
  }
  return mulmat_real(gconj(gtrans(p2)),p2);
}

/* Return the base change matrix giving the an LLL-reduced basis for the
 * maximal order of the nf given by x. Expressed in terms of the standard
 * HNF basis (as polynomials) given in base (ignored if x is an nf)
 */
GEN
LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec)
{
  GEN T2,p1, x = *ptx;
  int totally_real,n,i,j;

  if (typ(x) != t_POL)
  {
    p1=checknf(x); *ptx=x=(GEN)p1[1];
    base=(GEN)p1[7]; n=lgef(x)-3;
    totally_real = !signe(gmael(p1,2,2));
    T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2);
  }
  else
  {
    n=lgef(x)-3; totally_real = (!prec || sturm(x)==n);
    if (!totally_real)
      T2 = nf_get_T2(base,polr? polr: roots(x,prec));
    else
    { /* totally real */
      GEN ptrace=cgetg(n+2,t_VEC);
      long k;

      ptrace[2]=lstoi(n);
      for (k=2; k<=n; k++)
      { /* cf polsym */
	GEN y = x + (n-k+1);
	p1 = gmulsg(k-1,(GEN)y[2]);
	for (i=3; i<=k; i++)
	  p1 = gadd(p1,gmul((GEN)y[i],(GEN)ptrace[i]));
	ptrace[i] = lneg(p1);
      }
      T2=cgetg(n+1,t_MAT);
      for (i=1; i<=n; i++)
      {
	p1=cgetg(n+1,t_COL); T2[i]=(long)p1;
	for (j=1; j<i ; j++) p1[j] = coeff(T2,i,j);
	for (   ; j<=n; j++)
	{ /* cf quicktrace */
	  GEN p2 = gres(gmul((GEN)base[i],(GEN)base[j]),x);
          GEN p3 = gzero;
	  for (k=lgef(p2)-1; k>1; k--)
	    p3 = gadd(p3, gmul((GEN)p2[k],(GEN)ptrace[k]));
	  p1[j]=(long)p3;
	}
      }
    }
  }
  if (totally_real) return lllgramint(T2);
  for (i=1; ; i++)
  {
    if ((p1 = lllgramintern(T2,100,1,prec))) return p1;
    if (i == MAXITERPOL) err(accurer,"allpolred");
    prec=(prec<<1)-2; 
    if (DEBUGLEVEL) err(warnprec,"allpolred",prec);
    T2=nf_get_T2(base,roots(x,prec));
  }
}

/* x can be a polynomial, but also an nf or a bnf */
static GEN
allpolred(GEN x, GEN *pta, long code, long prec)
{
  GEN y,p1, base = NULL, polr = NULL;
  long av = avma;

  if (typ(x) == t_POL)
  {
    if (!signe(x)) return gcopy(x);
    check_pol_int(x);
    if (!gcmp1(leading_term(x)))
      err(impl,"allpolred for nonmonic polynomials");
    base = allbase4(x,code,&p1,NULL); /* p1 is junk */
  }
  else
  {
    long i = lg(x);
    if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
    { /* polynomial + integer basis */
      base=(GEN)x[2]; x=(GEN)x[1];
    }
  }
  p1 = LLL_nfbasis(&x,polr,base,prec);
  y = pols_for_polred(x,base,p1,pta);
  if (pta)
  {
    GEN *gptr[2]; gptr[0]=&y; gptr[1]=pta;
    gerepilemany(av,gptr,2); return y;
  }
  return gerepileupto(av,y);
}

GEN
polred0(GEN x,long flag, GEN p, long prec)
{
  GEN y;
  long smll = (flag & 1);

  if (p && !gcmp0(p)) smll=(long) p; /* factored polred */
  if (flag & 2) /* polred2 */
  {
    y=cgetg(3,t_MAT);
    y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec);
    return y;
  }
  return allpolred(x,NULL,smll,prec);
}

GEN
ordred(GEN x, long prec)
{
  GEN p1,p2,p3,base;
  long n=lgef(x)-3,i,j,av=avma,v = varn(x),tetpil;

  if (typ(x) != t_POL) err(typeer,"ordred");
  if (!signe(x)) return gcopy(x);
  if (!gcmp1((GEN)x[n+2])) err(impl,"ordred for nonmonic polynomials");

  p2=roots(x,prec); p3=cgetg(n+1,t_MAT);

  base=cgetg(n+1,t_VEC); /* power basis */
  for (i=1; i<=n; i++)
  {
    base[i]=lpowgs(polx[v],i-1);
    p1=cgetg(n+1,t_COL); p3[i]=(long)p1;
    for (j=1; j<=n; j++)
      p1[j]=lpuigs((GEN)p2[j],i-1);
  }
  p2 = mulmat_real(gconj(gtrans(p3)),p3);
  p1 = pols_for_polred(x,base,lllgram(p2,prec), NULL);

  tetpil=avma; return gerepile(av,tetpil, gcopy(p1));
}

GEN roots_to_pol_r1r2(GEN a, long r1, long v);
static GEN chk_basis_embed;
long chk_r1, chk_v;

static GEN
init_chk(GEN nf, GEN mat, GEN bound)
{
  GEN M = gmael(nf,5,1);
  long n = lg(nf[7])-1, prec,prec2;

  if (!bound)
  {
    chk_v = varn(nf[1]);
    chk_r1 = itos(gmael(nf,2,1));
    chk_basis_embed = gmul(M, mat);
    return gmul((GEN)nf[7],mat);
  }

  /* should be [max_k C^n_k (B/k) ^ k] ^ (1/2) */
  bound = gsqrt(gpowgs(bound, n), 3);
  prec2 = (1 + (gexpo(bound) >> TWOPOTBITS_IN_LONG));
  if (prec2 < 0) prec2 = 0;
  prec = 3 + prec2;
  prec2= (long)nfnewprec(nf,-1);
  if (DEBUGLEVEL)
    fprintferr("init_chk: estimated prec = %ld (initially %ld)\n",
                 prec, prec2);
  if (prec > prec2) return NULL;
  if (prec < prec2) chk_basis_embed = gprec_w(chk_basis_embed, prec);
  return (GEN)prec;
}

static GEN
checkgenerator(GEN x)
{
  long l,i,av = avma;
  GEN g = gmul(chk_basis_embed,x);
  GEN r = gneg((GEN)g[1]);
  long epsbit = 5 - bit_accuracy(gprecision(r));
  
  l = lg(g)-1;
  for (i=2; i<=l; i++)
    if (gexpo(gadd((GEN)g[i], r)) < epsbit) { avma=av; return NULL; }
  
  g = roots_to_pol_r1r2(g, chk_r1, chk_v);
  g = ground(g);
  if (lgef(modulargcd(g,derivpol(g))) > 3) { avma=av; return NULL; }
  if (DEBUGLEVEL>3) fprintferr("  generator: %Z\n",g);
  return g;
}

/* no garbage collecting, done in polredabs0 */
static GEN
findmindisc(GEN nf, GEN y, GEN a, GEN phimax, long flun)
{
  long i,k, c = lg(y);
  GEN v,dmin,z,beta,discs = cgetg(c,t_VEC);

  for (i=1; i<c; i++) discs[i] = labsi(discsr((GEN)y[i]));
  v = sindexsort(discs);
  k = v[1]; dmin = (GEN)discs[k]; z = (GEN)y[k]; beta = (GEN)a[k];
  for (i=2; i<c; i++)
  {
    k = v[i];
    if (!egalii((GEN)discs[k],dmin)) break;
    if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; }
  }
  if (flun & nf_RAW)
  {
    y=cgetg(3,t_VEC);
    y[1]=lcopy(z);
    y[2]=lcopy(beta);
  }
  else if (phimax)
  {
    y=cgetg(3,t_VEC);
    y[1]=lcopy(z);
    beta=polymodrecip(gmodulcp(beta,(GEN)nf[1]));
    y[2]=(long)poleval(phimax,beta);
  }
  else y = gcopy(z);
  return y;
}

/* no garbage collecting, done in polredabs0 */
static GEN
storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun)
{
  GEN p1,y,beta;

  if (flun & nf_RAW)
  {
    long i, c = lg(z);
    y=cgetg(c,t_VEC);
    for (i=1; i<c; i++)
    {
      p1=cgetg(3,t_VEC); y[i]=(long)p1;
      p1[1]=lcopy((GEN)z[i]);
      p1[2]=lcopy((GEN)a[i]);
    }
  }
  else if (phimax)
  {
    long i, c = lg(z);
    beta = new_chunk(c);
    for (i=1; i<c; i++)
      beta[i] = (long)polymodrecip(gmodulcp((GEN)a[i],(GEN)nf[1]));

    y=cgetg(c,t_VEC);
    for (i=1; i<c; i++)
    {
      p1=cgetg(3,t_VEC); y[i]=(long)p1;
      p1[1]=lcopy((GEN)z[i]);
      p1[2]=(long)poleval(phimax,(GEN)beta[i]);
    }
  }
  else y = gcopy(z);
  return y;
}

GEN
polredabs0(GEN x, long flun, long prec)
{
  long i,nv, av = avma;
  GEN nf,v,y,a,phimax;
  GEN (*storepols)(ANYARG);

  if ((ulong)flun >= 16) err(flagerr,"polredabs");
  nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec);
  if (lg(nf) == 3)
  {
    phimax = lift_to_pol((GEN)nf[2]);
    nf = (GEN)nf[1];
  }
  else
    phimax = (flun & nf_ORIG)? polx[0]: (GEN)NULL;
  prec = (long)nfnewprec(nf,0);

  nv = lgef(nf[1])==4? 1: 2;
  for (i=1; ; i++)
  {
    v = fincke_pohst(nf,NULL,stoi(5000),3,prec, &checkgenerator);
    if (v) break;
    if (i==MAXITERPOL) err(accurer,"polredabs0");
    prec = (prec<<1)-2; nf = nfnewprec(nf,prec);
    if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
  }
  a = (GEN)v[2];
  y = (GEN)v[1];
  nv = lg(a);
  for (i=1; i<nv; i++)
    if (canon_pol((GEN)y[i]) < 0 && phimax)
      a[i] = (long) gneg_i((GEN)a[i]);
  nv = remove_duplicates(y,a);

  if (DEBUGLEVEL)
    { fprintferr("%ld minimal vectors found.\n",nv-1); flusherr(); }
  if (nv >= 10000) flun &= (~nf_ALL); /* should not happen */
  storepols = (flun & nf_ALL)? storeallpols: findmindisc;

  if (DEBUGLEVEL) fprintferr("\n");
  if (nv==1)
  {
    x=(GEN)nf[1];
    y = cgetg(2,t_VEC); y[1]=(long)x;
    a = cgetg(2,t_VEC); a[1]=(long)polx[varn(x)];
  }
  return gerepileupto(av, storepols(nf,y,a,phimax,flun));
}

GEN
polredabsall(GEN x, long flun, long prec)
{
  return polredabs0(x, flun | nf_ALL, prec);
}

GEN
polredabs(GEN x, long prec)
{
  return polredabs0(x,nf_REGULAR,prec);
}

GEN
polredabs2(GEN x, long prec)
{
  return polredabs0(x,nf_ORIG,prec);
}

GEN
polredabsnored(GEN x, long prec)
{
  return polredabs0(x,nf_NORED,prec);
}

GEN
polred(GEN x, long prec)
{
  return allpolred(x,(GEN*)0,0,prec);
}

GEN
smallpolred(GEN x, long prec)
{
  return allpolred(x,(GEN*)0,1,prec);
}

GEN
factoredpolred(GEN x, GEN p, long prec)
{
  return allpolred(x,(GEN*)0,(long)p,prec);
}

GEN
polred2(GEN x, long prec)
{
  GEN y=cgetg(3,t_MAT);

  y[2]= (long) allpolred(x,(GEN*)(y+1),0,prec);
  return y;
}

GEN
smallpolred2(GEN x, long prec)
{
  GEN y=cgetg(3,t_MAT);

  y[2]= (long) allpolred(x,(GEN*)(y+1),1,prec);
  return y;
}

GEN
factoredpolred2(GEN x, GEN p, long prec)
{
  GEN y=cgetg(3,t_MAT);

  y[2]= (long) allpolred(x,(GEN*)(y+1),(long)p,prec);
  return y;
}

GEN makebasis(GEN nf,GEN pol);
/* relative polredabs. Returns
 * flag = 0: relative polynomial
 * flag = 1: relative polynomial + element
 * flag = 2: absolute polynomial */
GEN
rnfpolredabs(GEN nf, GEN relpol, long flag, long prec)
{
  GEN p1,bpol,rnf,elt,pol;
  long va, av = avma;

  if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
  nf=checknf(nf); va = varn(relpol);
  if (DEBUGLEVEL>1) timer2();
  p1 = makebasis(nf, unifpol(nf,relpol,1));
  rnf = (GEN)p1[3];
  if (DEBUGLEVEL>1) 
  {
    msgtimer("absolute basis");
    fprintferr("original absolute generator: %Z\n",p1[1]);
  }
  p1 = polredabs0(p1, nf_RAW | nf_NORED, prec);
  bpol = (GEN)p1[1];
  if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",bpol);
  if (flag==2) return gerepileupto(av,bpol);

  elt = rnfelementabstorel(rnf,(GEN)p1[2]);
  p1=cgetg(3,t_VEC);
#if 1
  pol = rnfcharpoly(nf,relpol,elt,va);
#else
{
  long i;
  p1=(GEN)nffactor(nf,bpol)[1];
  for (i=lg(p1)-1; i; i--)
    if(gcmp0(gsubst((GEN)p1[i],va,elt))) { pol=(GEN)p1[i]; break; }
  if (!i) err(bugparier,"rnfpolredabs (pol not found)");
}
#endif

  if (!flag) p1 = pol;
  else
  {
    p1[1]=(long)pol;
    p1[2]=(long)polymodrecip(elt);
  }
  return gerepileupto(av,p1);
}

/********************************************************************/
/**                                                                **/
/**                              MINIM                             **/
/**                                                                **/
/********************************************************************/
long addcolumntomatrix(long *V,long n,long r,GEN *INVP,long *L);
GEN gmul_mat_smallvec(GEN x, GEN y, long hx, long ly);

/* Minimal vectors for the integral definite quadratic form: a.
 * Result u:
 *   u[1]= Number of vectors of square norm <= BORNE
 *   u[2]= maximum norm found
 *   u[3]= list of vectors found (at most STOCKMAX)
 *
 *  If BORNE = gzero: Minimal non-zero vectors.
 *  flag = min_ALL,   as above
 *  flag = min_FIRST, exits when first suitable vector is found.
 *  flag = min_PERF,  only compute rank of the family of v.v~ (v min.)
 */
static GEN
minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
{
  GEN res,p1,u,r,liste,gnorme,gnorme_max,invp,V, *gptr[2];
  long n = lg(a), av0 = avma, av1,av,tetpil,lim, i,j,k,s,maxrank,*x;
  double p,borne,*v,*y,*z,**q, eps = 0.000001;

  switch(flag)
  {
    case min_FIRST: res = cgetg(3,t_VEC); break;
    case min_ALL: res = cgetg(4,t_VEC);
  }
  av=avma;

  x = (long*)    new_chunk(n);
  q = (double**) new_chunk(n);

  /* correct alignment for the following */
  s = avma % sizeof(double); avma -= s;
  if (avma<bot) err(errpile);

  s = (n * sizeof(double))/sizeof(long);
  y = (double*) new_chunk(s);
  z = (double*) new_chunk(s);
  v = (double*) new_chunk(s);
  for (j=1; j<n; j++) q[j] = (double*) new_chunk(s);
  av1=avma;

  u = lllgramint(a); a = qf_base_change(a,u,1);
  if (lg(a) != n)
    err(talker,"not a definite form in minim00");

  n--;
  a = gmul(a, realun(DEFAULTPREC)); r = sqred1(a);
  if (DEBUGLEVEL>4) { fprintferr("minim: r = "); outerr(r); }
  for (j=1; j<=n; j++)
  {
    v[j] = rtodbl(gcoeff(r,j,j));
    for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
  }

  if (flag==min_PERF || gcmp0(BORNE))
  {
    double c, b = rtodbl(gcoeff(a,1,1));

    for (i=2; i<=n; i++)
      { c=rtodbl(gcoeff(a,i,i)); if (c<b) b=c; }
    borne = b+eps;
    BORNE = ground(dbltor(borne));
    gnorme_max = NULL;
  }
  else
  {
    BORNE = gfloor(BORNE);
    borne = gtodouble(BORNE)+eps;
    gnorme_max = gzero;
  }

  switch(flag)
  {
    case min_ALL:
      maxrank=itos(STOCKMAX);
      liste = new_chunk(1+maxrank);
      break;
    case min_PERF:
      BORNE = gerepileupto(av1,BORNE);
      maxrank = (n*(n+1))>>1;
      liste = new_chunk(1+maxrank); V = new_chunk(1+maxrank);
      for (i=1; i<=maxrank; i++) liste[i]=0;
  }

  s=0; av1=avma; lim = stack_lim(av1,1);
  k = n; y[n] = z[n] = 0;
  x[n] = (long) sqrt(borne/v[n]+eps);
  if (flag == min_PERF) invp = idmat(maxrank);
  for(;;)
  {
    do
    {
      if (k>1)
      {
	k--; z[k]=0;
	for (j=k+1; j<=n; j++) z[k] += q[k][j]*x[j];
	p = x[k+1]+z[k+1];
	y[k] = y[k+1] + p*p*v[k+1];
	x[k] = (long) floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
      }
      for(;;)
      {
	p=x[k]+z[k];
	if (y[k] + p*p*v[k] <= borne+eps) break;
	k++; x[k]--;
      }
    }
    while (k>1);
    if (! x[1] && y[1]<=eps) break;
    p = x[1]+z[1]; gnorme = ground( dbltor(y[1]+ p*p*v[1]) );
    if (gnorme_max)
      { if (gcmp(gnorme,gnorme_max) > 0) gnorme_max=gnorme; }
    else
    {
      if (gcmp(gnorme,BORNE) < 0)
      {
	borne=gtodouble(gnorme); s=0;
        affii(gnorme,BORNE); avma=av1;
	if (flag == min_PERF) invp = idmat(maxrank);
      }
    }

    switch(flag)
    {
      case min_ALL:
        s++;
        if (s<=maxrank)
        {
          p1 = new_chunk(n+1); liste[s] = (long) p1;
          for (i=1; i<=n; i++) p1[i] = x[i];
        }
        break;

      case min_FIRST:
        if (! gnorme_max || gcmp(gnorme,BORNE)>0) break;

        tetpil=avma; gnorme = icopy(gnorme); r = gmul_mat_smallvec(r,x,n,n);
        gptr[0]=&gnorme; gptr[1]=&r; gerepilemanysp(av,tetpil,gptr,2);
        res[1]=(long)gnorme;
        res[2]=(long)r; return res;

      case min_PERF:
      {
        long av2=avma, I=1, newran;

        for (i=1; i<=n; i++)
          for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
        newran = addcolumntomatrix(V,maxrank,s,&invp,liste);
        if (newran == s)
        {
          avma=av2;
          if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
        }
        else
        {
          if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
          s = newran;
          if (s == maxrank)
          {
            if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
            avma=av0; return stoi(s);
          }

          if (low_stack(lim, stack_lim(av1,1)))
          {
            if(DEBUGMEM>1) err(warnmem,"minim00");
            if (DEBUGLEVEL>1)
            { 	
              fprintferr("\ngerepile in qfperfection. rank>=%ld\n",s);
              flusherr();
            }
            tetpil=avma; invp = gerepile(av1,tetpil,gcopy(invp));
          }
        }
      }
    }
    x[1]--;
  }
  switch(flag)
  {
    case min_FIRST:
      avma=av0; return cgetg(1,t_VEC);
    case min_PERF:
      if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
      avma=av0; return stoi(s);
  }
  k = min(s,maxrank);

  tetpil = avma; p1=cgetg(k+1,t_MAT);
  for (j=1; j<=k; j++)
    p1[j] = (long) gmul_mat_smallvec(u,(GEN)liste[j],n,n);
  liste = p1;
  r = gnorme_max? gnorme_max: BORNE;

  r=icopy(r); gptr[0]=&r; gptr[1]=&liste;
  gerepilemanysp(av,tetpil,gptr,2);
  res[1]=lstoi(s<<1);
  res[2]=(long)r;
  res[3]=(long)liste; return res;
}

GEN
minim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
{
  switch(flag)
  {
    case 0: return minim00(a,borne,stockmax,min_ALL);
    case 1: return minim00(a,borne,gzero   ,min_FIRST);
    case 2: return fincke_pohst(a,borne,stockmax,0,prec,NULL);
    default: err(flagerr,"qfminim");
  }
  return NULL; /* not reached */
}

GEN
minim(GEN a, GEN borne, GEN stockmax)
{
  return minim00(a,borne,stockmax,min_ALL);
}

GEN
minim2(GEN a, GEN borne, GEN stockmax)
{
  return minim00(a,borne,stockmax,min_FIRST);
}

GEN
perf(GEN a)
{
  return minim00(a,gzero,gzero,min_PERF);
}

/* programme general pour les formes quadratiques definies positives
 * quelconques (a coeffs reels). On doit avoir BORNE != 0; on n'entre dans
 * cette fonction qu'a partir de fincke_pohst (la reduction LLL n'est donc
 * pas faite ici). Si flag >= 2, on s'arrete des que stockmax est atteint.
 * Si flag&1 == 1, pas d'erreur dans sqred1intern. */
static GEN
smallvectors(GEN a, GEN BORNE, GEN STOCKMAX, long flag, long prec,
             GEN (*check)(GEN))
{
  long av = avma,av1,av2,lim,N,n,i,j,k,s,stockmax,epsbit,checkcnt = 0;
  GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms;

  N=lg(a); n=N-1; stockmax=itos(STOCKMAX);
  if (DEBUGLEVEL)
    fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3));

  lim=stack_lim(av,1); q=sqred1intern(a,flag&1);
  if (q == NULL) { avma=av; return NULL; }
  if (DEBUGLEVEL>5) fprintferr("q = %Z",q);
  epsbit = bit_accuracy(prec) >> 1;
  eps=realun(prec); setexpo(eps,-epsbit);
  alpha = dbltor(0.95);
  normax1=gzero;
  borne1=gadd(BORNE,eps);
  borne2=mpmul(borne1,alpha);
  v=cgetg(N,t_VEC);

  av2=avma;
  S = cgetg(stockmax+1,t_MAT);
  norms = cgetg(stockmax+1,t_VEC);
  x=cgetg(N,t_COL);
  y=cgetg(N,t_COL);
  z=cgetg(N,t_COL);
  for (i=1; i<N; i++) { v[i] = coeff(q,i,i); x[i]=y[i]=z[i] = zero; }

  x[n] = lmpent(mpsqrt(gdiv(borne1,(GEN)v[n])));
  if (DEBUGLEVEL>3) { fprintferr("\nx[%ld] = %Z\n",n,x[n]); flusherr(); }

  s=0; k=n;
  for(;;)
  {
    do
    {
      int fl=0;
      if (k>1)
      {
        av1=avma; k--;
        p1 = mpmul(gcoeff(q,k,k+1),(GEN)x[k+1]);
	for (j=k+2; j<N; j++)
	  p1 = mpadd(p1, mpmul(gcoeff(q,k,j),(GEN)x[j]));
        z[k] = lpileupto(av1,p1);

        av1=avma;
	p1 = gsqr(mpadd((GEN)x[k+1],(GEN)z[k+1]));
	y[k] = lpileupto(av1, mpadd((GEN)y[k+1], mpmul(p1,(GEN)v[k+1])));

        av1=avma;
	p1 = mpsub(borne1, (GEN)y[k]);
	if (signe(p1) < 0) { avma=av1; fl = 1; }
        else
        {
          p1 = mpadd(eps,mpsub(mpsqrt(gdiv(p1,(GEN)v[k])), (GEN)z[k]));
          avma=av1; x[k] = lmpent(p1); /* safe */
        }
      }
      else if (check)
      { /* don't waste time on the [x,0,...0] */
        for (i=k+1; i<N; i++)
          if (signe(x[i])) break;
        if (i == N && signe(x[1])) x[1] = un;
      }
      for(;;)
      {
        av1=avma;
	if (!fl)
	{
	  p1 = mpmul((GEN)v[k], gsqr(mpadd((GEN)x[k],(GEN)z[k])));
	  if (mpcmp(mpsub(mpadd(p1,(GEN)y[k]), borne1),
	            gmul2n(mpabs(p1),-epsbit)) <= 0) { avma=av1; break; }
	}
        avma=av1; k++; fl=0;
        x[k]=laddis((GEN)x[k],-1);
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
	GEN *gptr[8];
        int a = 4;
	if(DEBUGMEM>1) err(warnmem,"smallvectors");
	gptr[0]=&x; gptr[1]=&y; gptr[2]=&z; gptr[3]=&normax1;
	if (stockmax)
        {
          for (i=s+1; i<=stockmax; i++) S[i]=zero;
          gptr[4]=&S; a++;
        }
        if (check)
        {
          gptr[5]=&borne1; gptr[6]=&borne2; gptr[7]=&norms;
          for (i=s+1; i<=stockmax; i++) norms[i]=zero;
          a+=3;
        }
	gerepilemany(av2,gptr,a);
      }
      if (DEBUGLEVEL>3)
      {
	if (DEBUGLEVEL>5) fprintferr("%ld ",k);
	if (k==n) fprintferr("\nx[%ld] = %Z\n",n,x[n]);
	flusherr();
      }
    }
    while (k>1);
    if (!signe(x[1]) && gexpo((GEN)y[1]) <= -epsbit) break;

    av1=avma; p1 = gsqr(mpadd((GEN)x[1],(GEN)z[1]));
    norme1 = mpadd((GEN)y[1], mpmul(p1, (GEN)v[1]));
    if (mpcmp(norme1,borne1) > 0) avma=av1;
    else
    {
      norme1 = gerepileupto(av1,norme1);
      if (check)
      {
        if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
        {
          if (check(x))
          {
            borne1 = mpadd(norme1,eps);
            borne2 = mpmul(borne1,alpha);
            s = 0; checkcnt = 0;
          }
          else { checkcnt++ ; goto CONTINUE; }
        }
      }
      else if (mpcmp(norme1,normax1) > 0) normax1=norme1;
      if (++s <= stockmax)
      {
        norms[s] = (long)norme1;
        p1 = cgetg(N,t_COL); S[s] = (long)p1;
	for (i=1; i<N; i++) p1[i]=x[i];
        if (s == stockmax && (flag&2) && check)
        {
          long av1 = avma;
          GEN per = sindexsort(norms);
          if (DEBUGLEVEL) fprintferr("sorting...\n");
          for (i=1; i<=s; i++)
          {
            long k = per[i];
            if (check((GEN)S[k]))
            {
              S[1] = S[k]; avma = av1;
              borne1 = mpadd(norme1,eps);
              borne2 = mpmul(borne1,alpha);
              s = 1; checkcnt = 0; break;
            }
          }
          if (checkcnt) s = 0;
        }
      }
    }
CONTINUE:
    x[k] = laddis((GEN)x[k],-1);
  }
  if (s<stockmax) stockmax = s;
  if (check)
  {
    GEN per, alph,pols,p;
    if (DEBUGLEVEL) fprintferr("final sort & check...\n");
    setlg(norms,s+1); per = sindexsort(norms);
    alph = cgetg(s+1,t_VEC);
    pols = cgetg(s+1,t_VEC);
    for (j=0,i=1; i<=s; i++)
    {
      long k = per[i];
      if (j && mpcmp((GEN)norms[k], borne1) > 0) break;
      if ((p = check((GEN)S[k])))
      {
        if (!j) borne1 = gadd((GEN)norms[k],eps);
        j++; pols[j]=(long)p; alph[j]=S[k];
      }
    }
    u = cgetg(3,t_VEC);
    setlg(pols,j+1); u[1] = (long)pols;
    setlg(alph,j+1); u[2] = (long)alph; return u;
  }
  u=cgetg(4,t_VEC);
  u[1]=lstoi(s<<1);
  u[2]=(long)normax1; setlg(S,stockmax+1);
  u[3]=(long)S; return u;
}

/* return T2 norm of the polynomial defining nf */
static GEN 
get_Bnf(GEN nf)
{
  GEN p = gzero, r = (GEN)nf[6];
  long i, r1 = itos(gmael(nf,2,1)), ru = lg(r)-1;
  for (i=ru; i>0; i--)
  {
    if (i == r1) p = gmul2n(p, 1);
    p = gadd(p, gnorm((GEN)r[i]));
  }
  if (i == r1) p = gmul2n(p, 1);
  return p;
}

/* solve x~.a.x <= borne, a > 0. If check is non-NULL keep x only if check(x).
 * flag = 1, return NULL in case of precision problems (sqred1 or lllgram)
 *   raise an error otherwse.
 * flag = 2, return as soon as stockmax vectors found.
 * flag = 3, corresponds to 1+2 */
GEN
fincke_pohst(GEN a,GEN bound,GEN stockmax,long flag, long prec,
             GEN (*check)(GEN))
{
  long pr,av=avma,i,j,n;
  GEN B,nf,r,rinvtrans,v,v1,u,s,res,z,vnorm,sperm,perm,uperm,basis,gram;

  if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); }
  if (typ(a) == t_VEC) { nf = a; a = gmael(nf,5,3); } else nf = NULL;
  pr = gprecision(a);
  if (pr) prec = pr; else a = gmul(a,realun(prec));
  if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
  v1 = lllgramintern(a,4,flag&1, (prec<<1)-2);
  if (v1 == NULL) goto PRECPB;
  r = qf_base_change(a,v1,1);
  r = sqred1intern(r,flag&1);
  if (r == NULL) goto PRECPB;

  n = lg(a);
  for (i=1; i<n; i++)
  {
    GEN p1 = gsqrt(gcoeff(r,i,i), prec);
    coeff(r,i,i)=(long)p1;
    for (j=i+1; j<n; j++)
      coeff(r,i,j) = lmpmul(p1, gcoeff(r,i,j));
  }
  /* now r~ * r = a in LLL basis */
  rinvtrans = gtrans(invmat(r));
  if (DEBUGLEVEL>2)
    fprintferr("final LLL: prec = %ld, precision(rinvtrans) = %ld\n",
                prec,gprecision(rinvtrans));
  v = lllintern(rinvtrans,flag&1, (gprecision(rinvtrans)<<1)-2);
  if (v == NULL) goto PRECPB;
  rinvtrans = gmul(rinvtrans,v);

  u = invmat(gtrans(v)); s = gmul(r,u);
  u = gmul(v1, u);
  vnorm=cgetg(n,t_VEC);
  for (j=1; j<n; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
  perm = sindexsort(vnorm);
  sperm = cgetg(n,t_MAT);
  uperm = cgetg(n,t_MAT);
  for (i=1; i<n; i++) { uperm[n-i] = u[perm[i]]; sperm[n-i] = s[perm[i]]; }

  gram = gram_matrix(sperm);
  B = gcoeff(gram,n-1,n-1);
  if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB;

  if (check && nf) basis = init_chk(nf,uperm,NULL);
  if (!bound)
  { /* polred */
    GEN x = cgetg(n,t_COL);
   
    if (nf) bound = get_Bnf(nf);
    for (j=1; j<n; j++) x[j]=zero;
    for (j=2; j<n; j++)
    {
      x[j]=un; B = gcoeff(gram,j,j);
      if (!bound || mpcmp(B,bound) < 0)
        if (!check || check(x)) bound = B;
      x[j]=zero;
    }
  }

  if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();}
  if (check && nf) 
    if (! (prec = (long)init_chk(nf,uperm,bound))) goto PRECPB;
  i = check? 2: 1; if (i == n) i--;
  for (   ; i<n; i++)
  {
    res = smallvectors(gram,bound? bound: gceil(gcoeff(gram,i,i)),
                       stockmax,flag,prec,check);
    if (!res) goto PRECPB;
    if (!check || lg(res[2]) > 1) break;
    if (DEBUGLEVEL>2) fprintferr("  i = %ld failed\n",i);
  }
  if (check)
  {
    GEN p1 = (GEN)res[2];
    for (i=1; i<lg(p1); i++) p1[i] = lmul(basis, (GEN)p1[i]);
    return res;
  }

  if (DEBUGLEVEL>2) {fprintferr("leaving fincke_pohst\n"); flusherr();}
  z=cgetg(4,t_VEC);
  z[1]=lcopy((GEN)res[1]);
  z[2]=pr? lcopy((GEN)res[2]) : lround((GEN)res[2]);
  z[3]=lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
PRECPB:
  if (!(flag & 1))
    err(talker,"not a positive definite form in fincke_pohst");
  avma = av; return NULL;
}