/*******************************************************************/
/*                                                                 */
/*         CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN)          */
/*                    GENERAL NUMBER FIELDS                        */
/*                                                                 */
/*******************************************************************/
/* $Id: buch2.c,v 1.11 1999/09/24 16:19:25 karim Exp $ */
#include "pari.h"
#include "parinf.h"
long addcolumntomatrix(long *V, long n,long r,GEN *INVP,long *L);
double check_bach(double cbach, double B);
GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
GEN get_arch(GEN nf,GEN x,long prec);
GEN get_roots(GEN x,long r1,long ru,long prec);
long ideal_is_zk(GEN ideal,long N);
GEN idealpowred_prime(GEN nf, GEN vp, GEN n, long prec);
long int_elt_val(GEN nf, GEN x, GEN p, GEN b, long v);
GEN make_M(long n,long ru,GEN basis,GEN roo);
GEN make_MC(long n,long r1,long ru,GEN M);
GEN make_TI(GEN nf, GEN TI, GEN con);

#define SFB_MAX 2
#define SFB_STEP 2
#define MIN_EXTRA 1

#define RANDOM_BITS 4
static const int CBUCHG = (1<<RANDOM_BITS) - 1;
static const int randshift = BITS_IN_RANDOM-1 - RANDOM_BITS;
#undef RANDOM_BITS

static long KC,KCZ,KCZ2,MAXRELSUP;
static long primfact[500],expoprimfact[500];
static long *factorbase, *numfactorbase, *numideal;
static GEN *idealbase, vectbase, powsubfb;

/* factorbase[i]     i-th rational prime used in factor base
 * numfactorbase[i]  index k such that factorbase[k]=i (0 if i is not prime)
 *
 * vectbase          vector of all ideals in factorbase
 * vecbase o subfb = part of factorbase used to build random relations
 * powsubfb  array lg(subfb) x (CBUCHG+1) = all powers up to CBUCHG
 *
 * idealbase[i]      prime ideals above i in factorbase
 * numideal[i]       index k such that idealbase[k] = i.
 *
 * matcopy           all relations found (as long integers, not reduced)
 * cmptglob          lg(matcopy) = total number of relations found
 *
 * Use only non-inert primes, coprime to discriminant index F:
 *   KC = number of prime ideals in factor base (norm < Bach cst)
 *   KC2= number of prime ideals assumed to generate class group (>= KC)
 *
 *   KCZ = number of rational primes under ideal counted by KC
 *   KCZ2= same for KC2le nombre d'ideaux premiers utilises au total.
 */

/* x[0] = length(x) */
static long
ccontent(long* x)
{
  long i, s=labs(x[1]);
  for (i=2; i<=x[0] && s>1; i++) s = cgcd(x[i],s);
  return s;
}

static void
desallocate(long **matcopy)
{
  long i;
  free(numfactorbase); free(factorbase); free(numideal); free(idealbase);
  if (matcopy)
  {
    for (i=lg(matcopy)-1; i; i--) free(matcopy[i]);
    free(matcopy); matcopy = NULL;
  }
  powsubfb = NULL;
}

/* Return the list of indexes or the primes chosen for the subfactorbase.
 * Fill vperm (if !=0): primes ideals sorted by increasing norm (except the
 * ones in subfactorbase come first [dense rows come first for hnfspec])
 * ss = number of rational primes whose divisors are all in factorbase
 */
static GEN
subfactorbasegen(long N,long m,long minsfb,GEN vperm, long *ptss)
{
  long av = avma,i,j, lv=lg(vectbase),s=0,s1=0,n=0,ss=0,z=0;
  GEN y1,y2,subfb,perm,perm1,P,Q;
  double prod;

  (void)new_chunk(lv); /* room for subfb */
  y1 = cgetg(lv,t_COL);
  y2 = cgetg(lv,t_COL);
  for (i=1,P=(GEN)vectbase[i];;P=Q)
  { /* we'll sort ideals by norm (excluded ideals = "zero") */
    long e = itos((GEN)P[3]);
    long ef= e*itos((GEN)P[4]);
    
    s1 += ef;
    y2[i] = (long)powgi((GEN)P[1],(GEN)P[4]);
    /* take only unramified ideals */
    if (e>1) { y1[i]=zero; s=0; z++; } else { y1[i]=y2[i]; s += ef; }

    i++; Q = (GEN)vectbase[i];
    if (i == lv || !egalii((GEN)P[1], (GEN)Q[1]))
    { /* don't take all P above a given p (delete the last one) */
      if (s == N) { y1[i-1]=zero; z++; }
      if (s1== N) ss++;
      if (i == lv) break;
      s=0; s1=0;
    }
  }
  if (z+minsfb >= lv) return NULL;

  prod = 1.0;
  perm = sindexsort(y1) + z; /* skip "zeroes" (excluded ideals) */
  for(;;)
  {
    if (++n > minsfb && (z+n >= lv || prod > m + 0.5)) break;
    prod *= gtodouble((GEN)y1[perm[n]]);
  }
  if (prod < m) return NULL;
  n--;

  /* take the first (non excluded) n ideals (wrt norm), put them first, and
   * sort the rest by increasing norm */
  for (j=1; j<=n; j++) y2[perm[j]] = zero;
  perm1 = sindexsort(y2); avma = av;

  subfb = cgetg(n+1, t_VECSMALL);
  if (vperm)
  {
    for (j=1; j<=n; j++) vperm[j] = perm[j];
    for (   ; j<lv; j++) vperm[j] = perm1[j];
  }
  for (j=1; j<=n; j++) subfb[j] = perm[j];

  if (DEBUGLEVEL)
  {
    if (DEBUGLEVEL>3)
    {
      fprintferr("\n***** IDEALS IN FACTORBASE *****\n\n");
      for (i=1; i<=KC; i++) fprintferr("no %ld = %Z\n",i,vectbase[i]);
      fprintferr("\n***** IDEALS IN SUB FACTORBASE *****\n\n");
      P=cgetg(n+1,t_COL);
      for (j=1; j<=n; j++) P[j] = vectbase[subfb[j]];
      outerr(P);
      fprintferr("\n***** INITIAL PERMUTATION *****\n\n");
      fprintferr("vperm = %Z\n\n",vperm);
    }
    msgtimer("subfactorbase (%ld elements)",n);
  }
  *ptss = ss;
  return subfb;
}

static GEN
mulred(GEN nf,GEN x, GEN I, long prec,long precint)
{
  long av = avma;
  GEN p1, y = cgetg(3,t_VEC), z = cgetg(4,t_VEC);

  y[1] = (long)idealmulh(nf,I,(GEN)x[1]);
  y[2] = x[2];
  y = ideallllredall(nf,y,NULL,prec,precint);
  z[3]=(long)dethnf((GEN)y[1]);
  p1 = ideal_two_elt(nf,(GEN)y[1]);
  z[1]=p1[1];
  z[2]=p1[2]; y[1] = (long)z;
  return gerepileupto(av,gcopy(y));
}

/* Compute powers of prime ideals (P^0,...,P^a) in subfactorbase (assume a > 1)
 * powsubfb[j][i] contains P_i^j in LLL form + archimedean part
 */
static void
powsubfbgen(GEN nf,GEN subfb,long a,long prec,long precint)
{
  long i,j,n=lg(subfb),N=lgef(nf[1])-3,RU;
  GEN id, *pow;

  powsubfb = cgetg(n, t_VEC);
  if (DEBUGLEVEL)
  { fprintferr("Computing powers for sub-factor base:\n"); flusherr(); }
  RU=itos(gmael(nf,2,1)); RU = RU + (N-RU)/2;
  id=cgetg(3,t_VEC);
  id[1] = (long)idmat(N);
  id[2] = (long)zerocol(RU); settyp(id[2],t_VEC);

  for (i=1; i<n; i++)
  {
    GEN vp = (GEN)vectbase[subfb[i]];
    GEN z = cgetg(4,t_VEC);
    pow = (GEN*)cgetg(a+1,t_VEC);
    powsubfb[i] = (long)pow; pow[0]=id;
    pow[1]=cgetg(3,t_VEC);
    pow[1][1] = (long)z;
    z[1]=vp[1]; z[2]=vp[2]; z[3]=(long)powgi((GEN)vp[1], (GEN)vp[4]);
    pow[1][2] = id[2];
    vp = prime_to_ideal(nf,vp);
    for (j=2; j<=a; j++)
    {
      pow[j] = mulred(nf,pow[j-1],vp,prec,precint);
      if (DEBUGLEVEL>1) fprintferr(" %ld",j);
    }
    if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
  }
  if (DEBUGLEVEL)
  {
    if (DEBUGLEVEL>7)
    {
      fprintferr("**** POWERS IN SUB-FACTOR BASE ****\n\n");
      for (i=1; i<n; i++)
      {
        pow = (GEN*)powsubfb[i];
	fprintferr("powsubfb[%ld]:\n",i);
	for (j=0; j<=a; j++) fprintferr("^%ld = %Z\n", j,pow[j]);
	fprintferr("\n");
      }
    }
    msgtimer("powsubfbgen");
  }
}

/* Compute factorbase, numfactorbase, idealbase, vectbase, numideal.
 * n2: bound for norm of tested prime ideals (includes be_honest())
 * n : bound for prime ideals used to build relations (Bach cst) ( <= n2 )

 * Return prod_{p<=n2} (1-1/p) / prod_{Norm(P)<=n2} (1-1/Norm(P)),
 * close to residue of zeta_K at 1 = 2^r1 (2pi)^r2 h R / (w D)
 */
static GEN
factorbasegen(GEN nf,long n2,long n)
{
  byteptr delta=diffptr;
  long KC2,i,j,k,p,lon,ip,nor, N = lgef(nf[1])-3;
  GEN p2,p1,NormP,lfun;
  long prim[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0 };

  numfactorbase= (long*)gpmalloc(sizeof(long)*(n2+1));
  factorbase   = (long*)gpmalloc(sizeof(long)*(n2+1));
  numideal     = (long*)gpmalloc(sizeof(long)*(n2+1));
  idealbase    = (GEN *)gpmalloc(sizeof(GEN )*(n2+1));

  lfun=realun(DEFAULTPREC);
  p=*delta++; i=0; ip=0; KC=0;
  while (p<=n2)
  {
    long av = avma, av1;
    if (DEBUGLEVEL>=2) { fprintferr(" %ld",p); flusherr(); }
    prim[2] = p; p1 = primedec(nf,prim); lon=lg(p1);
    av1 = avma;
    divrsz(mulsr(p-1,lfun),p,lfun);
    if (itos(gmael(p1,1,4)) == N) /* p inert */
    {
      NormP = gpowgs(prim,N);
      if (!is_bigint(NormP) && (nor=NormP[2]) <= n2)
	divrsz(mulsr(nor,lfun),nor-1, lfun);
      avma = av1;
    }
    else
    {
      numideal[p]=ip;
      i++; numfactorbase[p]=i; factorbase[i]=p;
      for (k=1; k<lon; k++,ip++)
      {
	NormP = powgi(prim,gmael(p1,k,4));
	if (is_bigint(NormP) || (nor=NormP[2]) > n2) break;

        divrsz(mulsr(nor,lfun),nor-1, lfun);
      }
      /* keep all ideals with Norm <= n2 */
      avma = av1;
      if (k == lon)
        setisclone(p1); /* flag it: all prime divisors in factorbase */
      else
        { setlg(p1,k); p1 = gerepile(av,av1,gcopy(p1)); }
      idealbase[i] = p1;
    }
    if (!*delta) err(primer1);
    p += *delta++;
    if (KC == 0 && p>n) { KCZ=i; KC=ip; }
  }
  if (!KC) return NULL;
  KCZ2=i; KC2=ip; MAXRELSUP = min(50,4*KC);

  vectbase=cgetg(KC+1,t_COL);
  for (i=1; i<=KCZ; i++)
  {
    p1 = idealbase[i]; k=lg(p1);
    p2 = vectbase + numideal[factorbase[i]];
    for (j=1; j<k; j++) p2[j]=p1[j];
  }
  if (DEBUGLEVEL)
  {
    if (DEBUGLEVEL>1) fprintferr("\n");
    if (DEBUGLEVEL>6)
    {
      fprintferr("########## FACTORBASE ##########\n\n");
      fprintferr("KC2=%ld, KC=%ld, KCZ=%ld, KCZ2=%ld, MAXRELSUP=%ld\n",
                  KC2, KC, KCZ, KCZ2, MAXRELSUP);
      for (i=1; i<=KCZ; i++)
	fprintferr("++ idealbase[%ld] = %Z",i,idealbase[i]);
    }
    msgtimer("factor base");
  }
  return lfun;
}

/* can we factor I / m ? (m pseudo minimum, computed in ideallllredpart1) */
static long
factorisegen(GEN nf,GEN idealvec,long kcz,long limp)
{
  long i,j,n1,ip,v,p,k,lo,ifinal;
  GEN x,q,r,P,p1,listexpo;
  GEN I = (GEN)idealvec[1];
  GEN m = (GEN)idealvec[2];
  GEN Nm= (GEN)idealvec[3];

  x = divii(Nm, dethnf_i(I)); /* m in I, so NI | Nm */
  if (is_pm1(x)) { primfact[0]=0; return 1; }
  listexpo = new_chunk(kcz+1);
  for (i=1; ; i++)
  {
    p=factorbase[i]; q=dvmdis(x,p,&r);
    for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
    listexpo[i] = k;
    if (cmpis(q,p)<=0) break;
    if (i==kcz) return 0;
  }
  if (cmpis(x,limp) > 0) return 0;

  ifinal = i; lo = 0;
  for (i=1; i<=ifinal; i++)
  {
    k = listexpo[i];
    if (k)
    {
      p = factorbase[i]; p1 = idealbase[numfactorbase[p]];
      n1 = lg(p1); ip = numideal[p];
      for (j=1; j<n1; j++)
      {
        P = (GEN)p1[j];
	v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);
	if (v) /* hence < 0 */
	{
	  primfact[++lo]=ip+j; expoprimfact[lo]=v;
	  k += v * itos((GEN)P[4]);
	  if (!k) break;
	}
      }
      if (k) return 0;
    }
  }
  if (is_pm1(x)) { primfact[0]=lo; return 1; }

  p = itos(x); p1 = idealbase[numfactorbase[p]];
  n1 = lg(p1); ip = numideal[p];
  for (k=1,j=1; j<n1; j++)
  {
    P = (GEN)p1[j];
    v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);
    if (v)
    {
      primfact[++lo]=ip+j; expoprimfact[lo]=v;
      k += v*itos((GEN)P[4]);
      if (!k) { primfact[0]=lo; return 1; }
    }
  }
  return 0;
}

/* can we factor alpha ? */
static long
factorisealpha(GEN nf,GEN alpha,long kcz,long limp)
{
  long i,j,n1,ip,v,p,k,lo,ifinal;
  GEN x,q,r,P,p1,listexpo;

  x = absi(subres(gmul((GEN)nf[7],alpha), (GEN)nf[1]));
  if (is_pm1(x)) { primfact[0]=0; return 1; }
  listexpo = new_chunk(kcz+1);
  for (i=1; ; i++)
  {
    p=factorbase[i]; q=dvmdis(x,p,&r);
    for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
    listexpo[i] = k;
    if (cmpis(q,p)<=0) break;
    if (i==kcz) return 0;
  }
  if (cmpis(x,limp) > 0) return 0;

  ifinal=i; lo = 0;
  for (i=1; i<=ifinal; i++)
  {
    k = listexpo[i];
    if (k)
    {
      p = factorbase[i]; p1 = idealbase[numfactorbase[p]];
      n1 = lg(p1); ip = numideal[p];
      for (j=1; j<n1; j++)
      {
        P = (GEN)p1[j];
	v = int_elt_val(nf,alpha,(GEN)P[1],(GEN)P[5], k);
	if (v)
	{
	  primfact[++lo]=ip+j; expoprimfact[lo]=v;
	  k -= v * itos((GEN)P[4]);
	  if (!k) break;
	}
      }
      if (k) return 0;
    }
  }
  if (is_pm1(x)) { primfact[0]=lo; return 1; }

  p = itos(x); p1 = idealbase[numfactorbase[p]];
  n1 = lg(p1); ip = numideal[p];
  for (k=1,j=1; j<n1; j++)
  {
    P = (GEN)p1[j];
    v = int_elt_val(nf,alpha,(GEN)P[1],(GEN)P[5], k);
    if (v)
    {
      primfact[++lo]=ip+j; expoprimfact[lo]=v;
      k -= v*itos((GEN)P[4]);
      if (!k) { primfact[0]=lo; return 1; }
    }
  }
  return 0;
}

static GEN
cleancol(GEN x,long N,long PRECREG)
{
  long i,j,av,tetpil,tx=typ(x),R1,RU;
  GEN s,s2,re,p2,im,y;

  if (tx==t_MAT)
  {
    y=cgetg(lg(x),tx);
    for (j=1; j<lg(x); j++)
      y[j]=(long)cleancol((GEN)x[j],N,PRECREG);
    return y;
  }
  if (!is_vec_t(tx)) err(talker,"not a vector/matrix in cleancol");
  av = avma; RU=lg(x)-1; R1 = (RU<<1)-N;
  re=greal(x); s=(GEN)re[1]; for (i=2; i<=RU; i++) s=gadd(s,(GEN)re[i]);
  s=gdivgs(s,-N); if (N>R1) s2=gmul2n(s,1);
  p2=gmul2n(mppi(PRECREG),2); im=gimag(x);
  tetpil=avma; y=cgetg(RU+1,tx);
  for (i=1; i<=RU; i++)
  {
    GEN p1=cgetg(3,t_COMPLEX); y[i]=(long)p1;
    p1[1] = ladd((GEN)re[i], (i<=R1)?s:s2);
    p1[2] = lmod((GEN)im[i], p2);
  }
  return gerepile(av,tetpil,y);
}

#define RELAT 0
#define LARGE 1
#define PRECI 2
static GEN
not_given(long av, long flun, long reason)
{
  if (labs(flun)==2)
  {
    char *s=NULL;
    switch(reason)
    {
    case RELAT:
      s = "not enough relations for fundamental units, not given"; break;
    case LARGE:
      s = "fundamental units too large, not given"; break;
    case PRECI:
      s = "insufficient precision for fundamental units, not given"; break;
    }
    err(warner,s);
  }
  avma=av; return cgetg(1,t_MAT);
}

/* to check whether the exponential will get too big */
static long
expgexpo(GEN x)
{
  long i,j,e, E = -HIGHEXPOBIT;
  GEN p1;

  for (i=1; i<lg(x); i++)
    for (j=1; j<lg(x[1]); j++)
    {
      p1 = gmael(x,i,j);
      if (typ(p1)==t_COMPLEX) p1 = (GEN)p1[1];
      e = gexpo(p1); if (e>E) E=e;
    }
  return E;
}

static GEN
getfu(GEN nf,GEN *ptxarch,GEN reg,long flun,long *pte,long PRECREG)
{
  long av=avma,i,j,RU,N=lgef(nf[1])-3,e,R1,R2;
  GEN pol,p1,p2,p3,y,matep,s,xarch,vec;
  GEN *gptr[2];

  if (DEBUGLEVEL)
    { fprintferr("\n#### Computing fundamental units\n"); flusherr(); }
  R1=itos(gmael(nf,2,1)); R2=(N-R1)>>1; RU=R1+R2;
  if (RU==1) { *pte=BIGINT; return cgetg(1,t_MAT); }

  *pte = 0; xarch=*ptxarch;
  if (gexpo(reg)<-8) return not_given(av,flun,RELAT);

  matep=cgetg(RU,t_MAT);
  for (j=1; j<RU; j++)
  {
    s=gzero; for (i=1; i<=RU; i++) s=gadd(s,greal(gcoeff(xarch,i,j)));
    s=gdivgs(s,N);
    p1=cgetg(N+1,t_COL); matep[j]=(long)p1;
    for (i=1; i<=R1; i++)
      p1[i]=lsub(gcoeff(xarch,i,j),s);
    for (i=R1+1; i<=RU; i++)
    {
      p1[i]=lsub(gmul2n(gcoeff(xarch,i,j),-1),s);
      p1[i+R2]=lconj((GEN)p1[i]);
    }
  }
  p1 = lllintern(greal(matep),1,PRECREG);
  if (!p1) return not_given(av,flun,PRECI);
  p2 = gmul(matep,p1);
  if (expgexpo(p2) > 20) return not_given(av,flun,LARGE);
  matep=gexp(p2,PRECREG);
  xarch=gmul(xarch,p1);

  p1=gmael(nf,5,1);
  p2=cgetg(N+1,t_MAT);
  for (j=1; j<=N; j++)
  {
    p3=cgetg(N+1,t_COL); p2[j]=(long)p3;
    for (i=1; i<=R1; i++) p3[i]=coeff(p1,i,j);
    for (   ; i<=RU; i++)
    {
      p3[i]=coeff(p1,i,j);
      p3[i+R2]=lconj((GEN)p3[i]);
    }
  }
  y=greal(grndtoi(gauss(p2,matep),&e));
  if (e>=0) return not_given(av,flun,PRECI);
  *pte = -e; pol = (GEN) nf[1];
  p1 = cgetg(3,t_COMPLEX);
  p1[1] = zero; p1[2] = lmppi(PRECREG);  /* p1 = i * pi */
  if (R1<RU) p2 = gshift(p1,1);
  vec = cgetg(RU+1,t_COL);
  for (i=1; i<=R1; i++) vec[i]=(long)p1;
  for (   ; i<=RU; i++) vec[i]=(long)p2;
  p3=cgetg(N+1,t_COL);

  for (j=1; j<lg(y); j++)
  {
    p1=(GEN)y[j]; p2=ginvmod(gmul((GEN)nf[7],p1), pol);
    for (i=1; i<lgef(p2)-1; i++) p3[i]=p2[i+1];
    for (   ; i<=N; i++) p3[i]=zero;
    p2=gmul((GEN)nf[8],p3);
    if (gcmp(gnorml2(p2),gnorml2(p1))<0)
    {
      p1=p2; xarch[j]=lneg((GEN)xarch[j]);
    }
    i=N; while (i>=1 && gcmp0((GEN)p1[i])) i--;
    if (gsigne((GEN)p1[i])>=0) y[j]=(long)p1;
    else
    {
      y[j]=lneg(p1);
      xarch[j]=ladd((GEN)xarch[j],vec);
    }
  }
  p1=gmul((GEN)nf[7],y);
  for (j=1; j<lg(y); j++)
    if (!gcmp1(gabs(gnorm(gmodulcp((GEN)p1[j],pol)),0)))
      { *pte = 0; return not_given(av,flun,LARGE); }
  if (DEBUGLEVEL) msgtimer("getfu");
  *ptxarch=xarch; gptr[0]=ptxarch; gptr[1]=&y;
  gerepilemany(av,gptr,2); return y;
}
#undef RELAT
#undef LARGE
#undef PRECI

GEN
buchfu(GEN bnf)
{
  GEN nf,xarch,reg,res,fu,y;
  long av=avma,tetpil,c,RU;

  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
  RU=itos(gmael(nf,2,1))+itos(gmael(nf,2,2));
  res=(GEN)bnf[8];
  if (lg(res)==7 && lg(res[5])==RU)
  {
    y=cgetg(3,t_VEC); y[1]=lcopy((GEN)res[5]);
    y[2]=lcopy((GEN)res[6]); return y;
  }

  xarch=(GEN)bnf[3]; reg=(GEN)res[2];
  fu=getfu(nf,&xarch,reg,2,&c,gprecision(xarch));
  tetpil=avma; y=cgetg(3,t_VEC);
  y[1]=c?lmul((GEN)nf[7],fu):lcopy(fu); y[2]=lstoi(c);
  return gerepile(av,tetpil,y);
}

static long
factorgensimple(GEN nf,GEN ideal)
{
  long i,v,av1 = avma,lo;
  GEN x = dethnf_i(ideal);

  if (gcmp1(x)) { avma=av1; primfact[0]=0; return 1; }
  for (lo=0, i=1; i<lg(vectbase); i++)
  {
    GEN p1=(GEN)vectbase[i], p=(GEN)p1[1];
    if (!smodis(x,itos(p))) /* if p | x */
    {
      v=idealval(nf,ideal,p1);
      if (v)
      {
	lo++; primfact[lo]=i; expoprimfact[lo]=v;
	x = divii(x, gpuigs(p, v * itos((GEN)p1[4])));
	if (gcmp1(x)) { avma=av1; primfact[0]=lo; return 1; }
      }
    }
  }
  avma=av1; primfact[0]=lo; return 0;
}

static void
add_to_fact(long l, long v, long e)
{
  long i,lo;
  if (!e) return;
  for (i=1; i<=l && primfact[i] < v; i++)/*empty*/;
  if (i <= l && primfact[i] == v)
    expoprimfact[i] += e;
  else
  {
    lo = ++primfact[0];
    primfact[lo] = v;
    expoprimfact[lo] = e;
  }
}

/* factor x on vectbase (modulo principal ideals) */
static GEN
split_ideal(GEN nf, GEN x, GEN xar, long prec, GEN vperm)
{
  GEN id,vdir,x0,y,p1;
  long v1,v2,nbtest,bou,i, ru = lg(xar);
  int flag = (gexpo(gcoeff(x,1,1)) < 100);

  if (flag && factorgensimple(nf,x)) return xar;

  x0 = cgetg(3,t_VEC);
  x = gmul(x, lllintpartial(x));
  x0[1]=(long)x; x0[2]=(long)xar;
  y = ideallllred(nf,x0,NULL,prec);
  if (gcmp0((GEN)y[2])) flag = !flag;
  else
  {
    flag = 1; x=(GEN)y[1];
    x = hnfmod(x, detint(x));
  }
  if (flag && factorgensimple(nf,x)) return (GEN)y[2];

  vdir = cgetg(ru,t_VEC);
  for (i=2; i<ru; i++) vdir[i]=zero;
  for (i=1; i<ru; i++)
  {
    vdir[i]=lstoi(10);
    y = ideallllred(nf,x0,vdir,prec); x=(GEN)y[1];
    if (factorgensimple(nf,x)) return (GEN)y[2];
    vdir[i]=zero;
  }
  v1=itos((GEN)vperm[1]);
  v2=itos((GEN)vperm[2]);
  for(nbtest = 0;;)
  {
    long ex1 = mymyrand() >> randshift;
    long ex2 = mymyrand() >> randshift;
    id=idealpowred_prime(nf,(GEN)vectbase[v1],stoi(ex1),prec);
    p1=idealpowred_prime(nf,(GEN)vectbase[v2],stoi(ex2),prec);
    id = idealmulh(nf,idealmul(nf,x0,id),p1);
    for (i=1; i<ru; i++) vdir[i] = lstoi(mymyrand() >> randshift);
    for (bou=1; bou<ru; bou++)
    {
      if (bou>1)
      {
        for (i=1; i<ru; i++) vdir[i]=zero;
        vdir[bou]=lstoi(10);
      }
      nbtest++;
      y = ideallllred(nf,id,vdir,prec); x=(GEN)y[1];
      if (DEBUGLEVEL>2)
        fprintferr("nbtest = %ld, ideal = %Z\n",nbtest,(long)x);
      if (factorgensimple(nf,x))
      {
        long l = primfact[0];
        add_to_fact(l,v1,-ex1);
        add_to_fact(l,v2,-ex2); return (GEN)y[2];
      }
    }
  }
}

static void
get_split_expo(GEN xalpha, GEN yalpha, GEN vperm)
{
  long i, colW = lg(xalpha)-1;
  GEN vinvperm = new_chunk(lg(vectbase));
  for (i=1; i<lg(vectbase); i++) vinvperm[itos((GEN)vperm[i])]=i;
  for (i=1; i<=primfact[0]; i++)
  {
    long k = vinvperm[primfact[i]];
    long l = k - colW;
    if (l <= 0) xalpha[k]=lstoi(expoprimfact[i]);
    else        yalpha[l]=lstoi(expoprimfact[i]);
  }
}

static GEN
isprincipalall0(GEN bnf, GEN x, long prec, long flall)
{
  long i,j,colW,colB,N,R1,R2,RU,e,c;
  GEN xalpha,yalpha,u2,y,p1,p2,p3,p4,xar,gen,cyc,u1inv,xc,ex;
  GEN W       = (GEN)bnf[1];
  GEN B       = (GEN)bnf[2];
  GEN matunit = (GEN)bnf[3];
  GEN vperm   = (GEN)bnf[6];
  GEN nf      = (GEN)bnf[7];
  GEN RES     = (GEN)bnf[8];
  GEN clg2    = (GEN)bnf[9];

  vectbase = (GEN)bnf[5]; /* needed by factorgensimple */

  N=lgef(nf[1])-3;
  R1=itos(gmael(nf,2,1)); R2=(N-R1)>>1; RU=R1+R2;
  xc = content(x); if (!gcmp1(xc)) x=gdiv(x,xc);

  colW=lg(W)-1; colB=lg(B)-1;
  xar=cgetg(RU+1,t_VEC); for (i=1; i<=RU; i++) xar[i]=zero;
  p1 = split_ideal(nf,x,xar,prec,vperm);
  if (p1 != xar) xar = cleancol(p1,N,prec);

  xalpha=cgetg(colW+1,t_COL); for (i=1; i<=colW; i++) xalpha[i]=zero;
  yalpha=cgetg(colB+1,t_COL); for (i=1; i<=colB; i++) yalpha[i]=zero;
  get_split_expo(xalpha,yalpha,vperm);

  u1inv= (GEN)clg2[1]; /* 1/u1, we have u1*W*u2=diag(cyc_i) */
  u2   = (GEN)clg2[2];
  cyc = gmael(RES,1,2);
  gen = gmael(RES,1,3);

  p1 = gauss(u1inv, gsub(xalpha, gmul(B,yalpha)));
  p4 = cgetg(colW+colB+1,t_COL); c = lg(cyc)-1;
  ex = cgetg(c+1,t_COL);
  for (i=1; i<=c; i++)
    p4[i] = (long)truedvmdii((GEN)p1[i],(GEN)cyc[i],(GEN*)(ex+i));
  if (!(flall & nf_GEN)) return gcopy(ex);

{
  GEN col, baseclorig = (GEN)clg2[3];
  GEN M=gmael(nf,5,1), M2,s,s2;
  GEN Bc = dummycopy((GEN)bnf[4]);

  for (; i<=colW; i++) p4[i]=p1[i];
  for (; i<=colW+colB; i++) p4[i]=yalpha[i-colW];
  p2=cgetg(colW+1,t_MAT);
  for (i=1; i<=colW; i++) p2[i]=Bc[i];
  p3=gmul(p2,u2);
  for (i=1; i<=colW; i++) Bc[i]=p3[i];
  settyp(xar,t_COL); col=gsub(gmul(Bc,p4),xar);
  p4=cgetg(c+1,t_MAT);
  for (j=1; j<=c; j++)
  {
    GEN dmin,pmin,d;
    pmin = p2 = (GEN)baseclorig[j];
    dmin = dethnf((GEN)p2[1]);
    p1 = idealinv(nf,p2);
    p1[1]=(long)numer((GEN)p1[1]);

    d=dethnf((GEN)p1[1]);
    if (mpcmp(d,dmin) < 0) { pmin=p1; dmin=d; }
    p1 = ideallllred(nf,p1,NULL,prec);
    d = dethnf((GEN)p1[1]);
    if (mpcmp(d,dmin) < 0) pmin = p1;

    if (!gegal((GEN)pmin[1], (GEN)gen[j]))
      err(bugparier,"isprincipal(1)");
    p4[j]=lneg((GEN)pmin[2]); settyp(p4[j],t_COL);
  }
  if (c) col = gadd(col,gmul(p4,ex));
  col = cleancol(col,N,prec);

  if (RU > 1)
  {
    s=gzero; p4=cgetg(RU+1,t_MAT);
    for (j=1; j<RU; j++)
    {
      p2=cgetg(RU+1,t_COL); p4[j]=(long)p2;
      p1=gzero;
      for (i=1; i<RU; i++)
      {
        p2[i] = lreal(gcoeff(matunit,i,j));
        p1 = gadd(p1, gsqr((GEN)p2[i]));
      }
      p2[RU]=zero; if (gcmp(p1,s)>0) s=p1;
    }
    p2=cgetg(RU+1,t_COL); p4[RU]=(long)p2;
    for (i=1; i<RU; i++) p2[i]=lreal((GEN)col[i]);
    s=gsqrt(gmul2n(s,RU+1),prec);
    if (gcmpgs(s,100000000)<0) s=stoi(100000000);
    p2[RU]=(long)s;

    p4=(GEN)lll(p4,prec)[RU];
    if (signe(p4[RU]) < 0) p4 = gneg_i(p4);
    if (!gcmp1((GEN)p4[RU])) err(bugparier,"isprincipal(2)");
    setlg(p4,RU);
    col = gadd(col, gmul(matunit,p4));
  }

  s2 = gun;
  for (j=1; j<=c; j++)
    if (signe(ex[j]))
      s2 = mulii(s2, powgi(dethnf_i((GEN)gen[j]), (GEN)ex[j]));
  s = gdivgs(glog(gdiv(dethnf_i(x),s2),prec), N);
  p4 = cgetg(N+1,t_COL);
  for (i=1; i<=R1; i++) p4[i]=lexp(gadd(s,(GEN)col[i]),prec);
  for (   ; i<=RU; i++)
  {
    p4[i]=lexp(gadd(s,gmul2n((GEN)col[i],-1)),prec); ;
    p4[i+R2]=lconj((GEN)p4[i]);
  }
  M2=cgetg(N+1,t_MAT);
  for (j=1; j<=N; j++)
  {
    p1=cgetg(N+1,t_COL); M2[j]=(long)p1;
    for (i=1; i<=R1; i++) p1[i]=coeff(M,i,j);
    for (   ; i<=RU; i++)
    {
      p1[i]=coeff(M,i,j);
      p1[i+R2]=lconj((GEN)p1[i]);
    }
  }
  p1 = gdiv(grndtoi(gmul(s2,greal(gauss(M2,p4))),&e), s2);
  if (e < -5)
  {
    p3 = p1;
    if (!c) p3=idealhermite(nf,p3);
    else
      for (j=1; j<=c; j++)
        p3 = idealmul(nf,p3,idealpow(nf,(GEN)gen[j],(GEN)ex[j]));
    if (!gegal(x,p3)) e=0;
  }
  y=cgetg(4,t_VEC); y[1]=lcopy(ex);
  if (e < -5) { y[2]=lmul(xc,p1); y[3]=lstoi(-e); }
  else
  {
    if (flall & nf_FORCE)
    {
      if (DEBUGLEVEL)
        err(warner,"precision too low for generators, e = %ld",e);
      prec += (e >> TWOPOTBITS_IN_LONG) + (MEDDEFAULTPREC-2);
      return stoi(prec);
    }
    err(warner,"insufficient precision for generators, not given");
    y[2]=lgetg(1,t_COL); y[3]=zero;
  }
}
  return y;
}

static GEN
triv_gen(GEN nf, GEN x, long c, long flag)
{
  GEN y;
  if (!(flag & nf_GEN)) return cgetg(1,t_COL);
  y = cgetg(4,t_VEC);
  y[1] = (long)zerocol(c);
  y[2] = (long)algtobasis(nf,x);
  y[3] = lstoi(BIGINT); return y;
}

GEN
isprincipalall(GEN bnf,GEN x,long flag)
{
  long av = avma,c,pr, tx = typ(x);
  GEN nf;

  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
  if (tx == t_POLMOD)
  {
    if (!gegal((GEN)x[1],(GEN)nf[1]))
      err(talker,"not the same number field in isprincipal");
    x = (GEN)x[2]; tx = t_POL;
  }
  if (tx == t_POL)
  {
    if (gcmp0(x)) err(talker,"zero ideal in isprincipal");
    return triv_gen(nf, x, lg(mael3(bnf,8,1,2))-1, flag);
  }
  x = idealhermite(nf,x);
  if (lg(x)==1) err(talker,"zero ideal in isprincipal");
  if (lgef(nf[1])==4)
    return gerepileupto(av, triv_gen(nf, gcoeff(x,1,1), 0, flag));

  pr = gprecision(gmael(bnf,4,1)); /* precision of unit matrix */
  c = getrand();
  for (;;)
  {
    long av1 = avma;
    GEN y = isprincipalall0(bnf,x,pr,flag);
    if (typ(y) != t_INT) return gerepileupto(av,y);

    pr = itos(y); avma = av1;
    if (DEBUGLEVEL) err(warnprec,"isprincipalall0",pr);
    bnf = bnfnewprec(bnf,pr); setrand(c);
  }
}

GEN
isprincipal(GEN bnf,GEN x)
{
  return isprincipalall(bnf,x,nf_REGULAR);
}

GEN
isprincipalgen(GEN bnf,GEN x)
{
  return isprincipalall(bnf,x,nf_GEN);
}

GEN
isprincipalforce(GEN bnf,GEN x)
{
  return isprincipalall(bnf,x,nf_FORCE);
}

GEN
isprincipalgenforce(GEN bnf,GEN x)
{
  return isprincipalall(bnf,x,nf_GEN | nf_FORCE);
}

GEN
isunit(GEN bnf,GEN x)
{
  long av=avma,tetpil,tx = typ(x),i,R1,RU,n;
  GEN p1,logunit,y,ex,nf,z,pi2_sur_w,gn,emb;

  bnf = checkbnf(bnf); nf=(GEN)bnf[7];
  logunit=(GEN)bnf[3]; RU=lg(logunit);
  p1 = gmael(bnf,8,4); /* roots of 1 */
  gn = (GEN)p1[1]; n = itos(gn);
  z  = (GEN)p1[2];
  switch(tx)
  {
    case t_INT: case t_FRAC: case t_FRACN:
      if (!gcmp1(x) && !gcmp_1(x)) return cgetg(1,t_COL);
      y = zerocol(RU); i = (gsigne(x) > 0)? 0: n>>1;
      y[RU] = (long)gmodulss(i, n); return y;

    case t_POLMOD:
      if (!gegal((GEN)nf[1],(GEN)x[1]))
        err(talker,"not the same number field in isunit");
      x = (GEN)x[2]; /* fall through */
    case t_POL:
      p1 = x; x = algtobasis(bnf,x); break;

    case t_COL:
      if (lgef(nf[1])-2 == lg(x)) { p1 = basistoalg(nf,x); break; }

    default:
      err(talker,"not an algebraic number in isunit");
  }
  if (!gcmp1(denom(x))) { avma = av; return cgetg(1,t_COL); }
  if (typ(p1) != t_POLMOD) p1 = gmodulcp(p1,(GEN)nf[1]);
  p1 = gnorm(p1);
  if (!is_pm1(p1)) { avma = av; return cgetg(1,t_COL); }

  R1 = itos(gmael(nf,2,1)); p1 = cgetg(RU+1,t_COL);
  for (i=1; i<=R1; i++) p1[i] = un;
  for (   ; i<=RU; i++) p1[i] = deux;
  logunit = concatsp(logunit,p1);
  /* ex = fundamental units exponents */
  ex = ground(gauss(greal(logunit), get_arch_real(nf,x,&emb, MEDDEFAULTPREC)));
  if (!gcmp0((GEN)ex[RU]))
    err(talker,"insufficient precision in isunit");

  setlg(ex, RU);
  setlg(p1, RU); settyp(p1, t_VEC);
  for (i=1; i<RU; i++) p1[i] = coeff(logunit, 1, i);
  p1 = gneg(gimag(gmul(p1,ex))); if (!R1) p1 = gmul2n(p1, -1);
  p1 = gadd(garg((GEN)emb[1],DEFAULTPREC), p1);
  /* p1 = arg(the missing root of 1) */

  pi2_sur_w = divrs(mppi(DEFAULTPREC), n>>1);
  p1 = ground(gdiv(p1, pi2_sur_w));
  if (n > 2)
  {
    GEN ro = gmael(nf,6,1);
    GEN p2 = ground(gdiv(garg(poleval(z,ro), DEFAULTPREC), pi2_sur_w));
    p1 = mulii(p1,  mpinvmod(p2, gn));
  }

  tetpil = avma; y = cgetg(RU+1,t_COL);
  for (i=1; i<RU; i++) y[i] = lcopy((GEN)ex[i]);
  y[RU] = lmodulcp(p1, gn); return gerepile(av,tetpil,y);
}

GEN
signunits(GEN bnf)
{
  long av,i,j,R1,RU,mun;
  GEN matunit,y,p1,p2,nf,pi;

  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
  matunit = (GEN)bnf[3]; RU = lg(matunit);
  R1=itos(gmael(nf,2,1)); pi=mppi(MEDDEFAULTPREC);
  y=cgetg(RU,t_MAT); mun = lnegi(gun);
  for (j=1; j<RU; j++)
  {
    p1=cgetg(R1+1,t_COL); y[j]=(long)p1; av=avma;
    for (i=1; i<=R1; i++)
    {
      p2 = ground(gdiv(gimag(gcoeff(matunit,i,j)),pi));
      p1[i] = mpodd(p2)? mun: un;
    }
    avma=av;
  }
  return y;
}

static GEN
quad_form(GEN *cbase,GEN ideal,GEN T2vec,GEN prvec)
{
  long i;
  for (i=1; i<lg(T2vec); i++)
  {
    long prec = prvec[i];
    GEN p1,T2 = (GEN)T2vec[i];

    p1 = qf_base_change(T2,gmul(ideal,realun(prec)), 0);
    if ((*cbase=lllgramintern(p1,100,1,prec)) == NULL)
    {
      if (DEBUGLEVEL) err(warner, "prec too low in quad_form(1): %ld",prec);
      continue;
    }
    if (DEBUGLEVEL>6)
    {
      fprintferr(" input matrix for lllgram: %Z\n",(long)p1);
      fprintferr(" lllgram output (prec = %ld): %Z\n",prec,(long)*cbase);
    }
    p1 = sqred1intern(qf_base_change(p1,*cbase,1),1);
    if (p1) return p1;
    if (DEBUGLEVEL) err(warner, "prec too low in quad_form(2): %ld",prec);
  }
  return NULL;
}

/* y is a vector of LONG, of length ly. x is a hx x ly matrix */
GEN
gmul_mat_smallvec(GEN x, GEN y, long hx, long ly)
{
  GEN z=cgetg(hx+1,t_COL), p1,p2;
  long i,j,av,tetpil;

  for (i=1; i<=hx; i++)
  {
    p1=gzero; av=avma;
    for (j=1; j<=ly; j++)
    {
      p2=gmulgs(gcoeff(x,i,j),y[j]);
      tetpil=avma; p1=gadd(p1,p2);
    }
    z[i]=lpile(av,tetpil,p1);
  }
  return z;
}

static double
get_minkovski(long prec, long N, long R1, GEN D, GEN gborne)
{
  GEN p1,p2, pi = mppi(prec);
  double bound;

  p1 = gsqrt(gmulsg(N,gmul2n(pi,1)),prec);
  p1 = gdiv(p1, gexp(stoi(N),prec));
  p1 = gmulsg(N, gpui(p1, dbltor(2./(double)N),prec));
  p2 = gpui(gdivsg(4,pi), gdivgs(stoi(N-R1),N),prec);
  p1 = gmul(p1,p2);
  bound = gtodouble(gmul(p1, gpui(absi(D), dbltor(1./(double)N),prec)));
  bound = bound*gtodouble(gborne);
  if (DEBUGLEVEL)
  {
    fprintferr("Bound for norms = %.0f\n",bound); flusherr();
  }
  return bound;
}

static void
wr_rel(long *col)
{
  long i;
  fprintferr("\nrel = ");
  for (i=1; i<=KC; i++)
    if (col[i]) fprintferr("%ld^%ld ",i,col[i]);
  fprintferr("\n");
}

static long
small_norm_for_buchall(long t,long **mat,GEN matarch,long nbrel,long LIMC,
		       long PRECREG,GEN nf,GEN gborne,long nbrelpid,GEN invp,
		       long *L)
{
  long av=avma,av1,av2,av3,tetpil,limpile, *x,i,j,k,noideal,ran,keep_old_invp;
  long nbsmallnorm,nbsmallfact,R1,RU, N = lgef(nf[1])-3;
  double *y,*zz,**qq,*vv, MINKOVSKI_BOUND,IDEAL_BOUND,normideal,eps;
  GEN D,V,alpha,T2,ideal,rrr,cbase,T2vec,prvec;

  if (gsigne(gborne)<=0) return t;
  if (DEBUGLEVEL)
  {
    fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel);
    nbsmallnorm = nbsmallfact = 0; flusherr();
  }
  R1=itos(gmael(nf,2,1)); RU = R1 + itos(gmael(nf,2,2));
  D=(GEN)nf[3]; j=N+1; T2=gmael(nf,5,3);
  prvec=cgetg(3,t_VECSMALL);
  prvec[1]=(PRECREG>BIGDEFAULTPREC)? (PRECREG>>1)+1: DEFAULTPREC;
  prvec[2]=PRECREG;
  T2vec=cgetg(3,t_VEC);
  T2vec[1]=(long)gprec_w(T2,prvec[1]);
  T2vec[2]=(long)T2;
  x=(long*)gpmalloc(j*sizeof(long));
  y=(double*)gpmalloc(j*sizeof(double));
  zz=(double*)gpmalloc(j*sizeof(double));
  vv=(double*)gpmalloc(j*sizeof(double));
  qq=(double**)gpmalloc(j*sizeof(double*));
  for (k=1; k<=N; k++) qq[k]=(double*)gpmalloc(j*sizeof(double));

  V=new_chunk(KC+1); av1=avma;
  MINKOVSKI_BOUND = get_minkovski(DEFAULTPREC,N,R1,D,gborne);
  eps = 0.000001;
  for (noideal=1; noideal<=KC; noideal++)
  {
    long flbreak = 0, nbrelideal=0;

    ideal=(GEN)vectbase[KC+1-noideal];
    if (DEBUGLEVEL>1)
    {
      fprintferr("\n*** Ideal no %ld: S = %ld, ",noideal,t);
      fprintferr("prime = %ld, ",itos((GEN)ideal[1]));
      fprintferr("ideal = "); outerr(ideal);
    }
    normideal = gtodouble(powgi((GEN)ideal[1],(GEN)ideal[4]));
    IDEAL_BOUND = MINKOVSKI_BOUND*pow(normideal,2./(double)N);
    ideal = prime_to_ideal(nf,ideal);
    if ((rrr = quad_form(&cbase,ideal,T2vec,prvec)) == NULL)
      return -1; /* precision problem */

    for (k=1; k<=N; k++)
    {
      vv[k]=gtodouble(gcoeff(rrr,k,k));
      for (j=1; j<k; j++) qq[j][k]=gtodouble(gcoeff(rrr,j,k));
      if (DEBUGLEVEL>3) fprintferr("vv[%ld]=%.0f ",k,vv[k]);
    }
    if (DEBUGLEVEL>1)
    {
      if (DEBUGLEVEL>3) fprintferr("\n");
      fprintferr("IDEAL_BOUND = %.0f\n",IDEAL_BOUND); flusherr();
    }
    IDEAL_BOUND += eps; av2=avma; limpile = stack_lim(av2,1);
    x[0]=k=N; y[N]=zz[N]=0; x[N]= (long) sqrt(IDEAL_BOUND/vv[N]);
    for(;; x[1]--)
    {
      for(;;) /* looking for primitive element of small norm */
      {
	double p;

	if (k>1)
	{
	  /* We need to define `l' for NeXTgcc 2.5.8 */
	  long l=k-1;
	  zz[l]=0;
	  for (j=k; j<=N; j++) zz[l] += qq[l][j]*x[j];
	  p=x[k]+zz[k];
	  y[l]=y[k]+p*p*vv[k];
	  x[l]=(long) floor(sqrt((IDEAL_BOUND-y[l])/vv[l])-zz[l]);
          k=l;
	}
	for(;;)
	{
	  p=x[k]+zz[k];
	  if (y[k] + vv[k]*p*p <= IDEAL_BOUND) break;
	  k++; x[k]--;
	}
	if (k==1) /* element complete */
	{
	  if (!x[1] && y[1]<=eps) { flbreak=1; break; }
	  if (ccontent(x)==1) /* primitive */
	  {
	    if (DEBUGLEVEL>4)
            {
              fprintferr("** Found one element: AVMA = %ld\n",avma);
              flusherr();
            }
	    av3=avma; alpha=gmul(ideal,gmul_mat_smallvec(cbase,x,N,N));
	    j=N; while (j>=2 && !signe(alpha[j])) --j;
	    if (j!=1)
	    {
	      if (DEBUGLEVEL)
	      {
		if (DEBUGLEVEL>1)
		{
		  fprintferr(".");
		  if (DEBUGLEVEL>7)
		  {
		    GEN bq = gzero, cq;
		    outerr(gdiv(idealnorm(nf,alpha), idealnorm(nf,ideal)));
		    for (j=1; j<=N; j++)
		    {
		      cq=gzero;
		      for (i=j+1; i<=N; i++)
			cq=gadd(cq,gmulgs(gcoeff(rrr,j,i),x[i]));
		      cq=gaddgs(cq,x[j]);
		      bq=gadd(bq,gmul(gsqr(cq),gcoeff(rrr,j,j)));
		    }
		    outerr(bq);
		  }
		}
		nbsmallnorm++; flusherr();
	      }
	      if (factorisealpha(nf,alpha,KCZ,LIMC)) break; /* can factor it */
	    }
	    avma=av3;
	  }
	  x[1]--;
	}
      }
      if (flbreak) { flbreak=0; break; }

      if (t && t<KC) /* matrix empty or maximal rank */
      {
	for (i=1; i<=KC; i++) V[i]=0;
	for (i=1; i<=primfact[0]; i++) V[primfact[i]] = expoprimfact[i];
	keep_old_invp=0; ran=addcolumntomatrix(V,KC,t,&invp,L);
      }
      else { keep_old_invp=1; ran=t+1; }
      if (ran==t)
	{ if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); } }
      else
      {
	GEN p1, *newcol; /* record the new relation */
        long *colt;

	t=ran; newcol=(GEN*)matarch[t]; colt=mat[t];
        colt[0] = primfact[1]; /* for already_found_relation */
	for (j=1; j<=primfact[0]; j++)
	  colt[primfact[j]] = expoprimfact[j];

	p1=gmul(gmael(nf,5,1),alpha);
	for (j=1; j<=R1; j++)
	  gaffect(glog((GEN)p1[j],PRECREG), newcol[j]);
	for (   ; j<=RU; j++)
	  gaffect(gmul2n(glog((GEN)p1[j],PRECREG),1), newcol[j]);

	if (DEBUGLEVEL)
	{
	  if (DEBUGLEVEL==1) fprintferr("%4ld",t);
	  else
	  {
	    fprintferr("t = %ld. ",t);
	    if (DEBUGLEVEL>2) outerr(alpha);
            wr_rel(colt);
	  }
	  flusherr(); nbsmallfact++;
	}
	if (t>=nbrel) { flbreak=1; break; }
	nbrelideal++; if (nbrelideal==nbrelpid) break;
      }
      if (keep_old_invp)
	avma=av3;
      else if (low_stack(limpile, stack_lim(av2,1)))
      {
	if(DEBUGMEM>1) err(warnmem,"small_norm");
        tetpil=avma; invp=gerepile(av2,tetpil,gcopy(invp));
      }
      if (DEBUGLEVEL>4)
        { fprintferr("** Found one element: AVMA = %ld\n",avma); flusherr(); }
    }
    if (flbreak) break;
    tetpil=avma; invp=gerepile(av1,tetpil,gcopy(invp));
    if (DEBUGLEVEL>1) msgtimer("for this ideal");
  }
  if (DEBUGLEVEL)
  {
    fprintferr("\n");
    msgtimer("small norm relations");
    if (DEBUGLEVEL>1)
    {
      GEN p1,tmp=cgetg(t+1,t_MAT);

      fprintferr("Elements of small norm gave %ld relations.\n",t);
      fprintferr("Computing rank :"); flusherr();
      for(j=1;j<=t;j++)
      {
	p1=cgetg(KC+1,t_COL); tmp[j]=(long)p1;
	for(i=1;i<=KC;i++) p1[i]=lstoi(mat[j][i]);
      }
      tmp = (GEN)indexrank(tmp)[2]; k=lg(tmp)-1;
      fprintferr("rank = %ld; independent columns:\n",k);
      for (i=1; i<=k; i++) fprintferr("%4ld",itos((GEN)tmp[i]));
      fprintferr("\n");
    }
    if(nbsmallnorm)
      fprintferr("\nnb. fact./nb. small norm = %ld/%ld = %f\n",
               nbsmallfact,nbsmallnorm,((double)nbsmallfact)/nbsmallnorm);
    else
      fprintferr("\nnb. small norm = 0\n");
  }
  for (j=1; j<=N; j++) free(qq[j]);
  free(qq); free(x); free(y); free(zz); free(vv);
  avma=av; return t;
}

/* I assumed to be integral HNF */
static GEN
ideallllredpart1(GEN nf, GEN I, GEN matt2, long N, long PRECREGINT)
{
  GEN y,m,idealpro;

  if (!gcmp1(gcoeff(I,N,N))) { y=content(I); if (!gcmp1(y)) I=gdiv(I,y); }
  y = lllgramintern(qf_base_change(matt2,I,1),100,1,PRECREGINT+1);
  if (!y) return NULL;

  /* I, m, y integral */
  m = gmul(I,(GEN)y[1]);
  if (isnfscalar(m)) m = gmul(I,(GEN)y[2]);

  idealpro = cgetg(4,t_VEC);
  idealpro[1] = (long)I;
  idealpro[2] = (long)m; /* elt of small (weighted) T2 norm in I */
  idealpro[3] = labsi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */
  if (DEBUGLEVEL>5) fprintferr("\nidealpro = %Z\n");
  return idealpro;
}

static void
ideallllredpart2(GEN colarch, GEN nf, GEN arch, GEN gamma, long prec)
{
  GEN v = get_arch(nf,gamma,prec);
  long i;
  for (i=lg(v)-1; i; i--)
    gaffect(gadd((GEN)arch[i],gneg((GEN)v[i])), (GEN)colarch[i]);
}

static void
dbg_newrel(long jideal, long jdir, long phase, long cmptglob, long *col,
           GEN colarch, long lim)
{
  fprintferr("\n++++ cmptglob = %ld: new relation (need %ld)", cmptglob, lim);
  wr_rel(col);
  if (DEBUGLEVEL>3)
  {
    fprintferr("(jideal=%ld,jdir=%ld,phase=%ld)", jideal,jdir,phase);
    msgtimer("for this relation");
  }
  if (DEBUGLEVEL>6) fprintferr("archimedian part = %Z\n",colarch);
  flusherr() ;
}

static void
dbg_cancelrel(long i,long jideal,long jdir,long phase, long *col)
{
  fprintferr("rel. cancelled. phase %ld: %ld ",phase,i);
  if (DEBUGLEVEL>3) fprintferr("(jideal=%ld,jdir=%ld)",jideal,jdir);
  wr_rel(col); flusherr();
}

static void
dbg_outrel(long phase,long cmptglob, GEN vperm,long **ma,GEN maarch)
{
  long av,i,j;
  GEN p1,p2;

  if (phase == 0)
  {
    av=avma; p2=cgetg(cmptglob+1,t_MAT);
    for (j=1; j<=cmptglob; j++)
    {
      p1=cgetg(KC+1,t_COL); p2[j]=(long)p1;
      for (i=1; i<=KC; i++) p1[i]=lstoi(ma[j][i]);
    }
    fprintferr("\nRank = %ld, time = %ld\n",rank(p2),timer2());
    if (DEBUGLEVEL>3)
    {
      fprintferr("relations = \n");
      for (j=1; j <= cmptglob; j++) wr_rel(ma[j]);
      fprintferr("\nmatarch = %Z\n",maarch);
    }
    avma=av;
  }
  else if (DEBUGLEVEL>6)
  {
    fprintferr("before hnfadd:\nvectbase[vperm[]] = \n");
    fprintferr("[");
    for (i=1; i<=KC; i++)
    {
      bruterr((GEN)vectbase[vperm[i]],'g',-1);
      if (i<KC) fprintferr(",");
    }
    fprintferr("]~\n");
  }
  flusherr();
}

/* check if we already have a column mat[l] equal to mat[s] */
static long
already_found_relation(long **mat,long s)
{
  long l,bs,cl,*coll,*cols = mat[s];

  bs=1; while (bs<=KC && !cols[bs]) bs++;
  if (bs>KC) return s; /* zero relation */

#if 0
  /* Could check for colinearity and replace by gcd. Useless in practice */
  cs=cols[bs];
  for (l=s-1; l; l--)
  {
    coll=mat[l]; cl=coll[0]; /* = index of first non zero elt in coll */
    if (cl==bs)
    {
      long b=bs;
      cl=coll[cl];
      do b++; while (b<=KC && cl*cols[b] == cs*coll[b]);
      if (b>KC) return l;
    }
  }
#endif
  for (l=s-1; l; l--)
  {
    coll=mat[l]; cl=coll[0]; /* = index of first non zero elt in coll */
    if (cl==bs)
    {
      long b=bs;
      do b++; while (b<=KC && cols[b] == coll[b]);
      if (b>KC) return l;
    }
  }
  cols[0]=bs; return 0;
}

/* if phase != 1 re-initialize static variables. If <0 return immediately */
static long
random_relation(long phase,long cmptglob,long lim,long LIMC,long N,long RU,
                long PRECREG,long PRECREGINT,GEN nf,GEN subfb,GEN lmatt2,
		long **ma,GEN maarch,long *ex,GEN list_jideal)
{
  static long jideal, jdir;
  long i,av,av1,cptzer,nbmatt2,lgsub, jlist = 1, *col;
  GEN colarch,ideal,idealpro,P;

  if (phase != 1) { jideal=jdir=1; if (phase<0) return 0; }
  nbmatt2 = lg(lmatt2)-1;
  lgsub = lg(subfb);
  cptzer = 0;
  if (DEBUGLEVEL && list_jideal)
    fprintferr("looking hard for %Z\n",list_jideal);
  for (av = avma;;)
  {
    if (list_jideal && jlist < lg(list_jideal) && jdir <= nbmatt2) 
      jideal = list_jideal[jlist++];
    if (!list_jideal || jdir <= nbmatt2)
    {
      avma = av;
      P = prime_to_ideal(nf, (GEN)vectbase[jideal]);
    }
    ideal = P;
    do {
      for (i=1; i<lgsub; i++)
      {
        ex[i] = mymyrand()>>randshift;
        if (ex[i])
          ideal = idealmulh(nf,ideal, gmael(powsubfb,i,ex[i]));
      }
    }
    while (typ(ideal)==t_MAT); /* If ex  = 0, try another */
  
    if (phase != 1) jdir = 1; else phase = 2;
    for (av1 = avma; jdir <= nbmatt2; jdir++, avma = av1)
    { /* reduce along various directions */
      if (DEBUGLEVEL>2)
        fprintferr("phase=%ld,jideal=%ld,jdir=%ld,rand=%ld\n",
                   phase,jideal,jdir,getrand());
      idealpro = ideallllredpart1(nf,(GEN)ideal[1], (GEN)lmatt2[jdir],
                                  N, PRECREGINT);
      if (!idealpro) return -2;
      if (!factorisegen(nf,idealpro,KCZ,LIMC))
      {
        if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
        continue;
      }
      /* can factor ideal, record relation */
      col = ma[++cmptglob];
      for (i=1; i<lgsub; i++) col[subfb[i]] = -ex[i];
      for (i=1; i<=primfact[0]; i++) col[primfact[i]] += expoprimfact[i];
      col[jideal]--;
      i = already_found_relation(ma,cmptglob);
      if (i)
      { /* already known. Forget it */
        if (DEBUGLEVEL>1) dbg_cancelrel(i,jideal,jdir,phase,col);
        cmptglob--; for (i=1; i<=KC; i++) col[i]=0;
        if (++cptzer > MAXRELSUP)
        {
          if (list_jideal) { cptzer -= 10; break; }
          return -1;
        }
        continue;
      }

      /* Record archimedian part */
      cptzer=0; colarch = (GEN)maarch[cmptglob];
      ideallllredpart2(colarch,nf,(GEN)ideal[2],(GEN)idealpro[2],PRECREG);
      if (DEBUGLEVEL)
        dbg_newrel(jideal,jdir,phase,cmptglob,col,colarch,lim);

      /* Need more, try next P */
      if (cmptglob < lim) break;

      /* We have found enough. Return */
      if (phase)
      {
        jdir = 1;
        if (jideal == KC) jideal=1; else jideal++;
      }
      else if (DEBUGLEVEL>2)
        fprintferr("Upon exit: jideal=%ld,jdir=%ld\n",jideal,jdir);
      avma = av; return cmptglob;
    }
    if (!list_jideal)
    {
      if (jideal == KC) jideal=1; else jideal++;
    }
  }
}

static long
be_honest(GEN nf,GEN subfb,long RU,long PRECREGINT)
{
  long av,ex,i,j,k,iz,nbtest, N = lgef(nf[1])-3, lgsub = lg(subfb);
  GEN exu=new_chunk(RU+1), MCtw = cgetg(RU+1,t_MAT);
  GEN p1,p2,ideal,idealpro, MC = gmael(nf,5,2), M = gmael(nf,5,1);

  if (DEBUGLEVEL)
  {
    fprintferr("Be honest for primes from %ld to %ld\n",
		factorbase[KCZ+1],factorbase[KCZ2]);
    flusherr();
  }
  av=avma;
  for (iz=KCZ+1; iz<=KCZ2; iz++)
  {
    p1=idealbase[numfactorbase[factorbase[iz]]];
    if (DEBUGLEVEL>1) fprintferr("%ld ", factorbase[iz]);
    for (j=1; j<lg(p1); j++)
      for(nbtest=0;;)
      {
	ideal = prime_to_ideal(nf,(GEN)p1[j]);
	for (i=1; i<lgsub; i++)
	{
	  ex = mymyrand()>>randshift;
	  if (ex) ideal = idealmulh(nf,ideal,gmael3(powsubfb,i,ex,1));
	}
	for (k=1; k<=RU; k++)
	{
	  if (k==1)
            for (i=1; i<=RU; i++) exu[i] = mymyrand()>>randshift;
          else
	  {
	    for (i=1; i<=RU; i++) exu[i] = 0;
            exu[k] = 10;
	  }
          for (i=1; i<=RU; i++)
            MCtw[i] = exu[i]? lmul2n((GEN)MC[i],exu[i]<<1): MC[i];
          p2 = mulmat_real(MCtw,M);
          idealpro = ideallllredpart1(nf,ideal,p2,N,PRECREGINT);
          if (idealpro &&
	      factorisegen(nf,idealpro,iz-1,factorbase[iz-1])) break;
	  nbtest++; if (nbtest==20) return 0;
	}
	avma=av; if (k <= RU) break;
      }
  }
  if (DEBUGLEVEL)
  {
    if (DEBUGLEVEL>1) fprintferr("\n");
    msgtimer("be honest");
  }
  avma=av; return 1;
}

int
trunc_error(GEN x)
{
  return typ(x)==t_REAL && signe(x)
                        && (expo(x)>>TWOPOTBITS_IN_LONG) + 3 > lg(x);
}

/* xarch = complex logarithmic embeddings of units (u_j) found so far */
static GEN
compute_multiple_of_R(GEN xarch,long RU,long N,long PRECREG, GEN *ptsublambda)
{
  GEN v,mdet,Im_mdet,kR,sublambda,lambda,xreal;
  GEN *gptr[2];
  long av = avma, i,j, sreg = lg(xarch)-1, R1 = 2*RU - N;

  if (DEBUGLEVEL) { fprintferr("\n#### Computing regulator\n"); flusherr(); }
  /* xreal = (log |sigma_i(u_j)|) */
  xreal=greal(xarch); v=cgetg(RU+1,t_COL);
  for (i=1; i<=R1; i++) v[i]=un;
  for (   ; i<=RU; i++) v[i]=deux;
  mdet=cgetg(sreg+2,t_MAT); mdet[1]=(long)v;
  for (j=2; j<=sreg+1; j++) mdet[j]=xreal[j-1];
  /* det(Span(mdet)) = N * R */
  Im_mdet = imagereel(mdet,PRECREG);
  if (DEBUGLEVEL) msgtimer("imagereel");

  /* check we have full rank for units */
  if (lg(Im_mdet) != RU+1) { avma=av; return NULL; }
  /* integral multiple of R: the cols we picked form a Q-basis, they have an
   * index in the full lattice */
  kR = gdivgs(det2(Im_mdet), N);
  if (DEBUGLEVEL) msgtimer("detreel");
  /* R > 0.2 uniformly */
  if (gexpo(kR) < -3) { avma=av; return NULL; }

  kR = mpabs(kR);
  sublambda = cgetg(sreg+1,t_MAT);
  lambda = gauss(Im_mdet,xreal); /* rational entries */
  for (i=1; i<=sreg; i++)
  {
    GEN p1 = cgetg(RU,t_COL), p2 = (GEN)lambda[i];
    sublambda[i] = (long)p1;
    for (j=1; j<RU; j++)
    {
      p1[j] = p2[j+1];
      if (trunc_error((GEN)p1[j])) { *ptsublambda = NULL; return gzero; }
    }
  }
  if (DEBUGLEVEL) msgtimer("gauss & lambda");
  *ptsublambda = sublambda;
  gptr[0]=ptsublambda; gptr[1]=&kR;
  gerepilemany(av,gptr,2); return kR;
}

/* Assuming enough relations, c = Rz is close to an even integer, according
 * to Dirichlet's formula. Otherwise, close to a multiple.
 * Compute a tentative regulator (not a multiple this time) */
static GEN
compute_check(GEN sublambda, GEN z, GEN *parch, GEN *reg)
{
  long av = avma, av2, tetpil;
  GEN p1,c,den, R = *reg; /* multiple of regulator */

  if (DEBUGLEVEL) { fprintferr("\n#### Computing check\n"); flusherr(); }
  c = gmul(R,z);
  sublambda = bestappr(sublambda,c); den = denom(sublambda);
  if (gcmp(den,c) > 0)
  {
    if (DEBUGLEVEL) fprintferr("c = %Z\nden = %Z\n",c,den);
    avma=av; return NULL;
  }

  p1 = gmul(sublambda,den); tetpil=avma;
  *parch = lllint(p1);

  av2=avma; p1 = det2(gmul(sublambda,*parch));
  affrr(mpabs(gmul(R,p1)), R); avma=av2;

  if (DEBUGLEVEL) msgtimer("bestappr/regulator");
  *parch = gerepile(av,tetpil,*parch); return gmul(R,z);
}

/* U W V = D, Ui = U^(-1) */
GEN
compute_class_number(GEN W, GEN *D,GEN *Ui,GEN *V)
{
  GEN S = smith2(W);

  if (DEBUGLEVEL) { fprintferr("#### Computing class number\n"); flusherr(); }
  *Ui= ginv((GEN)S[1]);
  *V = (GEN)S[2];
  *D = (GEN)S[3];
  if (DEBUGLEVEL>=4) msgtimer("smith/class group");
  return dethnf_i(*D);
}

static void
class_group_gen(GEN nf,GEN cyc,GEN clh,GEN u1,GEN u2,GEN vperm,
                GEN *ptclg1,GEN *ptclg2, long prec)
{
  GEN basecl,baseclorig,I,J,p1,dmin,d, Vbase = vectbase;
  long i,j,s,inv, lo = lg(cyc), lo0 = lo;

  if (DEBUGLEVEL)
    { fprintferr("#### Computing class group generators\n"); flusherr(); }
  if (vperm)
  {
    s = lg(Vbase); Vbase = cgetg(s,t_VEC);
    for (i=1; i<s; i++) Vbase[i] = vectbase[vperm[i]];
  }
  if (typ(cyc) == t_MAT)
  { /* diagonal matrix */
    p1 = cgetg(lo,t_VEC);
    for (j=1; j<lo; j++)
    {
      p1[j] = coeff(cyc,j,j);
      if (gcmp1((GEN)p1[j])) break;
    }
    lo0 = lo; lo = j;
    cyc = p1; setlg(cyc, lo);
  }
  baseclorig = cgetg(lo,t_VEC); /* generators = Vbase * u1 (LLL-reduced) */
  basecl = cgetg(lo,t_VEC);
  for (j=1; j<lo; j++)
  {
    p1 = gcoeff(u1,1,j);
    I = idealpowred_prime(nf,(GEN)Vbase[1],p1,prec);
    if (signe(p1)<0) I[1] = lmul((GEN)I[1],denom((GEN)I[1]));
    for (i=2; i<lo0; i++)
    {
      p1=gcoeff(u1,i,j); s=signe(p1);
      if (s)
      {
	J = idealpowred_prime(nf,(GEN)Vbase[i],p1,prec);
        if (s<0) J[1] = lmul((GEN)J[1],denom((GEN)J[1]));
	I = idealmulh(nf,I,J);
	I = ideallllred(nf,I,NULL,prec);
      }
    }
    baseclorig[j]=(long)I; I=(GEN)I[1]; /* I = a generator, order cyc[j] */
    dmin = dethnf_i(I); J = idealinv(nf,I);
    J = gmul(J,denom(J));
    d = dethnf_i(J);
    /* check if J = denom * I^(-1) has smaller norm */
    if (cmpii(d,dmin) < 0) { inv=1; I=J; dmin=d; }
    else                   { inv=0; }
    /* try reducing (may _increase_ the norm) */
    J = ideallllred(nf,J,NULL,prec);
    d = dethnf_i(J);
    if (cmpii(d,dmin) < 0) { inv=1; I=J; }
    basecl[j] = (long)I;
    if (inv)
    {
      u1[j] = lneg((GEN)u1[j]);
      u2[j] = lneg((GEN)u2[j]);
    }
  }
  p1 = cgetg(4,t_VEC);
  p1[1]=(long)clh;
  p1[2]=(long)cyc;
  p1[3]=(long)basecl; *ptclg1 = p1;
  /* W*u2 = u1*diag(cyc) */
  p1 = cgetg(4,t_VEC);
  p1[1]=(long)u1;
  p1[2]=(long)u2;
  p1[3]=(long)baseclorig; *ptclg2 = p1;
  if (DEBUGLEVEL) msgtimer("classgroup generators");
}

static GEN
compute_matt2(long RU,GEN nf)
{
  GEN matt2, MCcopy, MCshif, M = gmael(nf,5,1), MC = gmael(nf,5,2);
  long i,j,k,n = min(RU,9), N = n*(n+1)/2, ind = 1;

  MCcopy=cgetg(RU+1,t_MAT); MCshif=cgetg(n+1,t_MAT);
  for (k=1; k<=RU; k++) MCcopy[k]=MC[k];
  for (k=1; k<=n; k++) MCshif[k]=lmul2n((GEN)MC[k],20);
  matt2=cgetg(N+1,t_VEC);
  for (j=1; j<=n; j++)
  {
    MCcopy[j]=MCshif[j];
    for (i=1; i<=j; i++)
    {
      MCcopy[i]=MCshif[i];
      matt2[ind++] = (long)mulmat_real(MCcopy,M);
      MCcopy[i]=MC[i];
    }
    MCcopy[j]=MC[j];
  }
  if (DEBUGLEVEL) msgtimer("weighted T2 matrices");
  return matt2;
}

/* no garbage collecting. destroys y */
static GEN
relationrank_partial(GEN ptinvp, GEN y, long k, long n)
{
  long i,j;
  GEN res=cgetg(n+1,t_MAT), p1;

  for (i=k+1; i<=n; i++) y[i] = ldiv(gneg_i((GEN)y[i]),(GEN)y[k]);
  for (j=1; j<=k; j++)
  {
    p1=cgetg(n+1,t_COL); res[j]=(long)p1;
    for (i=1; i<j; i++) p1[i]=zero;
    for (   ; i<k; i++) p1[i]=coeff(ptinvp,i,j);
    p1[k]=ldiv(gcoeff(ptinvp,k,j),(GEN)y[k]);
    if (j==k)
      for (i=k+1; i<=n; i++)
	p1[i]=lmul((GEN)y[i],gcoeff(ptinvp,k,k));
    else
      for (i=k+1; i<=n; i++)
	p1[i]=ladd(gcoeff(ptinvp,i,j), gmul((GEN)y[i], gcoeff(ptinvp,k,j)));
  }
  for (  ; j<=n; j++) res[j]=ptinvp[j];
  return res;
}

/* Programmes de calcul du rang d'une matrice A de M_{ n,r }(Q) avec rang(A)=
 * r <= n On transforme peu a peu la matrice  I dont les colonnes sont les
 * vecteurs de la base canonique de Q^n en une matrice de changement de base
 * P obtenue en prenant comme base les colonnes de A independantes et des
 * vecteurs de la base canonique. On rend P^(-1), L un vecteur ligne a n
 * composantes valant 0 ou 1 selon que le le vecteur correspondant de P est
 * e_i ou x_i (e_i vecteur de la base canonique, x_i i-eme colonne de A)
 */
static GEN
relationrank(long **mat,long n,long r,long *L)
{
  long av = avma,tetpil,i,j,lim;
  GEN ptinvp,y;

  if (r>n) err(talker,"incorrect matrix in relationrank");
  if (DEBUGLEVEL)
  {
    fprintferr("After trivial relations, cmptglob = %ld\n",r);
    msgtimer("mat & matarch");
  }
  lim=stack_lim(av,1); ptinvp=idmat(n);
  for (i=1; i<=r; i++)
  {
    j=1; y = gmul_mat_smallvec(ptinvp,mat[i],n,n);
    while (j<=n && (gcmp0((GEN)y[j]) || L[j])) j++;
    if (j>n && i==r) err(talker,"not a maximum rank matrix in relationrank");
    ptinvp = relationrank_partial(ptinvp,y,j,n); L[j]=1;
    if (low_stack(lim, stack_lim(av,1)))
    {
      if(DEBUGMEM>1) err(warnmem,"relationrank");
      tetpil=avma; ptinvp=gerepile(av,tetpil,gcopy(ptinvp));
    }
  }
  tetpil=avma; ptinvp=gerepile(av,tetpil,gcopy(ptinvp));
  if (DEBUGLEVEL>1)
    { fprintferr("\nRank of trivial relations matrix: %ld\n",r); flusherr(); }
  return ptinvp;
}

/* Etant donnes une matrice dans M_{ n,r }(Q), de rang maximum r < n, un
 * vecteur colonne V a n lignes, la matrice *INVP et le vecteur ligne *L
 * donnes par le programme relationrank() ci-dessus, on teste si le vecteur V
 * est lineairement independant des colonnes de la matrice; si la reponse est
 * non, on rend le rang de la matrice; si la reponse est oui, on rend le rang
 * de la matrice + 1, on met dans *INVP l'inverse de la nouvelle matrice
 * *INVP et dans *L le nouveau vecteur ligne *L
 */
long
addcolumntomatrix(long *V, long n,long r,GEN *INVP,long *L)
{
  long av = avma,i,k;
  GEN ptinvp,y;

  if (DEBUGLEVEL>4)
  {
    fprintferr("\n*** Entering addcolumntomatrix(). AVMA = %ld\n",avma);
    flusherr();
  }
  ptinvp=*INVP; y=gmul_mat_smallvec(ptinvp,V,n,n);
  if (DEBUGLEVEL>6)
  {
    fprintferr("vector = [\n");
    for (i=1; i<n; i++) fprintferr("%ld,",V[i]);
    fprintferr("%ld]~\n",V[n]); flusherr();
    fprintferr("vector in new basis = \n"); outerr(y);
    fprintferr("base change matrix = \n"); outerr(ptinvp);
    fprintferr("list = [");
    for (i=1; i<=n-1; i++) fprintferr("%ld,",L[i]);
    fprintferr("%ld]\n",L[n]); flusherr();
  }
  k=1; while (k<=n && (gcmp0((GEN)y[k]) || L[k])) k++;
  if (k>n) avma=av;
  else
  {
    *INVP = relationrank_partial(ptinvp,y,k,n);
    L[k]=1; r++;
  }
  if (DEBUGLEVEL>4)
  {
    fprintferr("*** Leaving addcolumntomatrix(). AVMA = %ld\n",avma);
    flusherr();
  }
  return r;
}

/* a usage special: uniquement pour passer du format smallbnf au format bnf
 * Ici, vectbase est deja permute, donc pas de vperm. A l'effet de
 * compute_class_number() suivi de class_group_gen().
 */
static void
classintern(GEN nf,GEN W,GEN *ptcl, GEN *ptcl2)
{
  long prec = (long)nfnewprec(nf,-1);
  GEN met,u1,u2, clh = compute_class_number(W,&met,&u1,&u2);
  class_group_gen(nf,met,clh,u1,u2,NULL,ptcl,ptcl2, prec);
}

static GEN
codeprime(GEN bnf, GEN pr)
{
  long j,av=avma,tetpil;
  GEN p,al,fa,p1;

  p=(GEN)pr[1]; al=(GEN)pr[2]; fa=primedec(bnf,p);
  for (j=1; j<lg(fa); j++)
    if (gegal(al,gmael(fa,j,2)))
    {
      p1=mulsi(lg(al)-1,p); tetpil=avma;
      return gerepile(av,tetpil,addsi(j-1,p1));
    }
  err(talker,"bug in codeprime/smallbuchinit");
  return NULL; /* not reached */
}

static GEN
decodeprime(GEN nf, GEN co)
{
  long n,indi,av=avma,tetpil;
  GEN p,rem,p1;

  n=lg(nf[7])-1; p=dvmdis(co,n,&rem); indi=itos(rem)+1;
  p1=primedec(nf,p); tetpil=avma;
  return gerepile(av,tetpil,gcopy((GEN)p1[indi]));
}

static GEN
makematal(GEN bnf)
{
  GEN W,B,pfb,vp,nf,ma,pr;
  long lm,lma,av=avma,tetpil,j,k;

  if (!gcmp0((GEN)bnf[10])) return (GEN)bnf[10];
  W=(GEN)bnf[1]; B=(GEN)bnf[2];
  pfb=(GEN)bnf[5]; vp=(GEN)bnf[6]; nf=(GEN)bnf[7];
  lm=(lg(W)>1)?lg(W[1])-1:0; lma=lm+lg(B);
  ma=cgetg(lma,t_MAT);
  for (j=1; j<lma; j++)
  {
    GEN ex = (j<=lm)? (GEN)W[j]: (GEN)B[j-lm];
    GEN id = (j<=lm)? gun: (GEN)pfb[itos((GEN)vp[j])];
    for (k=1; k<=lm; k++)
    {
      pr=(GEN)pfb[itos((GEN)vp[k])];
      id=idealmul(nf,id,idealpow(nf,pr,(GEN)ex[k]));
    }
    ma[j]=isprincipalgen(bnf,id)[2];
    if (lg(ma[j])==1)
      err(talker,"bnf not accurate enough to create a sbnf (makematal)");
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(ma));
}
			
GEN
smallbuchinit(GEN pol,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,long minsfb,long prec)
{
  long av=avma,tetpil,k;
  GEN y,bnf,pfb,vp,nf,mas,res,uni,v1,v2,v3;

  if (typ(pol)==t_VEC) bnf = checkbnf(pol);
  else
  {
    bnf=buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsfb,-3,prec);
    if (checkbnf(bnf) != bnf)
    {
      err(warner,"non-monic polynomial. Change of variables discarded");
      bnf = (GEN)bnf[1];
    }
  }
  pfb=(GEN)bnf[5]; vp=(GEN)bnf[6]; nf=(GEN)bnf[7];
  mas=(GEN)nf[5]; res=(GEN)bnf[8]; uni=(GEN)res[5];

  tetpil=avma;
  y=cgetg(13,t_VEC); y[1]=lcopy((GEN)nf[1]); y[2]=lcopy(gmael(nf,2,1));
  y[3]=lcopy((GEN)nf[3]); y[4]=lcopy((GEN)nf[7]);
  y[5]=lcopy((GEN)nf[6]); y[6]=lcopy((GEN)mas[5]);
  y[7]=lcopy((GEN)bnf[1]); y[8]=lcopy((GEN)bnf[2]);
  v1=cgetg(lg(pfb),t_VEC); y[9]=(long)v1;
  for (k=1; k<lg(pfb); k++)
    v1[k]=(long)codeprime(bnf,(GEN)pfb[itos((GEN)vp[k])]);
  v2=cgetg(3,t_VEC); y[10]=(long)v2;
  v2[1]=lcopy(gmael(res,4,1));
  v2[2]=(long)algtobasis(bnf,gmael(res,4,2));
  v3=cgetg(lg(uni),t_VEC); y[11]=(long)v3;
  for (k=1; k<lg(uni); k++)
    v3[k]=(long)algtobasis(bnf,(GEN)uni[k]);
  y[12]=gcmp0((GEN)bnf[10])? (long)makematal(bnf): lcopy((GEN)bnf[10]);
  return gerepile(av,tetpil,y);
}

static GEN
get_regulator(GEN mun,long prec)
{
  long av,tetpil;
  GEN p1;

  if (lg(mun)==1) return gun;
  av=avma; p1 = gtrans(greal(mun));
  setlg(p1,lg(p1)-1); p1 = det(p1);
  tetpil=avma; return gerepile(av,tetpil,gabs(p1,prec));
}

static GEN
get_mun(GEN funits, GEN ro, long ru, long r1, long prec)
{
  long j,k,av=avma,tetpil;
  GEN p1,p2, mun = cgetg(ru,t_MAT);

  for (k=1; k<ru; k++)
  {
    p1=cgetg(ru+1,t_COL); mun[k]=(long)p1;
    for (j=1; j<=ru; j++)
    {
      p2 = glog(poleval((GEN)funits[k],(GEN)ro[j]),prec);
      p1[j]=(j<=r1)? (long)p2: lmul2n(p2,1);
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(mun));
}

static GEN
get_mc(GEN nf, GEN alphs, long prec)
{
  GEN mc,p1,p2,p3,p4, bas = (GEN)nf[7], pol = (GEN)nf[1], ro = (GEN)nf[6];
  long ru = lg(ro), n = lgef(pol)-3, r1 = itos(gmael(nf,2,1));
  long j,k, la = lg(alphs);

  mc = cgetg(la,t_MAT);
  for (k=1; k<la; k++)
  {
    p4 = gmul(bas,(GEN)alphs[k]);
    p3 = gdivgs(glog(gabs(subres(pol,p4),prec),prec), n);
    p1 = cgetg(ru,t_COL); mc[k] = (long)p1;
    for (j=1; j<ru; j++)
    {
      p2 = gsub(glog(poleval(p4,(GEN)ro[j]),prec), p3);
      p1[j]=(j<=r1)? (long) p2: lmul2n(p2,1);
    }
  }
  return mc;
}

static void
my_class_group_gen(GEN bnf, GEN *ptcl, GEN *ptcl2)
{
  GEN nf=(GEN)bnf[7], Vbase=(GEN)bnf[5], vperm=(GEN)bnf[6], *gptr[2];
  long av = avma, i, lv = lg(Vbase);

  vectbase = cgetg(lv, t_VEC);
  for (i=1; i<lv; i++) vectbase[i] = Vbase[itos((GEN)vperm[i])];
  classintern(nf,(GEN)bnf[1],ptcl,ptcl2);
  gptr[0]=ptcl; gptr[1]=ptcl2; gerepilemany(av,gptr,2);
}

GEN
bnfnewprec(GEN bnf, long prec)
{
  GEN nf,ro,res,p1,funits,mun,matal,clgp,clgp2, y = cgetg(11,t_VEC);
  long r1,r2,ru,av;

  bnf = checkbnf(bnf); nf = nfnewprec((GEN)bnf[7],prec);
  if (prec <= 0) return nf;
  r1=itos(gmael(nf,2,1)); r2=itos(gmael(nf,2,2));
  ru = r1+r2;
  res=cgetg(7,t_VEC); p1=(GEN)bnf[8];
  funits = check_units(bnf,"bnfnewprec");
  ro=(GEN)nf[6];
  mun = get_mun(funits,ro,ru,r1,prec);
  res[2]=(long)get_regulator(mun,prec);
  res[3]=lcopy((GEN)p1[3]);
  res[4]=lcopy((GEN)p1[4]);
  res[5]=lcopy((GEN)p1[5]);
  res[6]=lcopy((GEN)p1[6]);

  y[1]=lcopy((GEN)bnf[1]);
  y[2]=lcopy((GEN)bnf[2]);
  y[3]=(long)mun;
  av = avma;
  matal = (GEN)bnf[10];
  if (gcmp0(matal))
  {
    if (DEBUGLEVEL) err(warner,"building matal and completing bnf");
    matal = gclone(makematal(bnf)); bnf[10] = (long)matal;
  }
  avma = av;
  y[4]=lpileupto(av, gcopy(get_mc(nf,matal,prec)));
  y[5]=lcopy((GEN)bnf[5]);
  y[6]=lcopy((GEN)bnf[6]);
  y[7]=(long)nf;
  y[8]=(long)res;
  my_class_group_gen(y,&clgp,&clgp2);
  res[1]=(long)clgp;
  y[9]=(long)clgp2;
  y[10]=(long)matal; return y;
}

GEN
bnfmake(GEN sbnf, long prec)
{
  long av = avma, j,k,n,r1,r2,ru,lpf;
  GEN p1,p2,pol,bas,ro,m,mul,pok,M,MC,T2,mas,T,TI,nf,mun,funits;
  GEN pfc,vp,mc,clgp,clgp2,res,y,W,mata,racu,reg;

  if (typ(sbnf)!=t_VEC || lg(sbnf)!=13)
    err(talker,"incorrect sbnf in bnfmake");
  pol=(GEN)sbnf[1]; bas=(GEN)sbnf[4]; n=lg(bas)-1;
  r1=itos((GEN)sbnf[2]); r2=(n-r1)/2; ru=r1+r2;
  ro=(GEN)sbnf[5];
  if (prec > gprecision(ro)) ro=get_roots(pol,r1,ru,prec);

  m=cgetg(n+1,t_MAT);
  for (k=1; k<=n; k++)
  {
    p1=cgetg(n+1,t_COL); m[k]=(long)p1; p2=(GEN)bas[k];
    for (j=1; j<=n; j++) p1[j]=(long)truecoeff(p2,j-1);
  }
  m=invmat(m);
  mul=cgetg(n*n+1,t_MAT);
  for (k=1; k<=n*n; k++)
  {
    pok = gres(gmul((GEN)bas[(k-1)%n+1], (GEN)bas[(long)((k-1)/n)+1]), pol);
    p1=cgetg(n+1,t_COL); mul[k]=(long)p1;
    for (j=1; j<=n; j++) p1[j]=(long)truecoeff(pok,j-1);
  }
  mul=gmul(m,mul);

  M  = make_M(n,ru,bas,ro);
  MC = make_MC(n,r1,ru,M);
  T2 = mulmat_real(MC,M);
  p1=mulmat_real(gconj(MC),M); T=ground(p1);
  if (gexpo(gnorml2(gsub(p1,T))) > -30)
    err(talker,"insufficient precision in bnfmake");
  TI=gmul((GEN)sbnf[3],invmat(T));

  mas=cgetg(8,t_VEC);
  nf=cgetg(10,t_VEC);
  p1=cgetg(3,t_VEC); p1[1]=lstoi(r1); p1[2]=lstoi(r2);
  nf[1]=sbnf[1]  ; nf[2]=(long)p1;  nf[3]=sbnf[3];
  nf[4]=ldet(m)  ; nf[5]=(long)mas; nf[6]=(long)ro;
  nf[7]=(long)bas; nf[8]=(long)m;   nf[9]=(long)mul;

  mas[1]=(long)M; mas[2]=(long)MC; mas[3]=(long)T2;
  mas[4]=(long)T; mas[5]=sbnf[6];  mas[6]=(long)TI;
  mas[7]=(long)make_TI(nf,TI,gun);

  funits=cgetg(ru,t_VEC); p1 = (GEN)sbnf[11];
  for (k=1; k < lg(p1); k++)
    funits[k] = lmul(bas,(GEN)p1[k]);
  mun = get_mun(funits,ro,ru,r1,prec);

  prec=gprecision(ro); if (prec<DEFAULTPREC) prec=DEFAULTPREC;
  mc = get_mc(nf, (GEN)sbnf[12], prec);

  pfc=(GEN)sbnf[9]; lpf=lg(pfc);
  vectbase=cgetg(lpf,t_COL); vp=cgetg(lpf,t_COL);
  for (j=1; j<lpf; j++)
  {
    vp[j]=lstoi(j);
    vectbase[j]=(long)decodeprime(nf,(GEN)pfc[j]);
  }
  classintern(nf,(GEN)sbnf[7], &clgp, &clgp2); /* uses vectbase */

  reg = get_regulator(mun,prec);
  p1=cgetg(3,t_VEC); racu=(GEN)sbnf[10];
  p1[1]=racu[1]; p1[2]=lmul(bas,(GEN)racu[2]);
  racu=p1;

  res=cgetg(7,t_VEC);
  res[1]=(long)clgp; res[2]=(long)reg;     res[3]=(long)dbltor(1.0);
  res[4]=(long)racu; res[5]=(long)funits;  res[6]=lstoi(1000);

  if (lg(sbnf[7])>1) { W=(GEN)sbnf[7]; mata=(GEN)sbnf[8]; }
  else
  {
    long la = lg(sbnf[12]);
    W=cgetg(1,t_MAT); mata=cgetg(la,t_MAT);
    for (k=1; k<la; k++) mata[k]=lgetg(1,t_COL);
  }
  y=cgetg(11,t_VEC);
  y[1]=(long)W; y[2]=(long)mata;     y[3]=(long)mun;
  y[4]=(long)mc;  y[5]=(long)vectbase; y[6]=(long)vp;
  y[7]=(long)nf;  y[8]=(long)res;      y[9]=(long)clgp2; y[10]=zero;
  return gerepileupto(av,gcopy(y));
}

static GEN
classgroupall(GEN P, GEN data, long flag, long prec)
{
  long court[3],doubl[4];
  long av=avma,flun,lx, minsfb=3,nbrelpid=4;
  GEN bach=doubl,bach2=doubl,RELSUP=court,borne=gun;

  if (!data) lx=1;
  else
  {
    lx = lg(data);
    if (typ(data)!=t_VEC || lx > 7)
      err(talker,"incorrect parameters in classgroup");
  }
  court[0]=evaltyp(t_INT) | evallg(3); affsi(5,court);
  doubl[0]=evaltyp(t_REAL)| evallg(4); affrr(dbltor(0.3),doubl);
  avma=av;
  switch(lx)
  {
    case 7: minsfb  = itos((GEN)data[6]);
    case 6: nbrelpid= itos((GEN)data[5]);
    case 5: borne  = (GEN)data[4];
    case 4: RELSUP = (GEN)data[3];
    case 3: bach2 = (GEN)data[2];
    case 2: bach  = (GEN)data[1];
  }
  switch(flag)
  {
    case 0: flun=-2; break;
    case 1: flun=-3; break;
    case 2: flun=-1; break;
    case 3: return smallbuchinit(P,bach,bach2,RELSUP,borne,nbrelpid,minsfb,prec);
    case 4: flun=2; break;
    case 5: flun=3; break;
    case 6: flun=0; break;
  }
  return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsfb,flun,prec);
}

GEN
bnfclassunit0(GEN P, long flag, GEN data, long prec)
{
  if (typ(P)==t_INT) return quadclassunit0(P,0,data,prec);
  if (flag < 0 || flag > 2) err(flagerr,"bnfclassunit");
  return classgroupall(P,data,flag+4,prec);
}

GEN
bnfinit0(GEN P, long flag, GEN data, long prec)
{
#if 0
  THIS SHOULD BE DONE...

  if (typ(P)==t_INT)
  {
    if (flag<4) err(impl,"specific bnfinit for quadratic fields");
    return quadclassunit0(P,0,data,prec);
  }
#endif
  if (flag < 0 || flag > 3) err(flagerr,"bnfinit");
  return classgroupall(P,data,flag,prec);
}

GEN
classgrouponly(GEN P, GEN data, long prec)
{
  GEN y,z;
  long av=avma,tetpil,i;

  if (typ(P)==t_INT)
  {
    z=quadclassunit0(P,0,data,prec); tetpil=avma;
    y=cgetg(4,t_VEC); for (i=1; i<=3; i++) y[i]=lcopy((GEN)z[i]);
    return gerepile(av,tetpil,y);
  }
  z=(GEN)classgroupall(P,data,6,prec)[1]; tetpil=avma;
  return gerepile(av,tetpil,gcopy((GEN)z[5]));
}

GEN
regulator(GEN P, GEN data, long prec)
{
  GEN z;
  long av=avma,tetpil;

  if (typ(P)==t_INT)
  {
    if (signe(P)>0)
    {
      z=quadclassunit0(P,0,data,prec); tetpil=avma;
      return gerepile(av,tetpil,gcopy((GEN)z[4]));
    }
    return gun;
  }
  z=(GEN)classgroupall(P,data,6,prec)[1]; tetpil=avma;
  return gerepile(av,tetpil,gcopy((GEN)z[6]));
}

#ifdef INLINE
INLINE
#endif
GEN
col_dup(long n, GEN col)
{
   GEN c = (GEN) gpmalloc(sizeof(long)*(n+1));
   memcpy(c,col,(n+1)*sizeof(long));
   return c;
}

#ifdef INLINE
INLINE
#endif
GEN
col_0(long n)
{
   GEN c = (GEN) gpmalloc(sizeof(long)*(n+1));
   long i;
   for (i=1; i<=n; i++) c[i]=0;
   return c;
}

static GEN
buchall_end(GEN nf,GEN CHANGE,long fl,long k, GEN fu, GEN clg1, GEN clg2,
            GEN reg, GEN c_1, GEN zu, GEN W, GEN B,
            GEN xarch, GEN matarch, GEN vectbase, GEN vperm)
{
  long l = labs(fl)>1? 11: fl? 9: 8;
  GEN p1,z, RES = cgetg(11,t_COL);

  setlg(RES,l);
  RES[5]=(long)clg1;
  RES[6]=(long)reg;
  RES[7]=(long)c_1;
  RES[8]=(long)zu;
  RES[9]=(long)fu;
  RES[10]=lstoi(k);
  if (fl>=0)
  {
    RES[1]=nf[1];
    RES[2]=nf[2]; p1=cgetg(3,t_VEC); p1[1]=nf[3]; p1[2]=nf[4];
    RES[3]=(long)p1;
    RES[4]=nf[7];
    z=cgetg(2,t_MAT); z[1]=lcopy(RES); return z;
  }
  z=cgetg(11,t_VEC);
  z[1]=(long)W;
  z[2]=(long)B;
  z[3]=(long)xarch;
  z[4]=(long)matarch;
  z[5]=(long)vectbase;
  z[6]=(long)vperm;
  z[7]=(long)nf; RES+=4; RES[0]=evaltyp(t_VEC) | evallg(l-4);
  z[8]=(long)RES;
  z[9]=(long)clg2;
  z[10]=zero; /* dummy: we MUST have lg(bnf) != lg(nf) */
  if (CHANGE) { p1=cgetg(3,t_VEC); p1[1]=(long)z; p1[2]=(long)CHANGE; z=p1; }
  return gcopy(z);
}

static GEN
buchall_for_degree_one_pol(GEN nf, GEN CHANGE, long flun)
{
  long av = avma, k = EXP220;
  GEN W,B,xarch,matarch,vectbase,vperm;
  GEN fu=cgetg(1,t_VEC), reg=gun, c_1=gun, zu=cgetg(3,t_VEC);
  GEN clg1=cgetg(4,t_VEC), clg2=cgetg(4,t_VEC);

  clg1[1]=un; clg1[2]=clg1[3]=clg2[3]=lgetg(1,t_VEC);
  clg2[1]=clg2[2]=lgetg(1,t_MAT);
  zu[1]=deux; zu[2]=lnegi(gun);
  W=B=xarch=matarch=cgetg(1,t_MAT);
  vectbase=cgetg(1,t_COL); vperm=cgetg(1,t_VEC);

  return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c_1,zu,W,B,xarch,matarch,vectbase,vperm));
}

GEN
buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,
        long minsfb,long flun,long prec)
{
  long av = avma,av0,av1,limpile,i,j,k,ss,cmptglob,lgsub;
  long N,R1,R2,RU,PRECREG,PRECREGINT,KCCO,KCCOPRO,RELSUP;
  long extrarel,nlze,sreg,nrelsup,nreldep,phase,slim,matcopymax;
  long first = 1, sfb_increase = 0, sfb_trials = 0;
  long **mat,**matcopy,*ex;
  double cbach,cbach2,drc,LOGD2,lim,LIMC,LIMC2;
  GEN p1,p2,lmatt2,fu,zu,nf,D,xarch,met,W,reg,lfun,z,clh,vperm,subfb;
  GEN B,C,u1,u2,c1,sublambda,pdep,parch,liste,invp,clg1,clg2;
  GEN CHANGE=NULL, extramat=NULL, extraC=NULL, list_jideal = NULL;

  if (DEBUGLEVEL) timer2();

  if (typ(P)==t_POL) nf = NULL;
  else
  {
    nf=checknf(P); P=(GEN)nf[1];
  }
  if (typ(gRELSUP)!=t_INT) gRELSUP=gtrunc(gRELSUP);
  RELSUP = itos(gRELSUP);
  if (RELSUP<=0) err(talker,"not enough relations in bnfxxx");

  /* Initializations */
  N=lgef(P)-3;
  if (!nf)
  {
    nf=initalgall0(P, flun>=0? nf_REGULAR: nf_DIFFERENT,
                   max(BIGDEFAULTPREC,prec));
    if (lg(nf)==3) /* P was a non-monic polynomial, nfinit changed it */
    {
      CHANGE=(GEN)nf[2]; nf=(GEN)nf[1];
    }
    if (DEBUGLEVEL) msgtimer("initalg");
  }
  if (N<=1) return buchall_for_degree_one_pol(nf,CHANGE,flun);
  zu=rootsof1(nf);
  zu[2] = lmul((GEN)nf[7],(GEN)zu[2]);
  if (DEBUGLEVEL) msgtimer("rootsof1");

  R1=itos(gmael(nf,2,1)); R2=(N-R1)>>1; RU=R1+R2;
  D=(GEN)nf[3]; drc=fabs(gtodouble(D));
  LOGD2=log(drc); LOGD2 = LOGD2*LOGD2;
  lim = exp(-(double)N) * sqrt(2*PI*N*drc) * pow(4/PI,(double)R2);
  if (lim < 3.) lim = 3.;
  cbach = min(12., gtodouble(gcbach)); cbach /= 2;
  cbach2 = gtodouble(gcbach2);
  if (DEBUGLEVEL)
  {
    fprintferr("N = %ld, R1 = %ld, R2 = %ld, RU = %ld\n",N,R1,R2,RU);
    fprintferr("D = %Z\n",D);
  }
  av0 = avma;
  matcopy = NULL;
  powsubfb = NULL;

INCREASEGEN:
  if (first) first = 0; else { desallocate(matcopy); avma = av0; }
  sfb_trials = sfb_increase = 0;
  cbach = check_bach(cbach,12.);
  nreldep = nrelsup = 0;
  LIMC = cbach*LOGD2; if (LIMC < 20.) LIMC = 20.;
  LIMC2=max(3. * N, max(cbach,cbach2)*LOGD2);
  if (LIMC2 < LIMC) LIMC2=LIMC;
  if (DEBUGLEVEL) { fprintferr("LIMC = %.1f, LIMC2 = %.1f\n",LIMC,LIMC2); }

  /* initialize factorbase, [sub]vperm */
  lfun = factorbasegen(nf,(long)LIMC2,(long)LIMC);
  if (!lfun) goto INCREASEGEN;

  vperm = cgetg(lg(vectbase), t_VECSMALL);
  subfb = subfactorbasegen(N,(long)min(lim,LIMC2), minsfb, vperm, &ss);
  if (!subfb) goto INCREASEGEN;
  lgsub = lg(subfb);
  ex = cgetg(lgsub,t_VECSMALL);

  PRECREGINT = DEFAULTPREC
             + ((expi(D)*(lgsub-2)+((N*N)>>2))>>TWOPOTBITS_IN_LONG);
  PRECREG = max(prec+1,PRECREGINT);
  KCCO = KC+RU-1 + max(ss,RELSUP);
  if (DEBUGLEVEL)
  {
    fprintferr("nbrelsup = %ld, ss = %ld, ",RELSUP,ss);
    fprintferr("KCZ = %ld, KC = %ld, KCCO = %ld \n",KCZ,KC,KCCO); flusherr();
  }
  mat=(long**)gpmalloc(sizeof(long*)*(KCCO+1));
  setlg(mat, KCCO+1);
  C = cgetg(KCCO+1,t_MAT);
  cmptglob=0;
  /* trivial relations */
  for (i=1; i<=KCZ; i++)
  {
    GEN P = idealbase[i];
    if (isclone(P))
    { /* all prime divisors in factorbase */
      unsetisclone(P); cmptglob++;
      mat[cmptglob] = p1 = col_0(KC);
      C[cmptglob] = (long)(p2 = cgetg(RU+1,t_COL));
      k = numideal[factorbase[i]];
      p1[0] = k+1; p1 += k; /* for already_found_relation */
      k = lg(P);
      for (j=1; j<k; j++) p1[j] = itos(gmael(P,j,3));
      for (j=1; j<=RU; j++) p2[j] = zero;
    }
  }
  /* initialize for other relations */
  for (i=cmptglob+1; i<=KCCO; i++)
  {
    mat[i] = col_0(KC);
    C[i] = (long) (p1 = cgetg(RU+1,t_COL));
    for (j=1; j<=RU; j++)
    {
      p2=cgetg(3,t_COMPLEX);
      p2[1]=lgetr(PRECREG);
      p2[2]=lgetr(PRECREG); p1[j]=(long)p2;
    }
  }
  av1 = avma; liste = new_chunk(KC+1);
  for (i=1; i<=KC; i++) liste[i]=0;
  invp = cmptglob? relationrank(mat,KC,cmptglob,liste): idmat(KC);

  /* relations through elements of small norm */
  cmptglob = small_norm_for_buchall(cmptglob,mat,C,KCCO,(long)LIMC,
                                    PRECREG,nf,gborne,nbrelpid,invp,liste);
  if (cmptglob < 0)
  {
    for (j=1; j<=KCCO; j++) free(mat[j]); free(mat);
    prec=(PRECREG<<1)-2;
    if (DEBUGLEVEL) err(warnprec,"buchall (small_norm)",prec);
    avma = av0; nf = nfnewprec(nf,prec); av0 = avma;
    cbach /= 2;
    goto INCREASEGEN;
  }
  avma = av1; limpile=stack_lim(av1,1);

  slim = KCCO; phase = 0;
  nlze = matcopymax = 0; /* for lint */
  lmatt2 = NULL;

  /* random relations */
  if (cmptglob == KCCO) /* enough relations, initialize nevertheless */
    ((void(*)(long))random_relation)(-1);
  else
  {
    GEN maarch;
    long **ma;

    if (DEBUGLEVEL)
      { fprintferr("\n#### Looking for random relations\n"); flusherr(); }
  LABELINT:
    if (sfb_increase)
    { /* increase subfactorbase */
      sfb_increase = 0;
      if (++sfb_trials >= SFB_MAX) goto INCREASEGEN;
      subfb = subfactorbasegen(N, (long)min(lim,LIMC2),
                                  lgsub-1+SFB_STEP, NULL, &ss);
      if (!subfb) goto INCREASEGEN;
      if (DEBUGLEVEL) fprintferr("*** Increasing subfactorbase\n");
      powsubfb = NULL;
      nreldep = nrelsup = 0;
      lgsub = lg(subfb);
    }

    if (phase == 0) { ma = mat; maarch = C; }
    else /* reduced the relation matrix at least once */
    {
      extrarel = nlze;
      if (extrarel < MIN_EXTRA) extrarel = MIN_EXTRA;
      slim = cmptglob+extrarel;
      setlg(extraC,extrarel+1);
      setlg(extramat,extrarel+1);
      if (slim > matcopymax)
      {
        matcopy = (long**)gprealloc(matcopy, (2*slim+1) * sizeof(long*),
                                             (matcopymax+1) * sizeof(long*));
        matcopymax = 2 * slim;
      }
      setlg(matcopy,slim+1);
      if (DEBUGLEVEL)
	fprintferr("\n(need %ld more relation%s)\n",
                    extrarel, (extrarel==1)?"":"s");
      for (j=cmptglob+1; j<=slim; j++) matcopy[j] = col_0(KC);
      maarch = extraC - cmptglob; /* start at 0, the others at cmptglob */
      ma = matcopy;
    }
    if (!lmatt2)
    {
      lmatt2 = compute_matt2(RU,nf);
      av1 = avma;
    }
    if (!powsubfb)
    {
      powsubfbgen(nf,subfb,CBUCHG+1,PRECREG,PRECREGINT);
      av1 = avma;
    }
    ss = random_relation(phase,cmptglob,slim,(long)LIMC,N,RU,PRECREG,
                         PRECREGINT,nf,subfb,lmatt2,ma,maarch,ex,list_jideal);
    if (ss < 0) /* could not find relations */
    {
      if (phase == 0) { for (j=1; j<=KCCO; j++) free(mat[j]); free(mat); }
      if (ss != -1)
      { /* precision problems */
        prec=(PRECREG<<1)-2;
	if (DEBUGLEVEL) err(warnprec,"buchall (random_relation)",prec);
	avma = av0; nf = nfnewprec(nf,prec);
        av0 = avma; cbach /= 2;
      }
      goto INCREASEGEN;
    }
    if (DEBUGLEVEL > 2) dbg_outrel(phase,cmptglob,vperm,ma,maarch);
    if (phase)
      for (j=1; j<=extrarel; j++)
      {
	long *c = matcopy[cmptglob+j];
	GEN  *g = (GEN*) extramat[j];
	for (k=1; k<=KC; k++) g[k] = stoi(c[vperm[k]]);
      }
    cmptglob = ss;
  }

  /* reduce relation matrices */
  if (phase == 0) /* never reduced before */
  {
    matcopymax = 10*KCCO + MAXRELSUP;
    matcopy = (long**)gpmalloc(sizeof(long*)*(matcopymax + 1));
    setlg(matcopy, KCCO+1);
    for (j=1; j<=KCCO; j++) matcopy[j] = col_dup(KC,mat[j]);
    W = hnfspec(mat,vperm,&pdep,&B,&C,lgsub-1);
    for (j=1; j<=KCCO; j++) free(mat[j]); free(mat);
    KCCOPRO = KCCO; phase = 1;
   /* keep some room for the extra relation. We will update matrix size when
    * extrarel goes down
    */
    nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
    if (nlze)
    {
      list_jideal = cgetg(nlze+1, t_VECSMALL);
      for (i=1; i<=nlze; i++) list_jideal[i] = vperm[i];
    }
    extrarel = nlze; /* in case the regulator is 0 */
    if (extrarel < MIN_EXTRA) extrarel = MIN_EXTRA;
    extramat =cgetg(extrarel+1,t_MAT);
    extraC=cgetg(extrarel+1,t_MAT);
    for (j=1; j<=extrarel; j++)
    {
      extramat[j]=lgetg(KC+1,t_COL);
      extraC[j]=lgetg(RU+1,t_COL);
      for (i=1; i<=RU; i++)
      {
	p1 = cgetg(3,t_COMPLEX); mael(extraC,j,i)=(long)p1;
	p1[1]=lgetr(PRECREG);
	p1[2]=lgetr(PRECREG);
      }
    }
  }
  else
  {
    list_jideal = NULL;
    if (low_stack(limpile, stack_lim(av1,1)))
    {
      GEN *gptr[6];
      if(DEBUGMEM>1) err(warnmem,"buchall");
      gptr[0]=&W; gptr[1]=&C; gptr[2]=&B; gptr[3]=&pdep;
      gptr[4]=&extramat; gptr[5]=&extraC;
      gerepilemany(av1,gptr,6);
    }
    W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);
    nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
    KCCOPRO += extrarel;
    if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto LABELINT; }
  }
  if (nlze) goto LABELINT; /* dependent rows */

  /* first attempt at regulator for the check */
  sreg = KCCOPRO - (lg(B)-1) - (lg(W)-1); /* = zc (cf hnffinal) */
  xarch = cgetg(sreg+1,t_MAT); /* cols corresponding to units */
  for (j=1; j<=sreg; j++) xarch[j] = C[j];
  reg = compute_multiple_of_R(xarch,RU,N,PRECREG,&sublambda);

  if (!reg)
  { /* not full rank for units */
    if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
    if (++nrelsup > MAXRELSUP) goto INCREASEGEN;
    nlze=MIN_EXTRA; goto LABELINT;
  }
  if (!sublambda)
  { /* anticipate precision problems */
    prec=(PRECREG<<1)-2;
    if (DEBUGLEVEL) err(warnprec,"buchall (bestappr)",prec);
    avma = av0; nf = nfnewprec(nf,prec);
    av0 = avma; cbach /= 2;
    goto INCREASEGEN;
  }

  /* class number */
  if (DEBUGLEVEL) fprintferr("\n");
  clh = compute_class_number(W,&met,&u1,&u2);

  /* check */
  z = mulrr(lfun,gmul(gmul2n(gpuigs(shiftr(mppi(DEFAULTPREC),1),-R2),-R1),
		      gsqrt(absi(D),DEFAULTPREC)));
  z = mulri(z,(GEN)zu[1]);
  /* z = Res (zeta_K, s = 1) * w D^(1/2) / [ 2^r1 (2pi)^r2 ] = h R */
  p1 = gmul2n(divir(clh,z), 1);
  /* c1 should be close to 2, and not much smaller */
  c1 = compute_check(sublambda,p1,&parch,&reg);
  if (!c1 || gcmpgs(gmul2n(c1,1),3) < 0)
  { /* precision problems */
    prec=(PRECREG<<1)-2;
    if (DEBUGLEVEL) err(warnprec,"buchall (compute_check)",prec);
    avma = av0; nf = nfnewprec(nf,prec);
    av0 = avma; cbach /= 2;
    goto INCREASEGEN;
  }
  if (gcmpgs(c1,3) > 0)
  {
    if (++nrelsup <= MAXRELSUP)
    {
      if (DEBUGLEVEL)
      {
        fprintferr("\n ***** check = %f\n",gtodouble(c1)/2);
        flusherr();
      }
      nlze=MIN_EXTRA; goto LABELINT;
    }
    if (cbach<11.99) { sfb_increase=1; goto LABELINT; }
    err(warner,"suspicious check. Try to increase extra relations");
  }

  /* Phase "be honest" */
  if (KCZ2 > KCZ)
  {
    if (!powsubfb)
      powsubfbgen(nf,subfb,CBUCHG+1,PRECREG,PRECREGINT);
    if (!be_honest(nf,subfb,RU,PRECREGINT)) goto INCREASEGEN;
  }

  /* regulator, roots of unity, fundamental units */
  if (flun < 0 || flun > 1)
  {
    xarch = cleancol(gmul(xarch,parch),N,PRECREG);
    if (DEBUGLEVEL) msgtimer("cleancol");
  }
  if (labs(flun) > 1)
  {
    fu = getfu(nf,&xarch,reg,flun,&k,PRECREG);
    if (k) fu = gmul((GEN)nf[7],fu);
    else if (labs(flun) > 2)
    {
      prec=(PRECREG<<1)-2;
      if (DEBUGLEVEL) err(warnprec,"buchall (getfu)",prec);
      avma = av0; nf = nfnewprec(nf,prec);
      av0 = avma; cbach /= 2;
      goto INCREASEGEN;
    }
  }

  /* class group generators */
  if (DEBUGLEVEL) fprintferr("\n");
  class_group_gen(nf,met,clh,u1,u2,vperm, &clg1, &clg2, PRECREGINT);

  /* cleanup */
  desallocate(matcopy);
  i = lg(C)-sreg; C += sreg; C[0] = evaltyp(t_MAT)|evallg(i);
  C = cleancol(C,N,PRECREG);
  settyp(vperm, t_COL);
  for (i=1; i<=KC; i++) vperm[i]=lstoi(vperm[i]);
  c1 = gdiv(gmul(reg,clh),z);

  return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c1,zu,W,B,xarch,C,vectbase,vperm));
}