===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/buch2.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/buch2.c	2001/10/02 11:17:03	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/buch2.c	2002/09/11 07:26:50	1.2
@@ -1,4 +1,4 @@
-/* $Id: buch2.c,v 1.1 2001/10/02 11:17:03 noro Exp $
+/* $Id: buch2.c,v 1.2 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -21,23 +21,28 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /*******************************************************************/
 #include "pari.h"
 #include "parinf.h"
+extern GEN gscal(GEN x,GEN y);
+extern GEN nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec);
+extern GEN get_nfindex(GEN bas);
+extern GEN sqred1_from_QR(GEN x, long prec);
+extern GEN computeGtwist(GEN nf, GEN vdir);
+extern GEN famat_mul(GEN f, GEN g);
+extern GEN famat_to_nf(GEN nf, GEN f);
+extern GEN famat_to_arch(GEN nf, GEN fa, long prec);
 extern GEN to_famat_all(GEN x, GEN y);
 extern int addcolumntomatrix(GEN V, GEN invp, GEN L);
 extern double check_bach(double cbach, double B);
 extern GEN gmul_mat_smallvec(GEN x, GEN y);
 extern GEN gmul_mati_smallvec(GEN x, GEN y);
+extern GEN get_arch(GEN nf,GEN x,long prec);
 extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
-extern GEN get_roots(GEN x,long r1,long ru,long prec);
-extern void get_nf_matrices(GEN nf, long small);
-extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *t, long v);
-extern GEN init_idele(GEN nf);
+extern GEN get_roots(GEN x,long r1,long prec);
+extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *t);
 extern GEN norm_by_embed(long r1, GEN x);
 extern void minim_alloc(long n,double ***q,long **x,double **y,double **z,double **v);
-extern GEN idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n);
 extern GEN arch_mul(GEN x, GEN y);
-extern GEN vecdiv(GEN x, GEN y);
-extern GEN vecmul(GEN x, GEN y);
-extern GEN mul_real(GEN x, GEN y);
+extern void wr_rel(GEN col);
+extern void dbg_rel(long s, GEN col);
 
 #define SFB_MAX 2
 #define SFB_STEP 2
@@ -48,132 +53,161 @@ 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 GEN *idealbase, vectbase, FB, numFB, powsubFB, numideal;
+/* used by factor[elt|gen|gensimple] to return factorizations of smooth elts
+ * HACK: MAX_FACTOR_LEN never checked, we assume default value is enough
+ * (since elts have small norm) */
+static long primfact[500], exprimfact[500];
 
-/* FB[i]     i-th rational prime used in factor base
- * numFB[i]  index k such that FB[k]=i (0 if i is not prime)
+/* a factor base contains only non-inert primes
+ * KC = # of P in factor base (p <= n, NP <= n2)
+ * KC2= # of P assumed to generate class group (NP <= n2)
  *
- * vectbase          vector of all ideals in FB
- * vecbase o subFB   part of FB used to build random relations
- * powsubFB  array lg(subFB) x (CBUCHG+1)   all powers up to CBUCHG
- *
- * idealbase[i]      prime ideals above i in FB
- * numideal[i]       index k such that idealbase[k] = i.
- *
- * mat            all relations found (as long integers, not reduced)
- * cglob          lg(mat) = total number of relations found so far
- *
- * 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 ideals counted by KC
- *   KCZ2= same for KC2
- */
+ * KCZ = # of rational primes under ideals counted by KC
+ * KCZ2= same for KC2 */
+typedef struct {
+  GEN FB; /* FB[i] = i-th rational prime used in factor base */
+  GEN LP; /* vector of all prime ideals in FB */
+  GEN *LV; /* LV[p] = vector of P|p, NP <= n2
+            * isclone() is set for LV[p] iff all P|p are in FB
+            * LV[i], i not prime or i > n2, is undefined! */
+  GEN iLP; /* iLP[p] = i such that LV[p] = [LP[i],...] */
+  long KC, KCZ, KCZ2;
+  GEN subFB; /* LP o subFB =  part of FB used to build random relations */
+  GEN powsubFB; /* array of [P^0,...,P^CBUCHG], P in LP o subFB */
+  GEN perm; /* permutation of LP used to represent relations [updated by
+               hnfspec/hnfadd: dense rows come first] */
+} FB_t;
 
 /* x a t_VECSMALL */
 static long
 ccontent(GEN x)
 {
-  long i, l = lg(x), s=labs(x[1]);
-  for (i=2; i<l && s>1; i++) s = cgcd(x[i],s);
+  long i, l = lg(x), s = labs(x[1]);
+  for (i=2; i<l && s!=1; i++) s = cgcd(x[i],s);
   return s;
 }
 
 static void
-desallocate(GEN **M)
+desallocate(GEN **M, GEN *first_nz)
 {
   GEN *m = *M;
   long i;
   if (m)
   {
     for (i=lg(m)-1; i; i--) free(m[i]);
-    free(m); *M = NULL;
+    free((void*)*M); *M = NULL;
+    free((void*)*first_nz); *first_nz = NULL;
   }
 }
 
-/* Return the list of indexes or the primes chosen for the subFB.
- * Fill vperm (if !=0): primes ideals sorted by increasing norm (except the
- * ones in subFB come first [dense rows come first for hnfspec])
- * ss = number of rational primes whose divisors are all in FB
- */
-static GEN
-subFBgen(long N,long m,long minsFB,GEN vperm, long *ptss)
+GEN
+cgetalloc(GEN x, size_t l, long t)
 {
-  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;
+  x = (GEN)gprealloc((void*)x, l * sizeof(long));
+  x[0] = evaltyp(t) | evallg(l); return x;
+}
 
-  (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]);
+static void
+reallocate(long max, GEN *M, GEN *first_nz)
+{
+  size_t l = max+1;
+  *M = cgetalloc(*M, l, t_VEC);
+  *first_nz = cgetalloc(*first_nz, l, t_VECSMALL);
+}
 
-    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; }
+/* don't take P|p if P ramified, or all other Q|p already there */
+static int
+ok_subFB(FB_t *F, long t, GEN D)
+{
+  GEN LP, P = (GEN)F->LP[t];
+  long p = itos((GEN)P[1]);
+  LP = F->LV[p];
+  return smodis(D,p) && (!isclone(LP) || t != F->iLP[p] + lg(LP)-1);
+}
 
-    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 = s1 = 0;
-    }
+/* set subFB, reset F->powsubFB
+ * Fill F->perm (if != NULL): primes ideals sorted by increasing norm (except
+ * the ones in subFB come first [dense rows for hnfspec]) */
+static int
+subFBgen(FB_t *F, GEN nf, double PROD, long minsFB)
+{
+  GEN y, perm, yes, no, D = (GEN)nf[3];
+  long i, j, k, iyes, ino, lv = F->KC + 1;
+  double prod;
+  const int init = (F->perm == NULL);
+  gpmem_t av;
+
+  if (init)
+  {
+    F->LP   = cgetg(lv, t_VEC);
+    F->perm = cgetg(lv, t_VECSMALL);
   }
-  if (z+minsFB >= lv) return NULL;
 
-  prod = 1.0;
-  perm = sindexsort(y1) + z; /* skip "zeroes" (excluded ideals) */
-  for(;;)
+  av = avma;
+  y = cgetg(lv,t_COL); /* Norm P */
+  for (k=0, i=1; i <= F->KCZ; i++)
   {
-    if (++n > minsFB && (z+n >= lv || prod > m + 0.5)) break;
-    prod *= gtodouble((GEN)y1[perm[n]]);
+    GEN LP = F->LV[F->FB[i]];
+    long l = lg(LP);
+    for (j = 1; j < l; j++)
+    {
+      GEN P = (GEN)LP[j];
+      k++;
+      y[k]     = (long)powgi((GEN)P[1], (GEN)P[4]);
+      F->LP[k] = (long)P; /* noop if init = 1 */
+    }
   }
-  if (prod < m) return NULL;
-  n--;
+  /* perm sorts LP by increasing norm */
+  perm = sindexsort(y);
+  no  = cgetg(lv, t_VECSMALL); ino  = 1;
+  yes = cgetg(lv, t_VECSMALL); iyes = 1;
+  prod = 1.0;
+  for (i = 1; i < lv; i++)
+  {
+    long t = perm[i];
+    if (!ok_subFB(F, t, D)) { no[ino++] = t; continue; }
 
-  /* 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)
+    yes[iyes++] = t;
+    prod *= (double)itos((GEN)y[t]);
+    if (iyes > minsFB && prod > PROD) break;
+  }
+  if (i == lv) return 0;
+  avma = av; /* HACK: gcopy(yes) still safe */
+  if (init)
   {
-    for (j=1; j<=n; j++) vperm[j] = perm[j];
-    for (   ; j<lv; j++) vperm[j] = perm1[j];
+    GEN p = F->perm;
+    for (j=1; j<iyes; j++)     p[j] = yes[j];
+    for (i=1; i<ino; i++, j++) p[j] =  no[i];
+    for (   ; j<lv; j++)       p[j] =  perm[j];
   }
-  for (j=1; j<=n; j++) subFB[j] = perm[j];
+  setlg(yes, iyes);
+  F->subFB = gcopy(yes);
+  F->powsubFB = NULL; return 1;
+}
+static int
+subFBgen_increase(FB_t *F, GEN nf, long step)
+{
+  GEN yes, D = (GEN)nf[3];
+  long i, iyes, lv = F->KC + 1, minsFB = lg(F->subFB)-1 + step;
 
-  if (DEBUGLEVEL)
+  yes = cgetg(minsFB+1, t_VECSMALL); iyes = 1;
+  for (i = 1; i < lv; i++)
   {
-    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");
-      outerr(vecextract_p(vectbase,subFB));
-      fprintferr("\n***** INITIAL PERMUTATION *****\n\n");
-      fprintferr("vperm = %Z\n\n",vperm);
-    }
-    msgtimer("sub factorbase (%ld elements)",n);
+    long t = F->perm[i];
+    if (!ok_subFB(F, t, D)) continue;
+
+    yes[iyes++] = t;
+    if (iyes > minsFB) break;
   }
-  powsubFB = NULL;
-  *ptss = ss; return subFB;
+  if (i == lv) return 0;
+  F->subFB = yes;
+  F->powsubFB = NULL; return 1;
 }
 
 static GEN
 mulred(GEN nf,GEN x, GEN I, long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN y = cgetg(3,t_VEC);
 
   y[1] = (long)idealmulh(nf,I,(GEN)x[1]);
@@ -185,22 +219,21 @@ mulred(GEN nf,GEN x, GEN I, long prec)
 
 /* Compute powers of prime ideals (P^0,...,P^a) in subFB (a > 1)
  * powsubFB[j][i] contains P_i^j in LLL form + archimedean part in
- * MULTIPLICATIVE form (logs are expensive)
- */
+ * MULTIPLICATIVE form (logs are expensive) */
 static void
-powsubFBgen(GEN nf,GEN subFB,long a,long prec)
+powsubFBgen(FB_t *F, GEN nf, long a, long prec)
 {
-  long i,j, n = lg(subFB), RU = lg(nf[6]);
+  long i,j, n = lg(F->subFB), RU = lg(nf[6]);
   GEN *pow, arch0 = cgetg(RU,t_COL);
   for (i=1; i<RU; i++) arch0[i] = un; /* 0 in multiplicative notations */
 
   if (DEBUGLEVEL) fprintferr("Computing powers for sub-factor base:\n");
-  powsubFB = cgetg(n, t_VEC);
+  F->powsubFB = cgetg(n, t_VEC);
   for (i=1; i<n; i++)
   {
-    GEN vp = (GEN)vectbase[subFB[i]];
+    GEN vp = (GEN)F->LP[ F->subFB[i] ];
     GEN z = cgetg(3,t_VEC); z[1]=vp[1]; z[2]=vp[2];
-    pow = (GEN*)cgetg(a+1,t_VEC); powsubFB[i] = (long)pow;
+    pow = (GEN*)cgetg(a+1,t_VEC); F->powsubFB[i] = (long)pow;
     pow[1]=cgetg(3,t_VEC);
     pow[1][1] = (long)z;
     pow[1][2] = (long)arch0;
@@ -215,273 +248,276 @@ powsubFBgen(GEN nf,GEN subFB,long a,long prec)
   if (DEBUGLEVEL) msgtimer("powsubFBgen");
 }
 
-/* Compute FB, numFB, idealbase, vectbase, numideal.
+/* Compute FB, LV, iLP + KC*. Reset perm
  * n2: bound for norm of tested prime ideals (includes be_honest())
- * n : bound for prime ideals used to build relations (Bach cst) ( <= n2 )
+ * n : bound for p, such that P|p (NP <= n2) used to build relations
 
  * 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)
- */
+ * close to residue of zeta_K at 1 = 2^r1 (2pi)^r2 h R / (w D) */
 static GEN
-FBgen(GEN nf,long n2,long n)
+FBgen(FB_t *F, GEN nf,long n2,long n)
 {
-  byteptr delta=diffptr;
-  long KC2,i,j,k,p,lon,ip,nor, N = degpol(nf[1]);
-  GEN p2,p1,NormP,lfun;
-  long prim[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0 };
+  byteptr delta = diffptr;
+  long i, p, ip;
+  GEN prim, Res;
 
-  numFB    = cgetg(n2+1,t_VECSMALL);
-  FB       = cgetg(n2+1,t_VECSMALL);
-  numideal = cgetg(n2+1,t_VECSMALL);
-  idealbase= (GEN*)cgetg(n2+1,t_VEC);
+  if (maxprime() < n2) err(primer1);
+  F->FB  = cgetg(n2+1, t_VECSMALL);
+  F->iLP = cgetg(n2+1, t_VECSMALL);
+  F->LV = (GEN*)new_chunk(n2+1);
 
-  lfun=realun(DEFAULTPREC);
-  p=*delta++; i=0; ip=0; KC=0;
-  while (p<=n2)
+  Res = realun(DEFAULTPREC);
+  prim = icopy(gun);
+  i = ip = 0;
+  F->KC = F->KCZ = 0;
+  for (p = 0;;) /* 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 */
+    gpmem_t av = avma, av1;
+    long k, l;
+    GEN P, a, b;
+
+    NEXT_PRIME_VIADIFF(p, delta);
+    if (!F->KC && p > n) { F->KCZ = i; F->KC = ip; }
+    if (p > n2) break;
+
+    if (DEBUGLEVEL>1) { fprintferr(" %ld",p); flusherr(); }
+    prim[2] = p; P = primedec(nf,prim); l = lg(P);
+
+    /* a/b := (p-1)/p * prod_{P|p, NP <= n2} NP / (NP-1) */
+    av1 = avma; a = b = NULL;
+    for (k=1; k<l; k++)
     {
-      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++; numFB[p]=i; FB[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;
+      GEN NormP = powgi(prim, gmael(P,k,4));
+      long nor;
+      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 FB */
-      else
-        { setlg(p1,k); p1 = gerepilecopy(av,p1); }
-      idealbase[i] = p1;
+      if (a) { a = mului(nor, a); b = mului(nor-1, b); }
+      else   { a = utoi(nor / p); b = utoi((nor-1) / (p-1)); }
     }
-    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);
+    if (a) affrr(divri(mulir(a,Res),b),   Res);
+    else   affrr(divrs(mulsr(p-1,Res),p), Res);
+    avma = av1;
+    if (l == 2 && itos(gmael(P,1,3)) == 1) continue; /* p inert */
 
-  setlg(FB,KCZ2);
-  setlg(numFB,KCZ2);
-  setlg(numideal,KCZ2);
-  setlg(idealbase,KCZ2);
-  vectbase=cgetg(KC+1,t_COL);
-  for (i=1; i<=KCZ; i++)
-  {
-    p1 = idealbase[i]; k=lg(p1);
-    p2 = vectbase + numideal[FB[i]];
-    for (j=1; j<k; j++) p2[j]=p1[j];
+    /* keep non-inert ideals with Norm <= n2 */
+    if (k == l)
+      setisclone(P); /* flag it: all prime divisors in FB */
+    else
+      { setlg(P,k); P = gerepilecopy(av,P); }
+    F->FB[++i]= p;
+    F->LV[p]  = P;
+    F->iLP[p] = ip; ip += k-1;
   }
+  if (! F->KC) return NULL;
+  setlg(F->FB, F->KCZ+1); F->KCZ2 = i;
   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]);
+      fprintferr("KC2=%ld, KC=%ld, KCZ=%ld, KCZ2=%ld\n",
+                  ip, F->KC, F->KCZ, F->KCZ2);
+      for (i=1; i<=F->KCZ; i++) fprintferr("++ LV[%ld] = %Z",i,F->LV[F->FB[i]]);
     }
     msgtimer("factor base");
   }
-  return lfun;
+  F->perm = NULL; return Res;
 }
 
-/* can we factor I / m ? (m pseudo minimum, computed in ideallllredpart1) */
-static long
-factorgen(GEN nf,GEN idealvec,long kcz,long limp)
+/*  SMOOTH IDEALS */
+static void
+store(long i, long e)
 {
-  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= absi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */
+  primfact[++primfact[0]] = i; /* index */
+  exprimfact[primfact[0]] = e; /* exponent */
+}
 
-  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++)
+/* divide out x by all P|p, where x as in can_factor().  k = v_p(Nx) */
+static int
+divide_p_elt(FB_t *F, long p, long k, GEN nf, GEN m)
+{
+  GEN P, LP = F->LV[p];
+  long j, v, l = lg(LP), ip = F->iLP[p];
+  for (j=1; j<l; j++)
   {
-    p=FB[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;
+    P = (GEN)LP[j];
+    v = int_elt_val(nf, m, (GEN)P[1], (GEN)P[5], NULL); /* v_P(m) */
+    if (!v) continue;
+    store(ip + j, v); /* v = v_P(m) > 0 */
+    k -= v * itos((GEN)P[4]);
+    if (!k) return 1;
   }
-  if (cmpis(x,limp) > 0) return 0;
-
-  ifinal = i; lo = 0;
-  for (i=1; i<=ifinal; i++)
+  return 0;
+}
+static int
+divide_p_id(FB_t *F, long p, long k, GEN nf, GEN I)
+{
+  GEN P, LP = F->LV[p];
+  long j, v, l = lg(LP), ip = F->iLP[p];
+  for (j=1; j<l; j++)
   {
-    k = listexpo[i];
-    if (k)
-    {
-      p = FB[i]; p1 = idealbase[numFB[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;
-    }
+    P = (GEN)LP[j];
+    v = idealval(nf,I, P);
+    if (!v) continue;
+    store(ip + j, v); /* v = v_P(I) > 0 */
+    k -= v * itos((GEN)P[4]);
+    if (!k) return 1;
   }
-  if (is_pm1(x)) { primfact[0]=lo; return 1; }
-
-  p = itos(x); p1 = idealbase[numFB[p]];
-  n1 = lg(p1); ip = numideal[p];
-  for (k=1,j=1; j<n1; j++)
+  return 0;
+}
+static int
+divide_p_quo(FB_t *F, long p, long k, GEN nf, GEN I, GEN m)
+{
+  GEN P, LP = F->LV[p];
+  long j, v, l = lg(LP), ip = F->iLP[p];
+  for (j=1; j<l; 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; }
-    }
+    P = (GEN)LP[j];
+    v = int_elt_val(nf, m, (GEN)P[1], (GEN)P[5], NULL); /* v_P(m) */
+    if (!v) continue;
+    v = idealval(nf,I, P) - v;
+    if (!v) continue;
+    store(ip + j, v); /* v = v_P(I / m) < 0 */
+    k += v * itos((GEN)P[4]);
+    if (!k) return 1;
   }
   return 0;
 }
 
-/* can we factor x ? Nx = norm(x) */
-static long
-factorelt(GEN nf,GEN cbase,GEN x,GEN Nx,long kcz,long limp)
+/* is *N > 0 a smooth rational integer wrt F ? (put the exponents in *ex) */
+static int
+smooth_int(FB_t *F, GEN *N, GEN *ex)
 {
-  long i,j,n1,ip,v,p,k,lo,ifinal;
-  GEN q,r,P,p1,listexpo;
+  GEN q, r, FB = F->FB;
+  const long KCZ = F->KCZ;
+  const long limp = FB[KCZ]; /* last p in FB */
+  long i, p, k;
 
-  if (is_pm1(Nx)) { primfact[0]=0; return 1; }
-  listexpo = new_chunk(kcz+1);
+  *ex = new_chunk(KCZ+1);
   for (i=1; ; i++)
   {
-    p=FB[i]; q=dvmdis(Nx,p,&r);
-    for (k=0; !signe(r); k++) { Nx=q; q=dvmdis(Nx,p,&r); }
-    listexpo[i] = k;
-    if (cmpis(q,p)<=0) break;
-    if (i==kcz) return 0;
+    p = FB[i]; q = dvmdis(*N,p,&r);
+    for (k=0; !signe(r); k++) { *N = q; q = dvmdis(*N, p, &r); }
+    (*ex)[i] = k;
+    if (cmpis(q,p) <= 0) break;
+    if (i == KCZ) return 0;
   }
-  if (cmpis(Nx,limp) > 0) return 0;
+  (*ex)[0] = i;
+  return (cmpis(*N,limp) <= 0);
+}
 
-  if (cbase) x = gmul(cbase,x);
-  ifinal=i; lo = 0;
-  for (i=1; i<=ifinal; i++)
-  {
-    k = listexpo[i];
-    if (k)
-    {
-      p = FB[i]; p1 = idealbase[numFB[p]];
-      n1 = lg(p1); ip = numideal[p];
-      for (j=1; j<n1; j++)
-      {
-        P = (GEN)p1[j];
-	v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, 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(Nx)) { primfact[0]=lo; return 1; }
+static int
+divide_p(FB_t *F, long p, long k, GEN nf, GEN I, GEN m)
+{
+  if (!m) return divide_p_id (F,p,k,nf,I);
+  if (!I) return divide_p_elt(F,p,k,nf,m);
+  return divide_p_quo(F,p,k,nf,I,m);
+}
 
-  p = itos(Nx); p1 = idealbase[numFB[p]];
-  n1 = lg(p1); ip = numideal[p];
-  for (k=1,j=1; j<n1; j++)
+
+/* Let x = m if I == NULL,
+ *         I if m == NULL,
+ *         I/m otherwise.
+ * Can we factor x ? N = Norm x > 0 */
+static long
+can_factor(FB_t *F, GEN nf, GEN I, GEN m, GEN N)
+{
+  GEN ex;
+  long i;
+  primfact[0] = 0;
+  if (is_pm1(N)) return 1;
+  if (!smooth_int(F, &N, &ex)) return 0;
+  for (i=1; i<=ex[0]; i++)
+    if (ex[i] && !divide_p(F, F->FB[i], ex[i], nf, I, m)) return 0;
+  return is_pm1(N) || divide_p(F, itos(N), 1, nf, I, m);
+}
+
+/* can we factor I/m ? [m in I from pseudomin] */
+static long
+factorgen(FB_t *F, GEN nf, GEN I, GEN m)
+{
+  GEN Nm = absi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */
+  GEN N  = diviiexact(Nm, dethnf_i(I)); /* N(m / I) */
+  return can_factor(F, nf, I, m, N);
+}
+
+/*  FUNDAMENTAL UNITS */
+
+/* a, m real. Return  (Re(x) + a) + I * (Im(x) % m) */
+static GEN
+addRe_modIm(GEN x, GEN a, GEN m)
+{
+  GEN re, im, z;
+  if (typ(x) == t_COMPLEX)
   {
-    P = (GEN)p1[j];
-    v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, k);
-    if (v)
+    re = gadd((GEN)x[1], a);
+    im = gmod((GEN)x[2], m);
+    if (gcmp0(im)) z = re;
+    else
     {
-      primfact[++lo]=ip+j; expoprimfact[lo]=v;
-      k -= v*itos((GEN)P[4]);
-      if (!k) { primfact[0]=lo; return 1; }
+      z = cgetg(3,t_COMPLEX);
+      z[1] = (long)re;
+      z[2] = (long)im;
     }
   }
-  return 0;
+  else
+    z = gadd(x, a);
+  return z;
 }
 
+/* clean archimedean components */
 static GEN
-cleancol(GEN x,long N,long PRECREG)
+cleanarch(GEN x, long N, long prec)
 {
-  long i,j,av,tetpil,tx=typ(x),R1,RU;
-  GEN s,s2,re,pi4,im,y;
+  long i, R1, RU, tx = typ(x);
+  GEN s, y, pi2;
 
-  if (tx==t_MAT)
+  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);
+    y = cgetg(lg(x), tx);
+    for (i=1; i < lg(x); i++)
+      y[i] = (long)cleanarch((GEN)x[i], N, prec);
     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]);
-  im = gimag(x); s = gdivgs(s,-N);
-  s2 = (N>R1)? gmul2n(s,1): NULL;
-  pi4=gmul2n(mppi(PRECREG),2);
-  tetpil=avma; y=cgetg(RU+1,tx);
-  for (i=1; i<=RU; i++)
+  if (!is_vec_t(tx)) err(talker,"not a vector/matrix in cleanarch");
+  RU = lg(x)-1; R1 = (RU<<1)-N;
+  s = gdivgs(sum(greal(x), 1, RU), -N); /* -log |norm(x)| / N */
+  y = cgetg(RU+1,tx);
+  pi2 = Pi2n(1, prec);
+  for (i=1; i<=R1; i++) y[i] = (long)addRe_modIm((GEN)x[i], s, pi2);
+  if (i <= RU)
   {
-    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], pi4);
+    GEN pi4 = Pi2n(2, prec), s2 = gmul2n(s, 1);
+    for (   ; i<=RU; i++) y[i] = (long)addRe_modIm((GEN)x[i], s2, pi4);
   }
-  return gerepile(av,tetpil,y);
+  return y;
 }
 
-#define RELAT 0
-#define LARGE 1
-#define PRECI 2
+enum { RELAT, LARGE, PRECI };
+
 static GEN
-not_given(long av, long flun, long reason)
+not_given(gpmem_t av, long fl, long reason)
 {
-  if (labs(flun)==2)
+  if (! (fl & nf_FORCE))
   {
     char *s;
     switch(reason)
     {
-      case RELAT: s = "not enough relations for fundamental units"; break;
       case LARGE: s = "fundamental units too large"; break;
       case PRECI: s = "insufficient precision for fundamental units"; break;
       default: s = "unknown problem with fundamental units";
     }
     err(warner,"%s, not given",s);
   }
-  avma=av; return cgetg(1,t_MAT);
+  avma = av; return cgetg(1,t_MAT);
 }
 
 /* check whether exp(x) will get too big */
 static long
 expgexpo(GEN x)
 {
-  long i,j,e, E = -HIGHEXPOBIT;
+  long i,j,e, E = - (long)HIGHEXPOBIT;
   GEN p1;
 
   for (i=1; i<lg(x); i++)
@@ -527,89 +563,83 @@ gauss_realimag(GEN x, GEN y)
   y = split_realimag(y,r1,r2); return gauss(M, y);
 }
 
-static GEN
-getfu(GEN nf,GEN *ptxarch,GEN reg,long flun,long *pte,long prec)
+GEN
+getfu(GEN nf,GEN *ptA,GEN reg,long fl,long *pte,long prec)
 {
-  long av = avma,e,i,j,R1,RU,N=degpol(nf[1]);
-  GEN p1,p2,u,y,matep,s,xarch,vec;
-  GEN *gptr[2];
+  long e, i, j, R1, RU, N=degpol(nf[1]);
+  gpmem_t av = avma;
+  GEN p1,p2,u,y,matep,s,A,vec;
 
   if (DEBUGLEVEL) fprintferr("\n#### Computing fundamental units\n");
   R1 = itos(gmael(nf,2,1)); RU = (N+R1)>>1;
   if (RU==1) { *pte=BIGINT; return cgetg(1,t_VEC); }
 
-  *pte = 0; xarch = *ptxarch;
-  if (gexpo(reg) < -8) return not_given(av,flun,RELAT);
-
+  *pte = 0; A = *ptA;
   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 = gzero; for (i=1; i<=RU; i++) s = gadd(s,greal(gcoeff(A,i,j)));
     s = gdivgs(s, -N);
     p1=cgetg(RU+1,t_COL); matep[j]=(long)p1;
-    for (i=1; i<=R1; i++) p1[i] = ladd(s, gcoeff(xarch,i,j));
-    for (   ; i<=RU; i++) p1[i] = ladd(s, gmul2n(gcoeff(xarch,i,j),-1));
+    for (i=1; i<=R1; i++) p1[i] = ladd(s, gcoeff(A,i,j));
+    for (   ; i<=RU; i++) p1[i] = ladd(s, gmul2n(gcoeff(A,i,j),-1));
   }
-  if (prec <= 0) prec = gprecision(xarch);
-  u = lllintern(greal(matep),1,prec);
-  if (!u) return not_given(av,flun,PRECI);
+  if (prec <= 0) prec = gprecision(A);
+  u = lllintern(greal(matep),100,1,prec);
+  if (!u) return not_given(av,fl,PRECI);
 
   p1 = gmul(matep,u);
-  if (expgexpo(p1) > 20) { *pte = BIGINT; return not_given(av,flun,LARGE); }
+  if (expgexpo(p1) > 20) { *pte = BIGINT; return not_given(av,fl,LARGE); }
   matep = gexp(p1,prec);
   y = grndtoi(gauss_realimag(nf,matep), &e);
   *pte = -e;
-  if (e>=0) return not_given(av,flun,PRECI);
+  if (e >= 0) return not_given(av,fl,PRECI);
   for (j=1; j<RU; j++)
     if (!gcmp1(idealnorm(nf, (GEN)y[j]))) break;
-  if (j < RU) { *pte = 0; return not_given(av,flun,PRECI); }
-  xarch = gmul(xarch,u);
+  if (j < RU) { *pte = 0; return not_given(av,fl,PRECI); }
+  A = gmul(A,u);
 
   /* y[i] are unit generators. Normalize: smallest L2 norm + lead coeff > 0 */
   y = gmul((GEN)nf[7], y);
-  vec = cgetg(RU+1,t_COL); p2 = mppi(prec);
-  p1 = pureimag(p2);
-  p2 = pureimag(gmul2n(p2,1));
-  for (i=1; i<=R1; i++) vec[i]=(long)p1;
-  for (   ; i<=RU; i++) vec[i]=(long)p2;
+  vec = cgetg(RU+1,t_COL);
+  p1 = PiI2n(0,prec); for (i=1; i<=R1; i++) vec[i] = (long)p1;
+  p2 = PiI2n(1,prec); for (   ; i<=RU; i++) vec[i] = (long)p2;
   for (j=1; j<RU; j++)
   {
-    p1 = (GEN)y[j]; p2 = ginvmod(p1, (GEN)nf[1]);
+    p1 = (GEN)y[j]; p2 = QX_invmod(p1, (GEN)nf[1]);
     if (gcmp(QuickNormL2(p2,DEFAULTPREC),
              QuickNormL2(p1,DEFAULTPREC)) < 0)
     {
-      xarch[j] = lneg((GEN)xarch[j]);
+      A[j] = lneg((GEN)A[j]);
       p1 = p2;
     }
     if (gsigne(leading_term(p1)) < 0)
     {
-      xarch[j] = ladd((GEN)xarch[j], vec);
+      A[j] = ladd((GEN)A[j], vec);
       p1 = gneg(p1);
     }
     y[j] = (long)p1;
   }
   if (DEBUGLEVEL) msgtimer("getfu");
-  *ptxarch=xarch; gptr[0]=ptxarch; gptr[1]=&y;
-  gerepilemany(av,gptr,2); return y;
+  gerepileall(av,2, &A, &y);
+  *ptA = A; return y;
 }
-#undef RELAT
-#undef LARGE
-#undef PRECI
 
 GEN
 buchfu(GEN bnf)
 {
-  long av = avma, c;
-  GEN nf,xarch,reg,res, y = cgetg(3,t_VEC);
+  long c;
+  gpmem_t av = avma;
+  GEN nf,A,reg,res, y = cgetg(3,t_VEC);
 
-  bnf = checkbnf(bnf); xarch = (GEN)bnf[3]; nf = (GEN)bnf[7];
+  bnf = checkbnf(bnf); A = (GEN)bnf[3]; nf = (GEN)bnf[7];
   res = (GEN)bnf[8]; reg = (GEN)res[2];
   if (lg(res)==7 && lg(res[5])==lg(nf[6])-1)
   {
     y[1] = lcopy((GEN)res[5]);
     y[2] = lcopy((GEN)res[6]); return y;
   }
-  y[1] = (long)getfu(nf,&xarch,reg,2,&c,0);
+  y[1] = (long)getfu(nf, &A, reg, nf_UNITS, &c, 0);
   y[2] = lstoi(c); return gerepilecopy(av, y);
 }
 
@@ -665,197 +695,280 @@ get_norm_fact(GEN gen, GEN ex, GEN *pd)
   *pd = d; return N;
 }
 
-/* try to split ideal / (some integer) on FB */
-static long
-factorgensimple(GEN nf,GEN ideal)
+static GEN
+get_pr_lists(GEN FB, long N, int list_pr)
 {
-  long N,i,v,av1 = avma,lo, L = lg(vectbase);
-  GEN x;
-  if (typ(ideal) != t_MAT) ideal = (GEN)ideal[1]; /* idele */
-  x = dethnf_i(ideal);
-  N = lg(ideal)-1;
-  if (gcmp1(x)) { avma=av1; primfact[0]=0; return 1; }
-  for (lo=0, i=1; i<L; i++)
-  {
-    GEN q,y, P=(GEN)vectbase[i], p=(GEN)P[1];
-    long vx0, vx, sum_ef, e,f,lo0,i0;
+  GEN pr, L;
+  long i, l = lg(FB), p, pmax;
 
-    vx = pvaluation(x,p,&y);
-    if (!vx) continue;
-    /* p | x = Nideal */
-    vx0 = vx; sum_ef = 0; lo0 =lo; i0 = i;
-    for(;;)
+  pmax = 0;
+  for (i=1; i<l; i++)
+  {
+    pr = (GEN)FB[i]; p = itos((GEN)pr[1]);
+    if (p > pmax) pmax = p;
+  }
+  L = cgetg(pmax+1, t_VEC);
+  for (p=1; p<=pmax; p++) L[p] = 0;
+  if (list_pr)
+  {
+    for (i=1; i<l; i++)
     {
-      e = itos((GEN)P[3]);
-      f = itos((GEN)P[4]); sum_ef += e * f;
-      v = idealval(nf,ideal,P);
-      if (v)
-      {
-        lo++; primfact[lo]=i; expoprimfact[lo]=v;
-        vx -= v * f; if (!vx) break;
-      }
-      i++; if (i == L) break;
-      P = (GEN)vectbase[i]; q = (GEN)P[1];
-      if (!egalii(p,q)) break;
+      pr = (GEN)FB[i]; p = itos((GEN)pr[1]);
+      if (!L[p]) L[p] = (long)cget1(N+1, t_VEC);
+      appendL((GEN)L[p], pr);
     }
-    if (vx)
-    { /* divisible by P | p not in FB */
-      long k,l,n;
-      k = N - sum_ef;
-      if (vx % k) break;
-      k = vx / k; /* ideal / p^k may factor on FB */
-      for (l = lo0+1; l <= lo; l++)
-      {
-        P = (GEN)vectbase[primfact[l]];
-        expoprimfact[l] -= k * itos((GEN)P[3]); /* may become 0 */
-      }
-      n = lo0+1;
-      for (l = i0; l < i; l++) /* update exponents for other P | p */
-      {
-        if (n <= lo && primfact[n] == l) { n++; continue; }
-        lo++; primfact[lo] = l; P = (GEN)vectbase[l];
-        expoprimfact[lo] = - k * itos((GEN)P[3]);
-      }
-      /* check that ideal / p^k / (tentative factorisation) is integral */
-      {
-        GEN z = ideal;
-        for (l = lo0+1; l <= lo; l++)
-        {
-          GEN n = stoi(- expoprimfact[l]);
-          z = idealmulpowprime(nf, z, (GEN)vectbase[primfact[l]], n);
-        }
-        z = gdiv(z, gpowgs(p, k));
-        if (!gcmp1(denom(z))) { avma=av1; return 0; }
-        ideal = z;
-      }
+    for (p=1; p<=pmax; p++)
+      if (L[p]) L[p] = (long)gen_sort((GEN)L[p],0,cmp_prime_over_p);
+  }
+  else
+  {
+    for (i=1; i<l; i++)
+    {
+      pr = (GEN)FB[i]; p = itos((GEN)pr[1]);
+      if (!L[p]) L[p] = (long)cget1(N+1, t_VECSMALL);
+      appendL((GEN)L[p], (GEN)i);
     }
-    if (gcmp1(y)) { avma=av1; primfact[0]=lo; return 1; }
-    x = y;
   }
-  avma=av1; return 0;
+  return L;
 }
 
-static void
-add_to_fact(long l, long v, long e)
+/* recover FB, LV, iLP, KCZ from Vbase */
+static GEN
+recover_partFB(FB_t *F, GEN Vbase, long N)
 {
-  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
+  GEN FB, LV, iLP, L = get_pr_lists(Vbase, N, 0);
+  long l = lg(L), p, ip, i;
+
+  i = ip = 0;
+  FB = cgetg(l, t_VECSMALL);
+  iLP= cgetg(l, t_VECSMALL);
+  LV = cgetg(l, t_VEC);
+#if 1 /* for debugging */
+  for (p=1;p<l;p++) FB[p]=iLP[p]=LV[p]=0;
+#endif
+  for (p = 2; p < l; p++)
   {
-    lo = ++primfact[0];
-    primfact[lo] = v;
-    expoprimfact[lo] = e;
+    if (!L[p]) continue;
+    FB[++i] = p;
+    LV[p] = (long)vecextract_p(Vbase, (GEN)L[p]);
+    iLP[p]= ip; ip += lg(L[p])-1;
   }
+  F->KCZ = i;
+  F->FB = FB; setlg(FB, i+1);
+  F->LV = (GEN*)LV;
+  F->iLP= iLP; return L;
 }
 
+static GEN
+init_famat(GEN x)
+{
+  GEN y = cgetg(3, t_VEC);
+  y[1] = (long)x;
+  y[2] = lgetg(1,t_MAT); return y;
+}
+
+/* add v^e to factorization */
 static void
-init_sub(long l, GEN perm, GEN *v, GEN *ex)
+add_to_fact(long v, long e)
 {
-  long i;
-  *v = cgetg(l,t_VECSMALL);
-  *ex= cgetg(l,t_VECSMALL);
-  for (i=1; i<l; i++) (*v)[i] = itos((GEN)perm[i]);
+  long i, l = primfact[0];
+  for (i=1; i<=l && primfact[i] < v; i++)/*empty*/;
+  if (i <= l && primfact[i] == v) exprimfact[i] += e; else store(v, e);
 }
 
-/* factor x = x0[1] on vectbase (modulo principal ideals).
- * We may have x0[2] = NULL (principal part not wanted), in which case,
- * don't compute archimedean parts */
+/* L (small) list of primes above the same p including pr. Return pr index */
+static int
+pr_index(GEN L, GEN pr)
+{
+  long j, l = lg(L);
+  GEN al = (GEN)pr[2];
+  for (j=1; j<l; j++)
+    if (gegal(al, gmael(L,j,2))) return j;
+  err(bugparier,"codeprime");
+  return 0; /* not reached */
+}
+
+static long
+Vbase_to_FB(FB_t *F, GEN pr)
+{
+  long p = itos((GEN)pr[1]);
+  return F->iLP[p] + pr_index(F->LV[p], pr);
+}
+
+/* return famat y (principal ideal) such that x / y is smooth [wrt Vbase] */
 static GEN
-split_ideal(GEN nf, GEN x0, long prec, GEN vperm)
+SPLIT(FB_t *F, GEN nf, GEN x, GEN Vbase)
 {
-  GEN p1,vdir,id,z,v,ex,y, x = (GEN)x0[1];
-  long nbtest_lim,nbtest,bou,i,ru, lgsub;
+  GEN vdir, id, z, ex, y, x0;
+  long nbtest_lim, nbtest, bou, i, ru, lgsub;
   int flag = (gexpo(gcoeff(x,1,1)) < 100);
 
-  y = x0;
-  if (flag && factorgensimple(nf,y)) return y;
+  /* try without reduction if x is small */
+  if (flag && can_factor(F, nf, x, NULL, dethnf_i(x))) return NULL;
 
-  y = ideallllred(nf,x0,NULL,prec);
-  if (flag &&  ((!x0[2] && gegal((GEN)y[1], (GEN)x[1]))
-             ||  (x0[2] && gcmp0((GEN)y[2])))) flag = 0; /* y == x0 */
-  if (flag && factorgensimple(nf,y)) return y;
+  /* if reduction fails (y scalar), do not retry can_factor */
+  y = idealred_elt(nf,x);
+  if ((!flag || !isnfscalar(y)) && factorgen(F, nf, x, y)) return y;
 
-  z = init_idele(nf); ru = lg(z[2]);
-  if (!x0[2]) { z[2] = 0; x0 = x; } /* stop cheating */
-  vdir = cgetg(ru,t_VEC);
-  for (i=2; i<ru; i++) vdir[i]=zero;
+  /* reduce in various directions */
+  ru = lg(nf[6]);
+  vdir = cgetg(ru,t_VECSMALL);
+  for (i=2; i<ru; i++) vdir[i]=0;
   for (i=1; i<ru; i++)
   {
-    vdir[i]=lstoi(10);
-    y = ideallllred(nf,x0,vdir,prec);
-    if (factorgensimple(nf,y)) return y;
-    vdir[i]=zero;
+    vdir[i] = 10;
+    y = ideallllred_elt(nf,x,vdir);
+    if (factorgen(F, nf, x, y)) return y;
+    vdir[i] = 0;
   }
-  nbtest = 0; nbtest_lim = (ru-1) << 2; lgsub = 3;
-  init_sub(lgsub, vperm, &v, &ex);
+
+  /* tough case, multiply by random products */
+  lgsub = 3;
+  ex = cgetg(lgsub, t_VECSMALL);
+  z  = init_famat(NULL);
+  x0 = init_famat(x);
+  nbtest = 1; nbtest_lim = 4;
   for(;;)
   {
-    int non0 = 0;
+    gpmem_t av = avma;
+    if (DEBUGLEVEL>2) fprintferr("# ideals tried = %ld\n",nbtest);
     id = x0;
     for (i=1; i<lgsub; i++)
     {
       ex[i] = mymyrand() >> randshift;
       if (ex[i])
-      { /* don't let id become too large as lgsub gets bigger: avoid
-           prec problems */
-        if (non0) id = ideallllred(nf,id,NULL,prec);
-        non0++;
-        z[1]=vectbase[v[i]]; p1=idealpowred(nf,z,stoi(ex[i]),prec);
-        id = idealmulh(nf,id,p1);
+      { /* avoid prec pb: don't let id become too large as lgsub increases */
+        if (id != x0) id = ideallllred(nf,id,NULL,0);
+        z[1] = Vbase[i];
+        id = idealmulh(nf, id, idealpowred(nf,z,stoi(ex[i]),0));
       }
     }
     if (id == x0) continue;
 
-    for (i=1; i<ru; i++) vdir[i] = lstoi(mymyrand() >> randshift);
+    for (i=1; i<ru; i++) vdir[i] = mymyrand() >> randshift;
     for (bou=1; bou<ru; bou++)
     {
-      if (bou>1)
+      y = ideallllred_elt(nf, (GEN)id[1], vdir);
+      if (factorgen(F, nf, (GEN)id[1], y))
       {
-        for (i=1; i<ru; i++) vdir[i]=zero;
-        vdir[bou]=lstoi(10);
+        for (i=1; i<lgsub; i++) 
+          if (ex[i]) add_to_fact(Vbase_to_FB(F,(GEN)Vbase[i]), -ex[i]);
+        return famat_mul((GEN)id[2], y);
       }
-      nbtest++;
-      y = ideallllred(nf,id,vdir,prec);
-      if (DEBUGLEVEL>2)
-        fprintferr("nbtest = %ld, ideal = %Z\n",nbtest,y[1]);
-      if (factorgensimple(nf,y))
-      {
-        long l = primfact[0];
-        for (i=1; i<lgsub; i++) add_to_fact(l,v[i],-ex[i]);
-        return y;
-      }
+      for (i=1; i<ru; i++) vdir[i] = 0;
+      vdir[bou] = 10;
     }
-    if (nbtest > nbtest_lim)
+    avma = av;
+    if (++nbtest > nbtest_lim)
     {
       nbtest = 0;
-      if (lgsub < 7)
+      if (++lgsub < 7)
       {
-        lgsub++; nbtest_lim <<= 2;
-        init_sub(lgsub, vperm, &v, &ex);
+        nbtest_lim <<= 1;
+        ex = cgetg(lgsub, t_VECSMALL);
       }
       else nbtest_lim = VERYBIGINT; /* don't increase further */
-      if (DEBUGLEVEL)
-        fprintferr("split_ideal: increasing factor base [%ld]\n",lgsub);
+      if (DEBUGLEVEL) fprintferr("SPLIT: increasing factor base [%ld]\n",lgsub);
     }
   }
 }
 
-static void
-get_split_expo(GEN xalpha, GEN yalpha, GEN vperm)
+static GEN
+split_ideal(GEN nf, GEN x, GEN Vbase)
 {
-  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;
+  FB_t F;
+  GEN L = recover_partFB(&F, Vbase, lg(x)-1);
+  GEN y = SPLIT(&F, nf, x, Vbase);
+  long p,j, i, l = lg(F.FB);
+
+  p = j = 0; /* -Wall */
   for (i=1; i<=primfact[0]; i++)
+  { /* decode index C = ip+j --> (p,j) */
+    long q,k,t, C = primfact[i];
+    for (t=1; t<l; t++)
+    {
+      q = F.FB[t];
+      k = C - F.iLP[q];
+      if (k <= 0) break;
+      p = q;
+      j = k;
+    }
+    primfact[i] = mael(L, p, j);
+  }
+  return y;
+}
+
+/* return sorted vectbase [sorted in bnf since version 2.2.4] */
+static GEN
+get_Vbase(GEN bnf)
+{
+  GEN vectbase = (GEN)bnf[5], perm = (GEN)bnf[6], Vbase;
+  long i, l, tx = typ(perm);
+
+  if (tx == t_INT) return vectbase;
+  /* old format */
+  l = lg(vectbase); Vbase = cgetg(l,t_VEC);
+  for (i=1; i<l; i++) Vbase[i] = vectbase[itos((GEN)perm[i])];
+  return Vbase;
+}
+
+/* all primes up to Minkowski bound factor on factorbase ? */
+void
+testprimes(GEN bnf, long bound)
+{
+  gpmem_t av0 = avma, av;
+  long p,i,nbideal,k,pmax;
+  GEN f, dK, p1, Vbase, vP, fb, nf=checknf(bnf);
+  byteptr d = diffptr;
+  FB_t F;
+
+  if (DEBUGLEVEL>1)
+    fprintferr("PHASE 1: check primes to Zimmert bound = %ld\n\n",bound);
+  dK= (GEN)nf[3];
+  f = (GEN)nf[4];
+  if (!gcmp1(f))
   {
-    long k = vinvperm[primfact[i]];
-    long l = k - colW;
-    if (l <= 0) xalpha[k]=lstoi(expoprimfact[i]);
-    else        yalpha[l]=lstoi(expoprimfact[i]);
+    GEN D = gmael(nf,5,5);
+    if (DEBUGLEVEL>1) fprintferr("**** Testing Different = %Z\n",D);
+    p1 = isprincipalall(bnf, D, nf_FORCE);
+    if (DEBUGLEVEL>1) fprintferr("     is %Z\n", p1);
   }
+  /* sort factorbase for tablesearch */
+  fb = gen_sort((GEN)bnf[5], 0, &cmp_prime_ideal);
+  p1 = gmael(fb, lg(fb)-1, 1); /* largest p in factorbase */
+  pmax = is_bigint(p1)? VERYBIGINT: itos(p1);
+  if ((ulong)bound > maxprime()) err(primer1);
+  Vbase = get_Vbase(bnf);
+  (void)recover_partFB(&F, Vbase, degpol(nf[1]));
+  av = avma;
+  for (p = 0; p < bound; )
+  {
+    NEXT_PRIME_VIADIFF(p, d);
+    if (DEBUGLEVEL>1) fprintferr("*** p = %ld\n",p);
+    vP = primedec(bnf, stoi(p));
+    nbideal = lg(vP)-1;
+    /* loop through all P | p if ramified, all but one otherwise */
+    if (!smodis(dK,p)) nbideal++;
+    for (i=1; i<nbideal; i++)
+    {
+      GEN P = (GEN)vP[i];
+      if (DEBUGLEVEL>1) fprintferr("  Testing P = %Z\n",P);
+      if (cmpis(idealnorm(bnf,P), bound) >= 0)
+      {
+        if (DEBUGLEVEL>1) fprintferr("    Norm(P) > Zimmert bound\n");
+        continue;
+      }
+      if (p <= pmax && (k = tablesearch(fb, P, &cmp_prime_ideal)))
+      { if (DEBUGLEVEL>1) fprintferr("    #%ld in factor base\n",k); }
+      else if (DEBUGLEVEL>1)
+        fprintferr("    is %Z\n", isprincipal(bnf,P));
+      else /* faster: don't compute result */
+        SPLIT(&F, nf, prime_to_ideal(nf,P), Vbase);
+    }
+    avma = av;
+  }
+  if (DEBUGLEVEL>1) { fprintferr("End of PHASE 1.\n\n"); flusherr(); }
+  avma = av0;
 }
 
 GEN
@@ -897,7 +1010,7 @@ red_mod_units(GEN col, GEN z, long prec)
   RU = lg(mat); x = cgetg(RU+1,t_COL);
   for (i=1; i<RU; i++) x[i]=lreal((GEN)col[i]);
   x[RU] = (long)N2;
-  x = lllintern(concatsp(mat,x), 1,prec);
+  x = lllintern(concatsp(mat,x),100, 1,prec);
   if (!x) return NULL;
   x = (GEN)x[RU];
   if (signe(x[RU]) < 0) x = gneg_i(x);
@@ -939,21 +1052,28 @@ get_Garch(GEN nf, GEN gen, GEN clg2, long prec)
 static GEN
 act_arch(GEN A, GEN x)
 {
-  GEN z, a;
-  long i,l = lg(A);
-  if (typ(A) == t_MAT)
+  GEN a;
+  long i,l = lg(A), tA = typ(A);
+  if (tA == t_MAT)
   {
     /* assume lg(x) >= l */
-    z = cgetg(l, t_VEC);
-    for (i=1; i<l; i++) z[i] = (long)act_arch((GEN)A[i], x);
-    return z;
+    a = cgetg(l, t_VEC);
+    for (i=1; i<l; i++) a[i] = (long)act_arch((GEN)A[i], x);
+    return a;
   }
-  /* A a t_COL of t_INT. Assume lg(A)==lg(x) */
-  l = lg(A); z = cgetg(l, t_VEC);
-  if (l==1) return z;
-  a = gmul((GEN)A[1], (GEN)x[1]);
-  for (i=2; i<l; i++)
-    if (signe(A[i])) a = gadd(a, gmul((GEN)A[i], (GEN)x[i]));
+  if (l==1) return cgetg(1, t_VEC);
+  if (tA == t_VECSMALL)
+  {
+    a = gmulsg(A[1], (GEN)x[1]);
+    for (i=2; i<l; i++)
+      if (A[i]) a = gadd(a, gmulsg(A[i], (GEN)x[i]));
+  }
+  else
+  { /* A a t_COL of t_INT. Assume lg(A)==lg(x) */
+    a = gmul((GEN)A[1], (GEN)x[1]);
+    for (i=2; i<l; i++)
+      if (signe(A[i])) a = gadd(a, gmul((GEN)A[i], (GEN)x[i]));
+  }
   settyp(a, t_VEC); return a;
 }
 
@@ -981,7 +1101,7 @@ isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN 
   N = degpol(nf[1]);
   R1 = itos(gmael(nf,2,1));
   RU = (N + R1)>>1;
-  col = cleancol(col,N,prec); settyp(col, t_COL);
+  col = cleanarch(col,N,prec); settyp(col, t_COL);
   if (RU > 1)
   { /* reduce mod units */
     GEN u, z = init_red_mod_units(bnf,prec);
@@ -1001,59 +1121,64 @@ isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN 
 static int
 fact_ok(GEN nf, GEN y, GEN C, GEN g, GEN e)
 {
+  gpmem_t av = avma;
   long i, c = lg(e);
   GEN z = C? C: gun;
   for (i=1; i<c; i++)
     if (signe(e[i])) z = idealmul(nf, z, idealpow(nf, (GEN)g[i], (GEN)e[i]));
   if (typ(z) != t_MAT) z = idealhermite(nf,z);
   if (typ(y) != t_MAT) y = idealhermite(nf,y);
-  return gegal(y, z);
+  i = gegal(y, z); avma = av; return i;
 }
 
 /* assume x in HNF. cf class_group_gen for notations */
 static GEN
-isprincipalall0(GEN bnf, GEN x, long *ptprec, long flag)
+_isprincipal(GEN bnf, GEN x, long *ptprec, long flag)
 {
   long i,lW,lB,e,c, prec = *ptprec;
   GEN Q,xar,Wex,Bex,U,y,p1,gen,cyc,xc,ex,d,col,A;
   GEN W       = (GEN)bnf[1];
   GEN B       = (GEN)bnf[2];
   GEN WB_C    = (GEN)bnf[4];
-  GEN vperm   = (GEN)bnf[6];
   GEN nf      = (GEN)bnf[7];
   GEN clg2    = (GEN)bnf[9];
   int old_format = (typ(clg2[2]) == t_MAT);
 
-  U = (GEN)clg2[1];
+  U = (GEN)clg2[1]; if (old_format) U = ginv(U);
   cyc = gmael3(bnf,8,1,2); c = lg(cyc)-1;
   gen = gmael3(bnf,8,1,3);
+  ex = cgetg(c+1,t_COL);
+  if (c == 0 && !(flag & (nf_GEN|nf_GENMAT|nf_GEN_IF_PRINCIPAL))) return ex;
 
-  vectbase = (GEN)bnf[5]; /* needed by factorgensimple */
-
   /* factor x */
-  xc = content(x); if (!gcmp1(xc)) x = gdiv(x,xc);
+  x = Q_primitive_part(x, &xc);
+  xar = split_ideal(nf,x,get_Vbase(bnf));
+  lW = lg(W)-1; Wex = vecsmall_const(lW, 0);
+  lB = lg(B)-1; Bex = vecsmall_const(lB, 0);
+  for (i=1; i<=primfact[0]; i++)
+  {
+    long k = primfact[i];
+    long l = k - lW;
+    if (l <= 0) Wex[k] = exprimfact[i];
+    else        Bex[l] = exprimfact[i];
+  }
 
-  p1 = init_idele(nf); p1[1] = (long)x;
-  if (!(flag & nf_GEN)) p1[2] = 0; /* don't compute archimedean part */
-  xar = split_ideal(nf,p1,prec,vperm);
-
-  lW = lg(W)-1; Wex = zerocol(lW);
-  lB = lg(B)-1; Bex = zerocol(lB); get_split_expo(Wex,Bex,vperm);
-
   /* x = g_W Wex + g_B Bex - [xar]
    *   = g_W A - [xar] + [C_B]Bex  since g_W B + g_B = [C_B] */
-  A = gsub(Wex, gmul(B,Bex));
-  if (old_format) U = ginv(U);
+  A = gsub(small_to_col(Wex), gmul_mati_smallvec(B,Bex));
   Q = gmul(U, A);
-  ex = cgetg(c+1,t_COL);
   for (i=1; i<=c; i++)
-    Q[i] = (long)truedvmdii((GEN)Q[i],(GEN)cyc[i],(GEN*)(ex+i));
-  if (!(flag & nf_GEN)) return gcopy(ex);
+    Q[i] = (long)truedvmdii((GEN)Q[i], (GEN)cyc[i], (GEN*)(ex+i));
+  if ((flag & nf_GEN_IF_PRINCIPAL))
+    { if (!gcmp0(ex)) return gzero; }
+  else if (!(flag & (nf_GEN|nf_GENMAT)))
+    return gcopy(ex);
 
   /* compute arch component of the missing principal ideal */
   if (old_format)
   {
     GEN Garch, V = (GEN)clg2[2];
+    Bex = small_to_col(Bex);
     p1 = c? concatsp(gmul(V,Q), Bex): Bex;
     col = act_arch(p1, WB_C);
     if (c)
@@ -1066,16 +1191,22 @@ isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla
   { /* g A = G Ur A + [ga]A, Ur A = D Q + R as above (R = ex)
            = G R + [GD]Q + [ga]A */
     GEN ga = (GEN)clg2[2], GD = (GEN)clg2[3];
-    col = act_arch(Bex, WB_C + lW);
+    if (lB) col = act_arch(Bex, WB_C + lW); else col = zerovec(1); /* nf=Q ! */
     if (lW) col = gadd(col, act_arch(A, ga));
     if (c)  col = gadd(col, act_arch(Q, GD));
   }
-  col = gsub(col, (GEN)xar[2]);
+  if (xar) col = gadd(col, famat_to_arch(nf, xar, prec));
 
   /* find coords on Zk; Q = N (x / \prod gj^ej) = N(alpha), denom(alpha) | d */
   Q = gdiv(dethnf_i(x), get_norm_fact(gen, ex, &d));
   col = isprincipalarch(bnf, col, Q, gun, d, &e);
   if (col && !fact_ok(nf,x, col,gen,ex)) col = NULL;
+  if (!col && !gcmp0(ex))
+  {
+    p1 = isprincipalfact(bnf, gen, gneg(ex), x, flag);
+    col = (GEN)p1[2];
+    e = itos((GEN)p1[3]);
+  }
   if (!col)
   {
     *ptprec = prec + (e >> TWOPOTBITS_IN_LONG) + (MEDDEFAULTPREC-2);
@@ -1088,8 +1219,11 @@ isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla
     e = 0;
   }
   y = cgetg(4,t_VEC);
+  if (!xc) xc = gun;
+  col = e? gmul(xc,col): cgetg(1,t_COL);
+  if (flag & nf_GEN_IF_PRINCIPAL) return col;
   y[1] = lcopy(ex);
-  y[2] = e? lmul(xc,col): lgetg(1,t_COL);
+  y[2] = (long)col;
   y[3] = lstoi(-e); return y;
 }
 
@@ -1097,18 +1231,20 @@ static GEN
 triv_gen(GEN nf, GEN x, long c, long flag)
 {
   GEN y;
-  if (!(flag & nf_GEN)) return zerocol(c);
+  if (flag & nf_GEN_IF_PRINCIPAL)
+    return (typ(x) == t_COL)? gcopy(x): algtobasis(nf,x);
+  if (!(flag & (nf_GEN|nf_GENMAT))) return zerocol(c);
   y = cgetg(4,t_VEC);
   y[1] = (long)zerocol(c);
-  x = (typ(x) == t_COL)? gcopy(x): algtobasis(nf,x);
-  y[2] = (long)x;
+  y[2] = (long)((typ(x) == t_COL)? gcopy(x): 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);
+  long c, pr, tx = typ(x);
+  gpmem_t av = avma;
   GEN nf;
 
   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
@@ -1118,26 +1254,26 @@ isprincipalall(GEN bnf,GEN x,long flag)
       err(talker,"not the same number field in isprincipal");
     x = (GEN)x[2]; tx = t_POL;
   }
-  if (tx == t_POL || tx == t_COL)
+  if (tx == t_POL || tx == t_COL || tx == t_INT || tx == t_FRAC)
   {
     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)
+  if (degpol(nf[1]) == 1)
     return gerepileupto(av, triv_gen(nf, gcoeff(x,1,1), 0, flag));
 
   pr = prec_arch(bnf); /* precision of unit matrix */
   c = getrand();
   for (;;)
   {
-    long av1 = avma;
-    GEN y = isprincipalall0(bnf,x,&pr,flag);
+    gpmem_t av1 = avma;
+    GEN y = _isprincipal(bnf,x,&pr,flag);
     if (y) return gerepileupto(av,y);
 
-    if (DEBUGLEVEL) err(warnprec,"isprincipalall0",pr);
-    avma = av1; bnf = bnfnewprec(bnf,pr); setrand(c);
+    if (DEBUGLEVEL) err(warnprec,"isprincipal",pr);
+    avma = av1; bnf = bnfnewprec(bnf,pr); (void)setrand(c);
   }
 }
 
@@ -1145,9 +1281,10 @@ isprincipalall(GEN bnf,GEN x,long flag)
 GEN
 isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag)
 {
-  long av = avma, l = lg(e), i,prec,c;
+  long l = lg(e), i, prec, c;
+  gpmem_t av = avma;
   GEN id,id2, nf = checknf(bnf), z = NULL; /* gcc -Wall */
-  int gen = flag & (nf_GEN | nf_GENMAT);
+  int gen = flag & (nf_GEN|nf_GENMAT);
 
   prec = prec_arch(bnf);
   if (gen)
@@ -1163,28 +1300,30 @@ isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag
       id2 = idealpowred(bnf,z, (GEN)e[i],prec);
       id = id? idealmulred(nf,id,id2,prec): id2;
     }
-  if (id == C)
+  if (id == C) /* e = 0 */
   {
-    if (!C) id = gun;
-    return isprincipalall(bnf, id, flag);
+    if (!C) return isprincipalall(bnf, gun, flag);
+    C = idealhermite(nf,C); id = z;
+    if (gen) id[1] = (long)C; else id = C;
   }
   c = getrand();
   for (;;)
   {
-    long av1 = avma;
-    GEN y = isprincipalall0(bnf, gen? (GEN)id[1]: id,&prec,flag);
+    gpmem_t av1 = avma;
+    GEN y = _isprincipal(bnf, gen? (GEN)id[1]: id,&prec,flag);
     if (y)
     {
-      if (gen && typ(y) == t_VEC)
+      GEN u = (GEN)y[2];
+      if (!gen || typ(y) != t_VEC) return gerepileupto(av,y);
+      if (flag & nf_GENMAT)
+        y[2] = (isnfscalar(u) && gcmp1((GEN)u[1]))? id[2]
+                                                  :(long)arch_mul((GEN)id[2],u);
+      else
       {
-        GEN u = lift_intern(basistoalg(nf, (GEN)y[2]));
-        if (flag & nf_GENMAT)
-          y[2] = (gcmp1(u)&&lg(id[2])>1)? id[2]: (long)arch_mul((GEN)id[2], u);
-        else
-          y[2] = (long)algtobasis(nf, gmul((GEN)id[2], u));
-        y = gcopy(y);
+        u = lift_intern(basistoalg(nf, u));
+        y[2] = (long)algtobasis(nf, gmul((GEN)id[2], u));
       }
-      return gerepileupto(av,y);
+      return gerepilecopy(av, y);
     }
 
     if (flag & nf_GIVEPREC)
@@ -1193,15 +1332,15 @@ isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag
         err(warner,"insufficient precision for generators, not given");
       avma = av; return stoi(prec);
     }
-    if (DEBUGLEVEL) err(warnprec,"isprincipalall0",prec);
-    avma = av1; bnf = bnfnewprec(bnf,prec); setrand(c);
+    if (DEBUGLEVEL) err(warnprec,"isprincipal",prec);
+    avma = av1; bnf = bnfnewprec(bnf,prec); (void)setrand(c);
   }
 }
 
 GEN
 isprincipal(GEN bnf,GEN x)
 {
-  return isprincipalall(bnf,x,nf_REGULAR);
+  return isprincipalall(bnf,x,0);
 }
 
 GEN
@@ -1222,91 +1361,114 @@ isprincipalgenforce(GEN bnf,GEN x)
   return isprincipalall(bnf,x,nf_GEN | nf_FORCE);
 }
 
+static GEN
+rational_unit(GEN x, long n, long RU)
+{
+  GEN y;
+  if (!gcmp1(x) && !gcmp_1(x)) return cgetg(1,t_COL);
+  y = zerocol(RU); 
+  y[RU] = (long)gmodulss((gsigne(x)>0)? 0: n>>1, n);
+  return y;
+}
+
+/* if x a famat, assume it is an algebraic integer (very costly to check) */
 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;
+  long tx = typ(x), i, R1, RU, n, prec;
+  gpmem_t av = avma;
+  GEN p1, v, rlog, logunit, ex, nf, z, pi2_sur_w, gn, emb;
 
   bnf = checkbnf(bnf); nf=(GEN)bnf[7];
-  logunit=(GEN)bnf[3]; RU=lg(logunit);
+  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];
+  z  = algtobasis(nf, (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;
+      return rational_unit(x, n, RU);
 
-    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_MAT: /* famat */
+      if (lg(x) != 3 || lg(x[1]) != lg(x[2]))
+        err(talker, "not a factorization matrix in isunit");
+      break;
 
     case t_COL:
-      if (lgef(nf[1])-2 == lg(x)) { p1 = basistoalg(nf,x); break; }
+      if (degpol(nf[1]) != lg(x)-1)
+        err(talker,"not an algebraic number in isunit");
+      break;
 
-    default:
-      err(talker,"not an algebraic number in isunit");
+    default: x = algtobasis(nf, x);
+      break;
   }
-  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); }
+  /* assume a famat is integral */
+  if (tx != t_MAT && !gcmp1(denom(x))) { avma = av; return cgetg(1,t_COL); }
+  if (isnfscalar(x)) return gerepileupto(av, rational_unit((GEN)x[1],n,RU));
 
-  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);
+  R1 = nf_get_r1(nf); v = cgetg(RU+1,t_COL);
+  for (i=1; i<=R1; i++) v[i] = un;
+  for (   ; i<=RU; i++) v[i] = deux;
+  logunit = concatsp(logunit, v);
   /* ex = fundamental units exponents */
+  rlog = greal(logunit);
+  prec = nfgetprec(nf);
+  for (i=1;;)
   {
-    GEN rx, rlog = greal(logunit);
-    long e, prec = nfgetprec(nf);
-    for (i=1;;)
+    GEN logN, rx = get_arch_real(nf,x,&emb, MEDDEFAULTPREC);
+    long e;
+    if (rx)
     {
-      rx = get_arch_real(nf,x,&emb, MEDDEFAULTPREC);
-      if (rx)
+      logN = sum(rx, 1, RU); /* log(Nx), should be ~ 0 */
+      if (gexpo(logN) > -20)
       {
-        ex = grndtoi(gauss(rlog, rx), &e);
-        if (gcmp0((GEN)ex[RU]) && e < -4) break;
+        long p = 2 + max(1, (nfgetprec(nf)-2) / 2);
+        if (typ(logN) != t_REAL || gprecision(rx) > p)
+          { avma = av; return cgetg(1,t_COL); } /* not a precision problem */
+        rx = NULL;
       }
-
-      if (++i > 4) err(precer,"isunit");
+    }
+    if (rx)
+    {
+      ex = grndtoi(gauss(rlog, rx), &e);
+      if (gcmp0((GEN)ex[RU]) && e < -4) break;
+    }
+    if (i == 1)
+      prec = MEDDEFAULTPREC + (gexpo(x) >> TWOPOTBITS_IN_LONG);
+    else
+    {
+      if (i > 4) err(precer,"isunit");
       prec = (prec-1)<<1;
-      if (DEBUGLEVEL) err(warnprec,"isunit",prec);
-      nf = nfnewprec(nf, prec);
     }
+    i++;
+    if (DEBUGLEVEL) err(warnprec,"isunit",prec);
+    nf = nfnewprec(nf, prec);
   }
 
   setlg(ex, RU);
-  setlg(p1, RU); settyp(p1, t_VEC);
-  for (i=1; i<RU; i++) p1[i] = coeff(logunit, 1, i);
+  p1 = row_i(logunit,1, 1,RU-1);
   p1 = gneg(gimag(gmul(p1,ex))); if (!R1) p1 = gmul2n(p1, -1);
-  p1 = gadd(garg((GEN)emb[1],DEFAULTPREC), p1);
+  p1 = gadd(garg((GEN)emb[1],prec), p1);
   /* p1 = arg(the missing root of 1) */
 
-  pi2_sur_w = divrs(mppi(DEFAULTPREC), n>>1);
+  pi2_sur_w = divrs(mppi(prec), n>>1); /* 2pi / n */
   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));
+    GEN ro = gmul(row(gmael(nf,5,1), 1), z);
+    GEN p2 = ground(gdiv(garg(ro, prec), 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);
+  ex[RU] = lmodulcp(p1, gn);
+  setlg(ex, RU+1); return gerepilecopy(av, ex);
 }
 
 GEN
 signunits(GEN bnf)
 {
-  long av,i,j,R1,RU,mun;
+  long i, j, R1, RU, mun;
+  gpmem_t av;
   GEN matunit,y,p1,p2,nf,pi;
 
   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
@@ -1326,57 +1488,30 @@ signunits(GEN bnf)
   return y;
 }
 
-/* LLL-reduce ideal and return T2 | ideal */
+/* LLL-reduce ideal and return Cholesky for T2 | ideal */
 static GEN
-red_ideal(GEN *ideal,GEN T2vec,GEN prvec)
+red_ideal(GEN *ideal,GEN Gvec,GEN prvec)
 {
   long i;
-  for (i=1; i<lg(T2vec); i++)
+  for (i=1; i<lg(Gvec); i++)
   {
-    GEN p1,base, T2 = (GEN)T2vec[i];
+    GEN p1,u, G = (GEN)Gvec[i];
     long prec = prvec[i];
 
-    p1 = qf_base_change(T2, *ideal, 1);
-    base = lllgramintern(p1,100,1,prec);
-    if (base)
+    p1 = gmul(G, *ideal);
+    u = lllintern(p1,100,1,prec);
+    if (u)
     {
-      p1 = sqred1intern(qf_base_change(p1,base,1),1);
-      if (p1) { *ideal = gmul(*ideal,base); return p1; }
+      p1 = gmul(p1, u);
+      p1 = sqred1_from_QR(p1, prec);
+      if (p1) { *ideal = gmul(*ideal,u); return p1; }
     }
     if (DEBUGLEVEL) err(warner, "prec too low in red_ideal[%ld]: %ld",i,prec);
   }
   return NULL;
 }
 
-static double
-get_minkovski(long N, long R1, GEN D, GEN gborne)
-{
-  const long prec = DEFAULTPREC;
-  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 *= 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 void
 set_log_embed(GEN col, GEN x, long R1, long prec)
 {
   long i, l = lg(x);
@@ -1385,11 +1520,11 @@ set_log_embed(GEN col, GEN x, long R1, long prec)
 }
 
 static void
-set_fact(GEN c)
+set_fact(GEN c, GEN first_nz, long cglob)
 {
   long i;
-  c[0] = primfact[1]; /* for already_found_relation */
-  for (i=1; i<=primfact[0]; i++) c[primfact[i]] = expoprimfact[i];
+  first_nz[cglob] = primfact[1];
+  for (i=1; i<=primfact[0]; i++) c[primfact[i]] = exprimfact[i];
 }
 
 static void
@@ -1399,73 +1534,81 @@ unset_fact(GEN c)
   for (i=1; i<=primfact[0]; i++) c[primfact[i]] = 0;
 }
 
-#define MAXTRY 20
+/* as small_to_mat, with explicit dimension d */
+static GEN
+small_to_mat_i(GEN z, long d)
+{
+  long i,j, l = lg(z);
+  GEN x = cgetg(l,t_MAT);
+  for (j=1; j<l; j++)
+  {
+    GEN c = cgetg(d+1, t_COL), cz = (GEN)z[j];
+    x[j] = (long)c;
+    for (i=1; i<=d; i++) c[i] = lstoi(cz[i]);
+  }
+  return x;
+}
 
-/* return -1 in case of precision problems. t = current number of relations */
+/* return -1 in case of precision problems. t = current # of relations */
 static long
-small_norm_for_buchall(long cglob,GEN *mat,GEN matarch,long LIMC, long PRECREG,
-                       GEN nf,GEN gborne,long nbrelpid,GEN invp,GEN L,GEN LLLnf)
+small_norm_for_buchall(long cglob,GEN *mat,GEN first_nz, GEN matarch,
+                       long LIMC, long PRECREG, FB_t *F,
+                       GEN nf,long nbrelpid,GEN invp,GEN L)
 {
+  const int maxtry_DEP  = 20;
+  const int maxtry_FACT = 500;
   const double eps = 0.000001;
-  double *y,*z,**q,*v, MINKOVSKI_BOUND,BOUND;
-  ulong av = avma, av1,av2,limpile;
-  long i,j,k,noideal, nbrel = lg(mat)-1;
-  long alldep = 0, nbsmallnorm,nbfact,R1, N = degpol(nf[1]);
-  GEN x,xembed,M,T2,r,cbase,invcbase,T2vec,prvec;
+  double *y,*z,**q,*v, BOUND;
+  gpmem_t av = avma, av1, av2, limpile;
+  long j,k,noideal, nbrel = lg(mat)-1;
+  long nbsmallnorm,nbfact,R1, N = degpol(nf[1]);
+  GEN x,xembed,M,G,r,Gvec,prvec;
 
-  if (gsigne(gborne)<=0) return cglob;
   if (DEBUGLEVEL)
     fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel);
   xembed = NULL; /* gcc -Wall */
   nbsmallnorm = nbfact = 0;
-  R1 = itos(gmael(nf,2,1));
-  T2 = (GEN)LLLnf[1];
-  M  = (GEN)LLLnf[2];
-  cbase   =(GEN)LLLnf[3];
-  invcbase=(GEN)LLLnf[4];
+  R1 = nf_get_r1(nf);
+  M = gmael(nf,5,1);
+  G = gmael(nf,5,2);
 
-  prvec = cgetg(3,t_VECSMALL);
-  prvec[1] = MEDDEFAULTPREC;
-  prvec[2] = PRECREG;
-  T2vec = cgetg(3,t_VEC);
-  T2vec[1] = (long)gprec_w(T2,prvec[1]);
-  T2vec[2] = (long)T2;
+  prvec = cgetg(3,t_VECSMALL); Gvec = cgetg(3,t_VEC);
+  prvec[1] = MEDDEFAULTPREC;   Gvec[1] = (long)gprec_w(G,prvec[1]);
+  prvec[2] = PRECREG;          Gvec[2] = (long)G;
   minim_alloc(N+1, &q, &x, &y, &z, &v);
   av1 = avma;
-  MINKOVSKI_BOUND = get_minkovski(N,R1,(GEN)nf[3],gborne);
-  for (noideal=KC; noideal; noideal--)
+  for (noideal=F->KC; noideal; noideal--)
   {
-    long nbrelideal=0, dependent = 0;
-    GEN IDEAL, ideal = (GEN)vectbase[noideal];
-    GEN normideal = idealnorm(nf,ideal);
+    gpmem_t av0 = avma;
+    long nbrelideal=0, dependent = 0, try_factor = 0, oldcglob = cglob;
+    GEN IDEAL, ideal = (GEN)F->LP[noideal];
 
-    if (alldep > 2*N) break;
-
     if (DEBUGLEVEL>1) fprintferr("\n*** Ideal no %ld: %Z\n", noideal, ideal);
     ideal = prime_to_ideal(nf,ideal);
-    IDEAL = invcbase? gmul(invcbase, ideal): ideal;
-    IDEAL = gmul(IDEAL, lllint(IDEAL)); /* should be almost T2-reduced */
-    r = red_ideal(&IDEAL,T2vec,prvec);
+    IDEAL = lllint_ip(ideal,4); /* should be almost T2-reduced */
+    r = red_ideal(&IDEAL,Gvec,prvec);
     if (!r) return -1; /* precision problem */
 
     for (k=1; k<=N; k++)
     {
       v[k]=gtodouble(gcoeff(r,k,k));
       for (j=1; j<k; j++) q[j][k]=gtodouble(gcoeff(r,j,k));
-      if (DEBUGLEVEL>3) fprintferr("v[%ld]=%.0f ",k,v[k]);
+      if (DEBUGLEVEL>3) fprintferr("v[%ld]=%.4g ",k,v[k]);
     }
 
-    BOUND = MINKOVSKI_BOUND * pow(gtodouble(normideal),2./(double)N);
+    BOUND = v[2] + v[1] * q[1][2] * q[1][2];
+    if (BOUND < v[1]) BOUND = v[1];
+    BOUND *= 2; /* at most twice larger than smallest known vector */
     if (DEBUGLEVEL>1)
     {
       if (DEBUGLEVEL>3) fprintferr("\n");
-      fprintferr("BOUND = %.0f\n",BOUND); flusherr();
+      fprintferr("BOUND = %.4g\n",BOUND); flusherr();
     }
-    BOUND += eps; av2=avma; limpile = stack_lim(av2,1);
+    BOUND *= 1 + eps; av2=avma; limpile = stack_lim(av2,1);
     k = N; y[N]=z[N]=0; x[N]= (long) sqrt(BOUND/v[N]);
     for(;; x[1]--)
     {
-      ulong av3 = avma;
+      gpmem_t av3 = avma;
       double p;
       GEN col;
 
@@ -1489,16 +1632,17 @@ small_norm_for_buchall(long cglob,GEN *mat,GEN matarch
 	}
 	if (k==1) /* element complete */
 	{
-	  if (!x[1] && y[1]<=eps) goto ENDIDEAL;
+	  if (y[1]<=eps) goto ENDIDEAL; /* skip all scalars: [*,0...0] */
 	  if (ccontent(x)==1) /* primitive */
 	  {
             GEN Nx, gx = gmul_mati_smallvec(IDEAL,x);
-            long av4;
+            gpmem_t av4;
             if (!isnfscalar(gx))
             {
               xembed = gmul(M,gx); av4 = avma; nbsmallnorm++;
+              if (++try_factor > maxtry_FACT) goto ENDIDEAL;
               Nx = ground(norm_by_embed(R1,xembed)); setsigne(Nx, 1);
-              if (factorelt(nf,cbase,gx,Nx,KCZ,LIMC)) { avma = av4; break; }
+              if (can_factor(F, nf, NULL, gx, Nx)) { avma = av4; break; }
               if (DEBUGLEVEL > 1) { fprintferr("."); flusherr(); }
             }
 	    avma = av3;
@@ -1507,24 +1651,20 @@ small_norm_for_buchall(long cglob,GEN *mat,GEN matarch
 	}
       }
       cglob++; col = mat[cglob];
-      set_fact(col);
-      if (cglob > 1 && cglob <= KC && ! addcolumntomatrix(col,invp,L))
+      set_fact(col, first_nz, cglob);
+      /* make sure we get maximal rank first, then allow all relations */
+      if (cglob > 1 && cglob <= F->KC && ! addcolumntomatrix(col,invp,L))
       { /* Q-dependent from previous ones: forget it */
         cglob--; unset_fact(col);
         if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
-        if (++dependent > MAXTRY) { alldep++; break; }
+        if (++dependent > maxtry_DEP) break;
         avma = av3; continue;
       }
       /* record archimedean part */
       set_log_embed((GEN)matarch[cglob], xembed, R1, PRECREG);
-      alldep = dependent = 0;
+      dependent = 0;
 
-      if (DEBUGLEVEL)
-      {
-        if (DEBUGLEVEL==1) fprintferr("%4ld",cglob);
-        else { fprintferr("cglob = %ld. ",cglob); wr_rel(mat[cglob]); }
-        flusherr(); nbfact++;
-      }
+      if (DEBUGLEVEL) { nbfact++; dbg_rel(cglob, mat[cglob]); }
       if (cglob >= nbrel) goto END; /* we have enough */
       if (++nbrelideal == nbrelpid) break;
 
@@ -1535,56 +1675,34 @@ small_norm_for_buchall(long cglob,GEN *mat,GEN matarch
       }
     }
 ENDIDEAL:
-    invp = gerepilecopy(av1, invp);
+    if (cglob == oldcglob) avma = av0; else invp = gerepilecopy(av1, invp);
     if (DEBUGLEVEL>1) msgtimer("for this ideal");
   }
 END:
   if (DEBUGLEVEL)
   {
-    fprintferr("\n");
-    msgtimer("small norm relations");
-    if (DEBUGLEVEL>1)
-    {
-      GEN p1,tmp=cgetg(cglob+1,t_MAT);
-
-      fprintferr("Elements of small norm gave %ld relations.\n",cglob);
-      fprintferr("Computing rank: "); flusherr();
-      for(j=1;j<=cglob;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)sindexrank(tmp)[2]; k=lg(tmp)-1;
-      fprintferr("%ld; independent columns: %Z\n",k,tmp);
-    }
-    if(nbsmallnorm)
-      fprintferr("\nnb. fact./nb. small norm = %ld/%ld = %.3f\n",
+    fprintferr("\n"); msgtimer("small norm relations");
+    fprintferr("  small norms gave %ld relations, rank = %ld.\n",
+               cglob, rank(small_to_mat_i((GEN)mat, F->KC)));
+    if (nbsmallnorm)
+      fprintferr("  nb. fact./nb. small norm = %ld/%ld = %.3f\n",
                   nbfact,nbsmallnorm,((double)nbfact)/nbsmallnorm);
-    else
-      fprintferr("\nnb. small norm = 0\n");
   }
   avma = av; return cglob;
 }
-#undef MAXTRY
 
-/* I assumed to be integral HNF, T2 a weighted T2 matrix */
+/* I assumed to be integral HNF, G the Cholesky form of a weighted T2 matrix.
+ * Return an irrational m in I with T2(m) small */
 static GEN
-ideallllredpart1(GEN I, GEN T2, long prec)
+pseudomin(GEN I, GEN G)
 {
-  GEN y,m,idealpro;
-
-  y = lllgramintern(qf_base_change(T2,I,1),100,1,prec+1);
+  GEN m, y = lllintern(gmul(G, I), 100,1, 0);
   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(3,t_VEC);
-  idealpro[1] = (long)I;
-  idealpro[2] = (long)m; /* irrational element of small T2 norm in I */
-  if (DEBUGLEVEL>5) fprintferr("\nidealpro = %Z\n",idealpro);
-  return idealpro;
+  if (DEBUGLEVEL>5) fprintferr("\nm = %Z\n",m);
+  return m;
 }
 
 static void
@@ -1605,70 +1723,52 @@ dbg_newrel(long jideal, long jdir, long phase, long cg
 static void
 dbg_cancelrel(long jideal,long jdir,long phase, long *col)
 {
-  fprintferr("rel. cancelled. phase %ld: %ld ",phase);
+  fprintferr("rel. cancelled. phase %ld: ",phase);
   if (DEBUGLEVEL>3) fprintferr("(jideal=%ld,jdir=%ld)",jideal,jdir);
   wr_rel(col); flusherr();
 }
 
 static void
-dbg_outrel(long phase,long cglob, GEN vperm,GEN *mat,GEN maarch)
+dbg_outrel(long cglob, GEN *mat,GEN maarch)
 {
-  long i,j;
-  GEN p1,p2;
+  gpmem_t av = avma;
+  GEN p1;
+  long j;
 
-  if (phase == 0)
+  p1 = cgetg(cglob+1, t_MAT);
+  for (j=1; j<=cglob; j++) p1[j] = (long)small_to_col(mat[j]);
+  fprintferr("\nRank = %ld\n", rank(p1));
+  if (DEBUGLEVEL>3)
   {
-    ulong av = avma; p2=cgetg(cglob+1,t_MAT);
-    for (j=1; j<=cglob; j++)
-    {
-      p1=cgetg(KC+1,t_COL); p2[j]=(long)p1;
-      for (i=1; i<=KC; i++) p1[i]=lstoi(mat[j][i]);
-    }
-    fprintferr("\nRank = %ld\n", rank(p2));
-    if (DEBUGLEVEL>3)
-    {
-      fprintferr("relations = \n");
-      for (j=1; j <= cglob; j++) wr_rel(mat[j]);
-      fprintferr("\nmatarch = %Z\n",maarch);
-    }
-    avma = av;
+    fprintferr("relations = \n");
+    for (j=1; j <= cglob; j++) wr_rel(mat[j]);
+    fprintferr("\nmatarch = %Z\n",maarch);
   }
-  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();
+  avma = av; flusherr();
 }
 
-/* Check if we already have a column mat[l] equal to mat[s]
+/* Check if we already have a column mat[i] equal to mat[s]
  * General check for colinearity useless since exceedingly rare */
 static long
-already_found_relation(GEN *mat,long s)
+already_found_relation(GEN *mat, long s, GEN first_nz)
 {
-  GEN coll, cols = mat[s];
-  long l,bs;
+  GEN cols = mat[s];
+  long i, bs, l = lg(cols);
 
-  bs = 1; while (bs<=KC && !cols[bs]) bs++;
-  if (bs > KC) return s; /* zero relation */
+  bs = 1; while (bs < l && !cols[bs]) bs++;
+  if (bs == l) return s; /* zero relation */
 
-  for (l=s-1; l; l--)
+  for (i=s-1; i; i--)
   {
-    coll = mat[l];
-    if (bs == coll[0]) /* = index of first non zero elt in coll */
+    if (bs == first_nz[i]) /* = index of first non zero elt in cols */
     {
+      GEN coll = mat[i];
       long b = bs;
-      do b++; while (b<=KC && cols[b] == coll[b]);
-      if (b > KC) return l;
+      do b++; while (b < l && cols[b] == coll[b]);
+      if (b == l) return i;
     }
   }
-  cols[0] = bs; return 0;
+  first_nz[s] = bs; return 0;
 }
 
 /* I integral ideal in HNF form */
@@ -1676,71 +1776,79 @@ static GEN
 remove_content(GEN I)
 {
   long N = lg(I)-1;
-  if (!gcmp1(gcoeff(I,N,N))) { GEN y=content(I); if (!gcmp1(y)) I=gdiv(I,y); }
+  if (!gcmp1(gcoeff(I,N,N))) I = Q_primpart(I);
   return I;
 }
 
 /* if phase != 1 re-initialize static variables. If <0 return immediately */
 static long
-random_relation(long phase,long cglob,long LIMC,long PRECLLL,long PRECREG,
-                GEN nf,GEN subFB,GEN vecT2,GEN *mat,GEN matarch,GEN list_jideal)
+random_relation(long phase,long cglob,long LIMC,long PRECREG,long MAXRELSUP,
+                GEN nf,GEN vecG,GEN *mat,GEN first_nz,GEN matarch,
+                GEN list_jideal, FB_t *F)
 {
   static long jideal, jdir;
-  long lim,i,av,av1,cptzer,nbT2,lgsub,r1, jlist = 1;
-  GEN arch,col,colarch,ideal,idealpro,P,ex;
+  long i, maxcglob, cptlist, cptzer, nbG, lgsub, r1, jlist = 1;
+  gpmem_t av, av1;
+  GEN arch,col,colarch,ideal,m,P,ex;
 
   if (phase != 1) { jideal=jdir=1; if (phase<0) return 0; }
 
   r1 = nf_get_r1(nf);
-  lim = lg(mat)-1; /* requested number of relations */
-  nbT2 = lg(vecT2)-1;
-  lgsub = lg(subFB); ex = cgetg(lgsub, t_VECSMALL);
-  cptzer = 0;
+  maxcglob = lg(mat)-1; /* requested # of relations */
+  nbG = lg(vecG)-1;
+  lgsub = lg(F->subFB); ex = cgetg(lgsub, t_VECSMALL);
+  cptzer = cptlist = 0;
   if (DEBUGLEVEL && list_jideal)
     fprintferr("looking hard for %Z\n",list_jideal);
   P = NULL; /* gcc -Wall */
   for (av = avma;;)
   {
-    if (list_jideal && jlist < lg(list_jideal) && jdir <= nbT2)
-      jideal = list_jideal[jlist++];
-    if (!list_jideal || jdir <= nbT2)
+    if (list_jideal && jlist < lg(list_jideal) && jdir <= nbG)
     {
+      jideal = list_jideal[jlist++]; cptlist = 0;
+    }
+    if (!list_jideal || jdir <= nbG)
+    {
       avma = av;
-      P = prime_to_ideal(nf, (GEN)vectbase[jideal]);
+      P = prime_to_ideal(nf, (GEN)F->LP[jideal]);
     }
+    else
+    {
+      if (++cptlist > 300) return -1;
+    }
     ideal = P;
     do {
       for (i=1; i<lgsub; i++)
       {
         ex[i] = mymyrand()>>randshift;
-        if (ex[i]) ideal = idealmulh(nf,ideal, gmael3(powsubFB,i,ex[i],1));
+        if (ex[i]) ideal = idealmulh(nf,ideal, gmael3(F->powsubFB,i,ex[i],1));
       }
     }
     while (ideal == P); /* If ex  = 0, try another */
     ideal = remove_content(ideal);
 
     if (phase != 1) jdir = 1; else phase = 2;
-    for (av1 = avma; jdir <= nbT2; jdir++, avma = av1)
+    for (av1 = avma; jdir <= nbG; 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(ideal,(GEN)vecT2[jdir],PRECLLL);
-      if (!idealpro) return -2;
-      if (!factorgen(nf,idealpro,KCZ,LIMC))
+      m = pseudomin(ideal,(GEN)vecG[jdir]);
+      if (!m) return -2;
+      if (!factorgen(F,nf,ideal,m))
       {
         if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
         continue;
       }
       /* can factor ideal, record relation */
       cglob++; col = mat[cglob];
-      set_fact(col); col[jideal]--;
-      for (i=1; i<lgsub; i++) col[subFB[i]] -= ex[i];
-      if (already_found_relation(mat, cglob))
+      set_fact(col, first_nz, cglob); col[jideal]--;
+      for (i=1; i<lgsub; i++) col[ F->subFB[i] ] -= ex[i];
+      if (already_found_relation(mat, cglob, first_nz))
       { /* Already known: forget it */
         if (DEBUGLEVEL>1) dbg_cancelrel(jideal,jdir,phase,col);
         cglob--; unset_fact(col); col[jideal] = 0;
-        for (i=1; i<lgsub; i++) col[subFB[i]] = 0;
+        for (i=1; i<lgsub; i++) col[ F->subFB[i] ] = 0;
 
         if (++cptzer > MAXRELSUP)
         {
@@ -1755,23 +1863,23 @@ random_relation(long phase,long cglob,long LIMC,long P
       for (i=1; i<lgsub; i++)
         if (ex[i])
         {
-          GEN p1 = gmael3(powsubFB,i,ex[i],2);
+          GEN p1 = gmael3(F->powsubFB,i,ex[i],2);
           arch = arch? vecmul(arch, p1): p1;
         }
       colarch = (GEN)matarch[cglob];
       /* arch = archimedean component (MULTIPLICATIVE form) of ideal */
-      arch = vecdiv(arch, gmul(gmael(nf,5,1), (GEN)idealpro[2]));
+      arch = vecdiv(arch, gmul(gmael(nf,5,1), m));
       set_log_embed(colarch, arch, r1, PRECREG);
-      if (DEBUGLEVEL) dbg_newrel(jideal,jdir,phase,cglob,col,colarch,lim);
+      if (DEBUGLEVEL) dbg_newrel(jideal,jdir,phase,cglob,col,colarch,maxcglob);
 
       /* Need more, try next P */
-      if (cglob < lim) break;
+      if (cglob < maxcglob) break;
 
       /* We have found enough. Return */
       if (phase)
       {
         jdir = 1;
-        if (jideal == KC) jideal=1; else jideal++;
+        if (jideal == F->KC) jideal=1; else jideal++;
       }
       if (DEBUGLEVEL>2)
         fprintferr("Upon exit: jideal=%ld,jdir=%ld\n",jideal,jdir);
@@ -1779,61 +1887,69 @@ random_relation(long phase,long cglob,long LIMC,long P
     }
     if (!list_jideal)
     {
-      if (jideal == KC) jideal=1; else jideal++;
+      if (jideal == F->KC) jideal=1; else jideal++;
     }
   }
 }
 
-static long
-be_honest(GEN nf,GEN subFB,long PRECLLL)
+/* remark: F->KCZ changes if be_honest() fails */
+static int
+be_honest(FB_t *F, GEN nf, long PRECLLL)
 {
-  ulong av;
-  GEN MC = gmael(nf,5,2), M = gmael(nf,5,1), D = (GEN)nf[3];
-  long ex,i,j,J,k,iz,nbtest, lgsub = lg(subFB), ru = lg(MC);
-  GEN P,ideal,idealpro, exu = cgetg(ru, t_VECSMALL), MCtw= cgetg(ru, t_MAT);
+  long ex, i, j, J, k, iz, nbtest, ru, lgsub = lg(F->subFB), KCZ0 = F->KCZ;
+  GEN G, M, P, ideal, m, vdir;
+  gpmem_t av;
 
+  if (F->KCZ2 <= F->KCZ) return 1;
   if (DEBUGLEVEL)
   {
-    fprintferr("Be honest for primes from %ld to %ld\n", FB[KCZ+1],FB[KCZ2]);
+    fprintferr("Be honest for primes from %ld to %ld\n",
+               F->FB[ F->KCZ+1 ], F->FB[ F->KCZ2 ]);
     flusherr();
   }
+  if (!F->powsubFB) powsubFBgen(F, nf, CBUCHG+1, 0);
+  M = gprec_w(gmael(nf,5,1), PRECLLL);
+  G = gprec_w(gmael(nf,5,2), PRECLLL);
+  ru = lg(nf[6]);
+  vdir = cgetg(ru, t_VECSMALL);
   av = avma;
-  for (iz=KCZ+1; iz<=KCZ2; iz++, avma = av)
+  for (iz=F->KCZ+1; iz<=F->KCZ2; iz++, avma = av)
   {
-    if (DEBUGLEVEL>1) fprintferr("%ld ", FB[iz]);
-    P = idealbase[numFB[FB[iz]]]; J = lg(P);
-    /* if unramified, check all but 1 */
-    if (J > 1 && !divise(D, gmael(P,1,1))) J--;
+    long p = F->FB[iz];
+    if (DEBUGLEVEL>1) fprintferr("%ld ", p);
+    P = F->LV[p]; J = lg(P);
+    /* all P|p in FB + last is unramified --> check all but last */
+    if (isclone(P) && itou(gmael(P,J-1,3)) == 1UL /* e = 1 */) J--;
+
     for (j=1; j<J; j++)
     {
       GEN ideal0 = prime_to_ideal(nf,(GEN)P[j]);
-      ulong av2 = avma;
+      gpmem_t av2 = avma;
       for(nbtest=0;;)
       {
 	ideal = ideal0;
 	for (i=1; i<lgsub; i++)
 	{
 	  ex = mymyrand()>>randshift;
-	  if (ex) ideal = idealmulh(nf,ideal,gmael3(powsubFB,i,ex,1));
+	  if (ex) ideal = idealmulh(nf,ideal,gmael3(F->powsubFB,i,ex,1));
 	}
         ideal = remove_content(ideal);
+        for (i=1; i<ru; i++) vdir[i] = mymyrand()>>randshift;
 	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];
-          idealpro = ideallllredpart1(ideal, mulmat_real(MCtw,M), PRECLLL);
-          if (idealpro && factorgen(nf,idealpro,iz-1,FB[iz-1])) break;
-	  nbtest++; if (nbtest==200) return 0;
+          m = pseudomin(ideal, computeGtwist(nf, vdir));
+          if (m && factorgen(F,nf,ideal,m)) break;
+	  for (i=1; i<ru; i++) vdir[i] = 0;
+          vdir[k] = 10;
 	}
 	avma = av2; if (k < ru) break;
+        if (++nbtest > 50)
+        {
+          err(warner,"be_honest() failure on prime %Z\n", P[j]);
+          return 0;
+        }
       }
+      F->KCZ++; /* SUCCESS, "enlarge" factorbase */
     }
   }
   if (DEBUGLEVEL)
@@ -1841,6 +1957,7 @@ be_honest(GEN nf,GEN subFB,long PRECLLL)
     if (DEBUGLEVEL>1) fprintferr("\n");
     msgtimer("be honest");
   }
+  F->KCZ = KCZ0;
   avma = av; return 1;
 }
 
@@ -1851,20 +1968,20 @@ trunc_error(GEN x)
                         && (expo(x)>>TWOPOTBITS_IN_LONG) + 3 > lg(x);
 }
 
-/* xarch = complex logarithmic embeddings of units (u_j) found so far */
+/* A = complex logarithmic embeddings of units (u_j) found so far */
 static GEN
-compute_multiple_of_R(GEN xarch,long RU,long N,GEN *ptsublambda)
+compute_multiple_of_R(GEN A,long RU,long N,GEN *ptlambda)
 {
-  GEN v,mdet,mdet_t,Im_mdet,kR,sublambda,lambda,xreal;
-  GEN *gptr[2];
-  long av = avma, i,j, sreg = lg(xarch)-1, R1 = 2*RU - N;
+  GEN T,v,mdet,mdet_t,Im_mdet,kR,xreal,lambda;
+  long i, zc = lg(A)-1, R1 = 2*RU - N;
+  gpmem_t av = avma;
 
-  if (DEBUGLEVEL) { fprintferr("\n#### Computing regulator\n"); flusherr(); }
-  xreal=greal(xarch); v=cgetg(RU+1,t_COL); /* xreal = (log |sigma_i(u_j)|) */
-  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 */
+  if (DEBUGLEVEL) fprintferr("\n#### Computing regulator multiple\n");
+  xreal = greal(A); /* = (log |sigma_i(u_j)|) */
+  T = cgetg(RU+1,t_COL);
+  for (i=1; i<=R1; i++) T[i] = un;
+  for (   ; i<=RU; i++) T[i] = deux;
+  mdet = concatsp(xreal,T); /* det(Span(mdet)) = N * R */
 
   i = gprecision(mdet); /* truncate to avoid "near dependant" vectors */
   mdet_t = (i <= 4)? mdet: gprec_w(mdet,i-1);
@@ -1874,55 +1991,74 @@ compute_multiple_of_R(GEN xarch,long RU,long N,GEN *pt
 
   Im_mdet = vecextract_p(mdet,v);
   /* integral multiple of R: the cols we picked form a Q-basis, they have an
-   * index in the full lattice */
+   * index in the full lattice. Last column is T */
   kR = gdivgs(det2(Im_mdet), N);
   /* 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; }
-    }
-  }
-  *ptsublambda = sublambda;
-  gptr[0]=ptsublambda; gptr[1]=&kR;
-  gerepilemany(av,gptr,2); return kR;
+  lambda = gauss(Im_mdet,xreal); /* approximate rational entries */
+  for (i=1; i<=zc; i++) setlg(lambda[i], RU);
+  gerepileall(av,2, &lambda, &kR);
+  *ptlambda = lambda; 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)
+bestappr_noer(GEN x, GEN k)
 {
-  long av = avma, av2, tetpil;
-  GEN p1,c,den, R = *reg; /* multiple of regulator */
+  GEN y;
+  CATCH(precer) { y = NULL; }
+  TRY { y = bestappr(x,k); } ENDCATCH;
+  return y;
+}
 
+/* Input:
+ * lambda = approximate rational entries: coords of units found so far on a
+ * sublattice of maximal rank (sublambda)
+ * *ptkR = regulator of sublambda = multiple of regulator of lambda
+ * Compute R = true regulator of lambda.
+ *
+ * If c := Rz ~ 1, by Dirichlet's formula, then lambda is the full group of
+ * units AND the full set of relations for the class group has been computed.
+ *
+ * In fact z is a very rough approximation and we only expect 0.75 < Rz < 1.5
+ *
+ * Output: *ptkR = R, *ptU = basis of fundamental units (in terms lambda) */
+static int
+compute_R(GEN lambda, GEN z, GEN *ptL, GEN *ptkR)
+{
+  gpmem_t av = avma;
+  long r;
+  GEN L,H,D,den,R;
+  double c;
+
   if (DEBUGLEVEL) { fprintferr("\n#### Computing check\n"); flusherr(); }
-  c = gmul(R,z);
-  sublambda = bestappr(sublambda,c); den = denom(sublambda);
-  if (gcmp(den,c) > 0)
+  D = gmul2n(gmul(*ptkR,z), 1); /* bound for denom(lambda) */
+  lambda = bestappr_noer(lambda,D);
+  if (!lambda)
   {
-    if (DEBUGLEVEL) fprintferr("c = %Z\nden = %Z\n",c,den);
-    avma=av; return NULL;
+    if (DEBUGLEVEL) fprintferr("truncation error in bestappr\n");
+    return PRECI;
   }
+  den = denom(lambda);
+  if (gcmp(den,D) > 0)
+  {
+    if (DEBUGLEVEL) fprintferr("D = %Z\nden = %Z\n",D,den);
+    return PRECI;
+  }
+  L = Q_muli_to_int(lambda, den);
+  H = hnfall_i(L, NULL, 1); r = lg(H)-1;
 
-  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);
+  /* tentative regulator */
+  R = mpabs( gmul(*ptkR, gdiv(dethnf_i(H), gpowgs(den, r))) );
+  c = gtodouble(gmul(R,z)); /* should be n (= 1 if we are done) */
+  if (DEBUGLEVEL)
+  {
+    msgtimer("bestappr/regulator");
+    fprintferr("\n ***** check = %f\n",c);
+  }
+  if (c < 0.75 || c > 1.5) { avma = av; return RELAT; }
+  *ptkR = R; *ptL = L; return LARGE;
 }
 
 /* find the smallest (wrt norm) among I, I^-1 and red(I^-1) */
@@ -1930,18 +2066,15 @@ static GEN
 inverse_if_smaller(GEN nf, GEN I, long prec)
 {
   GEN J, d, dmin, I1;
-  int inv;
+
   J = (GEN)I[1];
   dmin = dethnf_i(J); I1 = idealinv(nf,I);
   J = (GEN)I1[1]; J = gmul(J,denom(J)); I1[1] = (long)J;
-  d = dethnf_i(J);
-  if (cmpii(d,dmin) < 0) {inv=1; I=I1; dmin=d;}
-  else                    inv=0;
+  d = dethnf_i(J); if (cmpii(d,dmin) < 0) {I=I1; dmin=d;}
   /* try reducing (often _increases_ the norm) */
   I1 = ideallllred(nf,I1,NULL,prec);
   J = (GEN)I1[1];
-  d = dethnf_i(J);
-  if (cmpii(d,dmin) < 0) {inv=1; I=I1;}
+  d = dethnf_i(J); if (cmpii(d,dmin) < 0) I=I1;
   return I;
 }
 
@@ -1962,17 +2095,16 @@ setlg_col(GEN U, long l)
 
 /* compute class group (clg1) + data for isprincipal (clg2) */
 static void
-class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptclg1,GEN *ptclg2,long prec)
+class_group_gen(GEN nf,GEN W,GEN C,GEN Vbase,long prec, GEN nf0,
+                GEN *ptclg1,GEN *ptclg2)
 {
-  GEN z,G,Ga,ga,GD,cyc,X,Y,D,U,V,Ur,Ui,Uir,I,J,Vbase;
-  long i,j,s,lo,lo0;
+  GEN z,G,Ga,ga,GD,cyc,X,Y,D,U,V,Ur,Ui,Uir,I,J;
+  long i,j,lo,lo0;
 
   if (DEBUGLEVEL)
-    { fprintferr("\n#### Computing class group generators\n"); timer2(); }
-  z = smith2(W); /* U W V = D, D diagonal, G = g Ui (G=new gens, g=old)  */
-  U = (GEN)z[1]; Ui = ginv(U);
-  V = (GEN)z[2];
-  D = (GEN)z[3];
+    { fprintferr("\n#### Computing class group generators\n"); (void)timer2(); }
+  D = smithall(W,&U,&V); /* UWV = D, D diagonal, G = g Ui (G=new gens, g=old) */
+  Ui = ginv(U);
   lo0 = lo = lg(D);
  /* we could set lo = lg(cyc) and truncate all matrices below
   *   setlg_col(D && U && Y, lo) + setlg(D && V && X && Ui, lo)
@@ -1987,39 +2119,32 @@ class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl
   * P invertible diagonal matrix (\pm 1) which is only implicitly defined
   * G = g Uir P + [Ga],  Uir = Ui + WX
   * g = G P Ur  + [ga],  Ur  = U + DY */
-  Vbase = cgetg(lo0,t_VEC);
-  if (typ(vperm) == t_VECSMALL)
-    for (i=1; i<lo0; i++) Vbase[i] = vectbase[vperm[i]];
-  else
-    for (i=1; i<lo0; i++) Vbase[i] = vectbase[itos((GEN)vperm[i])];
-
   G = cgetg(lo,t_VEC);
   Ga= cgetg(lo,t_VEC);
-  z = init_idele(nf);
+  z = init_famat(NULL);
+  if (!nf0) nf0 = nf;
   for (j=1; j<lo; j++)
   {
     GEN p1 = gcoeff(Uir,1,j);
-    z[1]=Vbase[1]; I = idealpowred(nf,z,p1,prec);
-    if (signe(p1)<0) I[1] = lmul((GEN)I[1],denom((GEN)I[1]));
+    z[1]=Vbase[1]; I = idealpowred(nf0,z,p1,prec);
     for (i=2; i<lo0; i++)
     {
-      p1 = gcoeff(Uir,i,j); s = signe(p1);
-      if (s)
+      p1 = gcoeff(Uir,i,j);
+      if (signe(p1))
       {
-	z[1]=Vbase[i]; J = idealpowred(nf,z,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);
+	z[1]=Vbase[i]; J = idealpowred(nf0,z,p1,prec);
+	I = idealmulh(nf0,I,J);
+	I = ideallllred(nf0,I,NULL,prec);
       }
     }
-    J = inverse_if_smaller(nf, I, prec);
+    J = inverse_if_smaller(nf0, I, prec);
     if (J != I)
     { /* update wrt P */
       neg_row(Y ,j); V[j] = lneg((GEN)V[j]);
       neg_row(Ur,j); X[j] = lneg((GEN)X[j]);
     }
     G[j] = (long)J[1]; /* generator, order cyc[j] */
-    Ga[j]= (long)J[2];
+    Ga[j]= lneg(famat_to_arch(nf, (GEN)J[2], prec));
   }
   /* at this point Y = PY, Ur = PUr, V = VP, X = XP */
 
@@ -2056,36 +2181,50 @@ class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl
   if (DEBUGLEVEL) msgtimer("classgroup generators");
 }
 
-/* real(MC * M), where columns a and b of MC have been multiplied by 2^20+1 */
-static GEN
-shift_t2(GEN T2, GEN M, GEN MC, long a, long b)
+static void
+shift_embed(GEN G, GEN Gtw, long a, long r1, long r2)
 {
-  long i,j,N = lg(T2);
-  GEN z, t2 = cgetg(N,t_MAT);
-  for (j=1; j<N; j++)
+  long j, k, l = lg(G);
+  if (a <= r1)
+    for (j=1; j<l; j++) coeff(G,a,j) = coeff(Gtw,a,j);
+  else
   {
-    t2[j] = lgetg(N,t_COL);
-    for (i=1; i<=j; i++)
+    k = (a<<1) - r1;
+    for (j=1; j<l; j++)
     {
-      z = mul_real(gcoeff(MC,i,a), gcoeff(M,a,j));
-      if (a!=b) z = gadd(z, mul_real(gcoeff(MC,i,b), gcoeff(M,b,j)));
-      coeff(t2,j,i) = coeff(t2,i,j) = ladd(gcoeff(T2,i,j), gmul2n(z,20));
+      coeff(G,k-1,j) = coeff(Gtw,k-1,j);
+      coeff(G,k  ,j) = coeff(Gtw,k,  j);
     }
   }
-  return t2;
 }
 
+/* G where embeddings a and b are multiplied by 2^10 */
 static GEN
-compute_vecT2(long RU,GEN nf)
+shift_G(GEN G, GEN Gtw, long a, long b, long r1, long r2)
 {
-  GEN vecT2, M = gmael(nf,5,1), MC = gmael(nf,5,2), T2 = gmael(nf,5,3);
-  long i,j,n = min(RU,9), ind = 1;
+  GEN g = dummycopy(G);
+  shift_embed(g,Gtw,a,r1,r2);
+  shift_embed(g,Gtw,b,r1,r2); return g;
+}
 
-  vecT2=cgetg(1 + n*(n+1)/2,t_VEC);
-  for (j=1; j<=n; j++)
-    for (i=1; i<=j; i++) vecT2[ind++] = (long)shift_t2(T2,M,MC,i,j);
-  if (DEBUGLEVEL) msgtimer("weighted T2 matrices");
-  return vecT2;
+static GEN
+compute_vecG(GEN nf,long prec)
+{
+  GEN vecG, Gtw, M = gmael(nf,5,1), G = gmael(nf,5,2);
+  long r1,r2,i,j,ind, n = min(lg(M[1])-1, 9);
+
+  nf_get_sign(nf,&r1,&r2);
+  vecG=cgetg(1 + n*(n+1)/2,t_VEC);
+  if (nfgetprec(nf) > prec)
+  {
+    M = gprec_w(M,prec);
+    G = gprec_w(G,prec);
+  }
+  Gtw = gmul2n(G, 10);
+  for (ind=j=1; j<=n; j++)
+    for (i=1; i<=j; i++) vecG[ind++] = (long)shift_G(G,Gtw,i,j,r1,r2);
+  if (DEBUGLEVEL) msgtimer("weighted G matrices");
+  return vecG;
 }
 
 /* cf. relationrank()
@@ -2129,12 +2268,13 @@ addcolumntomatrix(GEN V, GEN invp, GEN L)
  * Starting from the identity (canonical basis of Q^n), we obtain a base
  * change matrix P by taking the independant columns of A and vectors from
  * the canonical basis. Update invp = 1/P, and L in {0,1}^n, with L[i] = 0
- * of P[i] = e_i, and 1 if P[i] = A_i (i-th column of A)
+ * if P[i] = e_i, and 1 if P[i] = A_i (i-th column of A)
  */
 static GEN
 relationrank(GEN *A, long r, GEN L)
 {
-  long i, n = lg(L)-1, av = avma, lim = stack_lim(av,1);
+  long i, n = lg(L)-1;
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   GEN invp = idmat(n);
 
   if (!r) return invp;
@@ -2152,32 +2292,51 @@ relationrank(GEN *A, long r, GEN L)
   return gerepilecopy(av, invp);
 }
 
+/* SMALLBUCHINIT */
+
 static GEN
-codeprime(GEN bnf, GEN pr)
+decode_pr_lists(GEN nf, GEN pfc)
 {
-  long j,av=avma,tetpil;
-  GEN p,al,fa,p1;
+  long i, p, pmax, n = degpol(nf[1]), l = lg(pfc);
+  GEN t, L;
 
-  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 */
+  pmax = 0;
+  for (i=1; i<l; i++)
+  {
+    t = (GEN)pfc[i]; p = itos(t) / n;
+    if (p > pmax) pmax = p;
+  }
+  L = cgetg(pmax+1, t_VEC);
+  for (p=1; p<=pmax; p++) L[p] = 0;
+  for (i=1; i<l; i++)
+  {
+    t = (GEN)pfc[i]; p = itos(t) / n;
+    if (!L[p]) L[p] = (long)primedec(nf, stoi(p));
+  }
+  return L;
 }
 
 static GEN
-decodeprime(GEN nf, GEN co)
+decodeprime(GEN T, GEN L, long n)
 {
-  long n,indi,av=avma;
-  GEN p,rem,p1;
+  long t = itos(T);
+  return gmael(L, t/n, t%n + 1);
+}
+static GEN
+codeprime(GEN L, long N, GEN pr)
+{
+  long p = itos((GEN)pr[1]);
+  return utoi( N*p + pr_index((GEN)L[p], pr)-1 );
+}
 
-  n=lg(nf[7])-1; p=dvmdis(co,n,&rem); indi=itos(rem)+1;
-  p1=primedec(nf,p);
-  return gerepilecopy(av,(GEN)p1[indi]);
+static GEN
+codeprimes(GEN Vbase, long N)
+{
+  GEN v, L = get_pr_lists(Vbase, N, 1);
+  long i, l = lg(Vbase);
+  v = cgetg(l, t_VEC);
+  for (i=1; i<l; i++) v[i] = (long)codeprime(L, N, (GEN)Vbase[i]);
+  return v;
 }
 
 /* v = bnf[10] */
@@ -2226,15 +2385,8 @@ makecycgen(GEN bnf)
       if (y && !fact_ok(nf,y,NULL,gen,(GEN)D[i])) y = NULL;
       if (y) { h[i] = (long)to_famat_all(y,gun); continue; }
     }
-    y = isprincipalfact(bnf, gen, (GEN)D[i], NULL,
-                        nf_GEN|nf_GENMAT|nf_FORCE|nf_GIVEPREC);
-    if (typ(y) != t_INT)
-      h[i] = y[2];
-    else
-    {
-      y = idealpow(nf, (GEN)gen[i], (GEN)cyc[i]);
-      h[i] = isprincipalgenforce(bnf,y)[2];
-    }
+    y = isprincipalfact(bnf, gen, (GEN)D[i], NULL, nf_GENMAT|nf_FORCE);
+    h[i] = y[2];
   }
   return h;
 }
@@ -2243,7 +2395,7 @@ makecycgen(GEN bnf)
 static GEN
 makematal(GEN bnf)
 {
-  GEN W,B,pFB,vp,nf,ma, WB_C;
+  GEN W,B,pFB,nf,ma, WB_C;
   long lW,lma,j,prec;
 
   ma = get_matal((GEN)bnf[10]);
@@ -2252,11 +2404,9 @@ makematal(GEN bnf)
   W   = (GEN)bnf[1];
   B   = (GEN)bnf[2];
   WB_C= (GEN)bnf[4];
-  vp  = (GEN)bnf[6];
   nf  = (GEN)bnf[7];
   lW=lg(W)-1; lma=lW+lg(B);
-  pFB=cgetg(lma,t_VEC); ma=(GEN)bnf[5]; /* reorder factor base */
-  for (j=1; j<lma; j++) pFB[j] = ma[itos((GEN)vp[j])];
+  pFB = get_Vbase(bnf);
   ma = cgetg(lma,t_MAT);
 
   prec = prec_arch(bnf);
@@ -2274,7 +2424,7 @@ makematal(GEN bnf)
       ma[j] = (long)y; continue;
     }
 
-    if (!y) y = isprincipalfact(bnf,pFB,ex,C, nf_GEN|nf_FORCE|nf_GIVEPREC);
+    if (!y) y = isprincipalfact(bnf,pFB,ex,C, nf_GENMAT|nf_FORCE|nf_GIVEPREC);
     if (typ(y) != t_INT)
     {
       if (DEBUGLEVEL>1) fprintferr("%ld ",j);
@@ -2284,7 +2434,7 @@ makematal(GEN bnf)
     prec = itos(y); j--;
     if (DEBUGLEVEL) err(warnprec,"makematal",prec);
     nf = nfnewprec(nf,prec);
-    bnf = bnfinit0(nf,1,NULL,prec); setrand(c);
+    bnf = bnfinit0(nf,1,NULL,prec); (void)setrand(c);
   }
   if (DEBUGLEVEL>1) fprintferr("\n");
   return ma;
@@ -2316,7 +2466,7 @@ check_and_build_cycgen(GEN bnf)
   GEN cycgen = get_cycgen((GEN)bnf[10]);
   if (!cycgen)
   {
-    long av = avma;
+    gpmem_t av = avma;
     if (DEBUGLEVEL) err(warner,"completing bnf (building cycgen)");
     bnfinsert(bnf, makecycgen(bnf), 2); avma = av;
     cycgen = get_cycgen((GEN)bnf[10]);
@@ -2330,168 +2480,120 @@ check_and_build_matal(GEN bnf)
   GEN matal = get_matal((GEN)bnf[10]);
   if (!matal)
   {
-    long av = avma;
+    gpmem_t av = avma;
     if (DEBUGLEVEL) err(warner,"completing bnf (building matal)");
     bnfinsert(bnf, makematal(bnf), 1); avma = av;
     matal = get_matal((GEN)bnf[10]);
   }
   return matal;
 }
-			
+
 GEN
 smallbuchinit(GEN pol,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,long minsFB,long prec)
 {
-  long av=avma,av1,k;
-  GEN y,bnf,pFB,vp,nf,mas,res,uni,v1,v2,v3;
+  GEN y, bnf, nf, res, p1;
+  gpmem_t av = avma;
 
   if (typ(pol)==t_VEC) bnf = checkbnf(pol);
   else
   {
-    bnf=buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsFB,-3,prec);
+    const long fl = nf_INIT | nf_UNITS | nf_FORCE;
+    bnf = buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsFB,fl,prec);
     bnf = checkbnf_discard(bnf);
   }
-  pFB=(GEN)bnf[5]; vp=(GEN)bnf[6]; nf=(GEN)bnf[7];
-  mas=(GEN)nf[5]; res=(GEN)bnf[8]; uni=(GEN)res[5];
+  nf  = (GEN)bnf[7];
+  res = (GEN)bnf[8];
 
-  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]);
-  av1 = avma;
-  y[12] = gcmp0((GEN)bnf[10])? lpilecopy(av1, makematal(bnf))
-                             : lcopy((GEN)bnf[10]);
-  return gerepileupto(av,y);
+  y = cgetg(13,t_VEC);
+  y[1] = nf[1];
+  y[2] = mael(nf,2,1);
+  y[3] = nf[3];
+  y[4] = nf[7];
+  y[5] = nf[6];
+  y[6] = mael(nf,5,5);
+  y[7] = bnf[1];
+  y[8] = bnf[2];
+  y[9] = (long)codeprimes((GEN)bnf[5], degpol(nf[1]));
+
+  p1 = cgetg(3, t_VEC);
+  p1[1] = mael(res,4,1);
+  p1[2] = (long)algtobasis(bnf,gmael(res,4,2));
+  y[10] = (long)p1;
+
+  y[11] = (long)algtobasis(bnf, (GEN)res[5]);
+  y[12] = gcmp0((GEN)bnf[10])? (long)makematal(bnf): bnf[10];
+  return gerepilecopy(av, y);
 }
 
 static GEN
-get_regulator(GEN mun,long prec)
+get_regulator(GEN mun)
 {
-  long av,tetpil;
-  GEN p1;
+  gpmem_t av = avma;
+  GEN A;
 
   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));
+  A = gtrans( greal(mun) );
+  setlg(A, lg(A)-1);
+  return gerepileupto(av, gabs(det(A), 0));
 }
 
+/* return corrected archimedian components for elts of x (vector)
+ * (= log(sigma_i(x)) - log(|Nx|) / [K:Q]) */
 static GEN
-log_poleval(GEN P, GEN *ro, long i, GEN nf, long prec0)
+get_archclean(GEN nf, GEN x, long prec, int units)
 {
-  GEN x = poleval(P, (GEN)(*ro)[i]);
-  long prec = prec0, k = 0;
-  while (gcmp0(x) || ((k = gprecision(x)) && k < DEFAULTPREC))
-  {
-    prec = (prec-2)<<1;
-    if (DEBUGLEVEL) err(warnprec,"log_poleval",prec);
-    *ro = get_roots((GEN)nf[1], itos(gmael(nf,2,1)),lg(*ro)-1,prec);
-    x = poleval(P, (GEN)(*ro)[i]);
-  }
-  if (k > prec0) x = gprec_w(x,prec0);
-  return glog(x, prec0);
-}
+  long k,N, la = lg(x);
+  GEN M = cgetg(la,t_MAT);
 
-/* return corrected archimedian components
- * (= log(sigma_i(x)) - log(|Nx|)/[K:Q]) for a (matrix of elements) */
-static GEN
-get_arch2_i(GEN nf, GEN a, long prec, int units)
-{
-  GEN M,N, ro = dummycopy((GEN)nf[6]);
-  long j,k, la = lg(a), ru = lg(ro), r1 = nf_get_r1(nf);
-
-  M = cgetg(la,t_MAT); if (la == 1) return M;
-  if (typ(a[1]) == t_COL) a = gmul((GEN)nf[7], a);
-  if (units) N = NULL; /* no correction necessary */
-  else
-  {
-    GEN pol = (GEN)nf[1];
-    N = cgetg(la,t_VEC);
-    for (k=1; k<la; k++) N[k] = (long)gabs(subres(pol,(GEN)a[k]),0);
-    N = gdivgs(glog(N,prec), - (degpol(pol))); /* - log(|norm|) / [K:Q] */
-  }
+  if (la == 1) return M;
+  N = degpol(nf[1]);
   for (k=1; k<la; k++)
   {
-    GEN p, c = cgetg(ru,t_COL); M[k] = (long)c;
-    for (j=1; j<ru; j++)
-    {
-      p = log_poleval((GEN)a[k],&ro,j,nf,prec);
-      if (N) p = gadd(p, (GEN)N[k]);
-      c[j]=(j<=r1)? (long) p: lmul2n(p,1);
-    }
+    GEN c = get_arch(nf, (GEN)x[k], prec);
+    if (!units) c = cleanarch(c, N, prec);
+    M[k] = (long)c;
   }
   return M;
 }
 
-static GEN
-get_arch2(GEN nf, GEN a, long prec, int units)
-{
-  long av = avma;
-  return gerepilecopy(av, get_arch2_i(nf,a,prec,units));
-}
-
 static void
-my_class_group_gen(GEN bnf, GEN *ptcl, GEN *ptcl2, long prec)
+my_class_group_gen(GEN bnf, long prec, GEN nf0, GEN *ptcl, GEN *ptcl2)
 {
-  GEN W=(GEN)bnf[1], C=(GEN)bnf[4], vperm=(GEN)bnf[6], nf=(GEN)bnf[7], *gptr[2];
-  long av = avma;
-
-  vectbase = (GEN)bnf[5]; /* needed by class_group_gen */
-  class_group_gen(nf,W,C,vperm,ptcl,ptcl2, prec);
-  gptr[0]=ptcl; gptr[1]=ptcl2; gerepilemany(av,gptr,2);
+  GEN W=(GEN)bnf[1], C=(GEN)bnf[4], nf=(GEN)bnf[7];
+  class_group_gen(nf,W,C,get_Vbase(bnf),prec,nf0, ptcl,ptcl2);
 }
 
 GEN
 bnfnewprec(GEN bnf, long prec)
 {
-  GEN nf,ro,res,p1,funits,mun,matal,clgp,clgp2,y;
-  long r1,r2,ru,pl1,pl2,prec1;
+  GEN nf0 = (GEN)bnf[7], nf, res, funits, mun, matal, clgp, clgp2, y;
+  gpmem_t av = avma;
+  long r1, r2, prec1;
 
   bnf = checkbnf(bnf);
   if (prec <= 0) return nfnewprec(checknf(bnf),prec);
-  y = cgetg(11,t_VEC);
-  funits = check_units(bnf,"bnfnewprec");
   nf = (GEN)bnf[7];
-  ro = (GEN)nf[6]; ru = lg(ro)-1;
-  r1 = nf_get_r1(nf);
-  r2 = (r1 + ru) >> 1;
-  pl1 = (ru == 1 && r1 == 0)? 0: gexpo(funits);
-  pl2 = gexpo(ro);
+  nf_get_sign(nf, &r1, &r2);
+  funits = algtobasis(nf, check_units(bnf,"bnfnewprec"));
+
   prec1 = prec;
-  prec += ((ru + r2 - 1) * (pl1 + (ru + r2) * pl2)) >> TWOPOTBITS_IN_LONG;
-  nf = nfnewprec((GEN)bnf[7],prec);
-  res = cgetg(7,t_VEC);
-  ro = (GEN)nf[6];
-  mun = get_arch2(nf,funits,prec,1);
+  if (r2 > 1 || r1 != 0)
+    prec += 1 + (gexpo(funits) >> TWOPOTBITS_IN_LONG);
+  nf = nfnewprec(nf0,prec);
+  mun = get_archclean(nf,funits,prec,1);
   if (prec != prec1) { mun = gprec_w(mun,prec1); prec = prec1; }
-  res[2]=(long)get_regulator(mun,prec);
-  p1 = (GEN)bnf[8];
-  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;
   matal = check_and_build_matal(bnf);
-  y[4]=(long)get_arch2(nf,matal,prec,0);
-  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,prec);
-  res[1]=(long)clgp;
-  y[9]=(long)clgp2;
-  y[10]=lcopy((GEN)bnf[10]); return y;
+
+  y = dummycopy(bnf);
+  y[3] = (long)mun;
+  y[4] = (long)get_archclean(nf,matal,prec,0);
+  y[7] = (long)nf;
+  my_class_group_gen(y,prec,nf0, &clgp,&clgp2);
+  res = dummycopy((GEN)bnf[8]);
+  res[1] = (long)clgp;
+  res[2] = (long)get_regulator(mun);
+  y[8] = (long)res;
+  y[9] = (long)clgp2; return gerepilecopy(av, y);
 }
 
 GEN
@@ -2505,72 +2607,82 @@ bnrnewprec(GEN bnr, long prec)
   return y;
 }
 
+static void
+nfbasic_from_sbnf(GEN sbnf, nfbasic_t *T)
+{
+  T->x    = (GEN)sbnf[1];
+  T->dK   = (GEN)sbnf[3];
+  T->bas  = (GEN)sbnf[4];
+  T->index= get_nfindex(T->bas);
+  T->r1   = itos((GEN)sbnf[2]);
+  T->dx   = NULL;
+  T->lead = NULL;
+  T->basden = NULL;
+}
+
+static GEN
+get_clfu(GEN clgp, GEN reg, GEN c1, GEN zu, GEN fu, long k)
+{
+  GEN z = cgetg(7, t_VEC);
+  z[1]=(long)clgp; z[2]=(long)reg; z[3]=(long)c1;
+  z[4]=(long)zu;   z[5]=(long)fu;  z[6]=lstoi(k); return z;
+}
+
 GEN
 bnfmake(GEN sbnf, long prec)
 {
-  long av = avma, j,k,n,r1,r2,ru,lpf;
-  GEN p1,x,bas,ro,nf,mun,funits,index;
-  GEN pfc,vp,mc,clgp,clgp2,res,y,W,racu,reg,matal;
+  long j, k, l, n;
+  gpmem_t av = avma;
+  GEN p1, bas, ro, nf, mun, fu, L;
+  GEN pfc, mc, clgp, clgp2, res, y, W, zu, reg, matal, Vbase;
+  nfbasic_t T;
 
-  if (typ(sbnf)!=t_VEC || lg(sbnf)!=13)
-    err(talker,"incorrect sbnf in bnfmake");
-  x=(GEN)sbnf[1]; bas=(GEN)sbnf[4]; n=lg(bas)-1;
-  r1=itos((GEN)sbnf[2]); r2=(n-r1)>>1; ru=r1+r2;
-  ro=(GEN)sbnf[5];
-  if (prec > gprecision(ro)) ro=get_roots(x,r1,ru,prec);
-  index = gun;
-  for (j=2; j<=n; j++) index = mulii(index, denom(leading_term(bas[j])));
+  if (typ(sbnf) != t_VEC || lg(sbnf) != 13) err(typeer,"bnfmake");
+  if (prec < DEFAULTPREC) prec = DEFAULTPREC;
 
-  nf=cgetg(10,t_VEC);
-  nf[1]=sbnf[1]; p1=cgetg(3,t_VEC); p1[1]=lstoi(r1); p1[2]=lstoi(r2);
-  nf[2]=(long)p1;
-  nf[3]=sbnf[3];
-  nf[4]=(long)index;
-  nf[6]=(long)ro;
-  nf[7]=(long)bas;
-  get_nf_matrices(nf, 0);
+  nfbasic_from_sbnf(sbnf, &T);
+  ro = (GEN)sbnf[5];
+  if (prec > gprecision(ro)) ro = get_roots(T.x,T.r1,prec);
+  nf = nfbasic_to_nf(&T, ro, prec);
+  bas = (GEN)nf[7];
 
-  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_arch2_i(nf,funits,prec,1);
+  p1 = (GEN)sbnf[11]; l = lg(p1); fu = cgetg(l, t_VEC);
+  for (k=1; k < l; k++) fu[k] = lmul(bas, (GEN)p1[k]);
+  mun = get_archclean(nf,fu,prec,1);
 
-  prec=gprecision(ro); if (prec<DEFAULTPREC) prec=DEFAULTPREC;
+  prec = gprecision(ro);
   matal = get_matal((GEN)sbnf[12]);
   if (!matal) matal = (GEN)sbnf[12];
-  mc = get_arch2_i(nf,matal,prec,0);
+  mc = get_archclean(nf,matal,prec,0);
 
-  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]);
-  }
+  pfc = (GEN)sbnf[9];
+  l = lg(pfc);
+  Vbase = cgetg(l,t_COL);
+  L = decode_pr_lists(nf, pfc);
+  n = degpol(nf[1]);
+  for (j=1; j<l; j++) Vbase[j] = (long)decodeprime((GEN)pfc[j], L, n);
   W = (GEN)sbnf[7];
-  class_group_gen(nf,W,mc,vp,&clgp,&clgp2, prec); /* uses vectbase */
+  class_group_gen(nf,W,mc,Vbase,prec,NULL, &clgp,&clgp2);
 
-  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;
+  reg = get_regulator(mun);
+  p1 = cgetg(3,t_VEC); zu = (GEN)sbnf[10];
+  p1[1] = zu[1];
+  p1[2] = lmul(bas,(GEN)zu[2]); zu = 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);
-
+  res = get_clfu(clgp,reg,dbltor(1.0),zu,fu,1000);
   y=cgetg(11,t_VEC);
-  y[1]=(long)W;   y[2]=sbnf[8];        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]=sbnf[12];
-  return gerepilecopy(av,y);
+  y[1]=(long)W;   y[2]=sbnf[8];     y[3]=(long)mun;
+  y[4]=(long)mc;  y[5]=(long)Vbase; y[6]=zero;
+  y[7]=(long)nf;  y[8]=(long)res;   y[9]=(long)clgp2;
+  y[10] = sbnf[12]; return gerepilecopy(av,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;
+  long fl, lx, minsFB=3, nbrelpid=4;
+  gpmem_t av=avma;
   GEN bach=doubl,bach2=doubl,RELSUP=court,borne=gun;
 
   if (!data) lx=1;
@@ -2594,17 +2706,17 @@ classgroupall(GEN P, GEN data, long flag, long prec)
   }
   switch(flag)
   {
-    case 0: flun=-2; break;
-    case 1: flun=-3; break;
-    case 2: flun=-1; break;
+    case 0: fl = nf_INIT | nf_UNITS; break;
+    case 1: fl = nf_INIT | nf_UNITS | nf_FORCE; break;
+    case 2: fl = nf_INIT | nf_ROOT1; 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;
+    case 4: fl = nf_UNITS; break;
+    case 5: fl = nf_UNITS | nf_FORCE; break;
+    case 6: fl = 0; break;
     default: err(flagerr,"classgroupall");
       return NULL; /* not reached */
   }
-  return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,flun,prec);
+  return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,fl,prec);
 }
 
 GEN
@@ -2633,7 +2745,7 @@ bnfinit0(GEN P, long flag, GEN data, long prec)
 GEN
 classgrouponly(GEN P, GEN data, long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN z;
 
   if (typ(P)==t_INT)
@@ -2648,7 +2760,7 @@ classgrouponly(GEN P, GEN data, long prec)
 GEN
 regulator(GEN P, GEN data, long prec)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN z;
 
   if (typ(P)==t_INT)
@@ -2673,7 +2785,7 @@ col_0(long n)
    GEN c = (GEN) gpmalloc(sizeof(long)*(n+1));
    long i;
    for (i=1; i<=n; i++) c[i]=0;
-   c[0] = evaltyp(t_VECSMALL) | evallg(n);
+   c[0] = evaltyp(t_VECSMALL) | evallg(n+1);
    return c;
 }
 
@@ -2687,129 +2799,163 @@ cgetc_col(long n, long prec)
 }
 
 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)
+buchall_end(GEN nf,GEN CHANGE,long fl,GEN res, GEN clg2, GEN W, GEN B,
+            GEN A, GEN matarch, GEN Vbase)
 {
-  long i, l = labs(fl)>1? 11: fl? 9: 8;
-  GEN p1,z, RES = cgetg(11,t_COL);
+  long l = (fl & nf_UNITS)? 7
+                          : (fl & nf_ROOT1)? 5: 4;
+  GEN p1, z;
 
-  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)
+  setlg(res, l);
+  if (! (fl & nf_INIT))
   {
-    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;
+    GEN x = cgetg(5, t_VEC);
+    x[1]=nf[1];
+    x[2]=nf[2]; p1=cgetg(3,t_VEC); p1[1]=nf[3]; p1[2]=nf[4];
+    x[3]=(long)p1;
+    x[4]=nf[7];
+    z=cgetg(2,t_MAT); z[1]=(long)concatsp(x, res); return z;
   }
   z=cgetg(11,t_VEC);
   z[1]=(long)W;
   z[2]=(long)B;
-  z[3]=(long)xarch;
+  z[3]=(long)A;
   z[4]=(long)matarch;
-  z[5]=(long)vectbase;
-  for (i=lg(vperm)-1; i>0; i--) vperm[i]=lstoi(vperm[i]);
-  settyp(vperm, t_VEC);
-  z[6]=(long)vperm;
-  z[7]=(long)nf; RES+=4; RES[0]=evaltyp(t_VEC) | evallg(l-4);
-  z[8]=(long)RES;
+  z[5]=(long)Vbase;
+  z[6]=zero;
+  z[7]=(long)nf;
+  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);
+  return 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);
+  gpmem_t av = avma;
+  GEN W,B,A,matarch,Vbase,res;
+  GEN fu=cgetg(1,t_VEC), R=gun, c1=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);
+  clg1[1]=un; clg1[2]=clg1[3]=clg2[2]=clg2[3]=lgetg(1,t_VEC);
+  clg2[1]=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);
+  W=B=A=matarch=cgetg(1,t_MAT);
+  Vbase=cgetg(1,t_COL);
 
-  return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c_1,zu,W,B,xarch,matarch,vectbase,vperm));
+  res = get_clfu(clg1, R, c1, zu, fu, EXP220);
+  return gerepilecopy(av, buchall_end(nf,CHANGE,flun,res,clg2,W,B,A,matarch,Vbase));
 }
 
-static GEN
-get_LLLnf(GEN nf, long prec)
+/* return (small set of) indices of columns generating the same lattice as x.
+ * Assume HNF(x) is inexpensive (few rows, many columns).
+ * Dichotomy approach since interesting columns may be at the very end */
+GEN
+extract_full_lattice(GEN x)
 {
-  GEN M  = gmael(nf,5,1);
-  GEN T2 = gmael(nf,5,3);
-  GEN cbase = lllgramintern(T2,100,1,prec);
-  GEN v = cgetg(5,t_VEC);
-  if (!cbase) return NULL;
-  if (gegal(cbase, idmat(lg(T2)-1))) cbase = NULL;
-  v[1] = (long) (cbase? qf_base_change(T2, cbase, 1): T2);
-  v[2] = (long) (cbase? gmul(M, cbase): M);
-  v[3] = (long) cbase;
-  v[4] = (long) (cbase? ZM_inv(cbase,gun): NULL); return v;
+  long dj, j, k, l = lg(x);
+  GEN h, h2, H, v;
+
+  if (l < 200) return NULL; /* not worth it */
+
+  v = cgetg(l, t_VECSMALL);
+  setlg(v, 1);
+  H = hnfall_i(x, NULL, 1);
+  h = cgetg(1, t_MAT);
+  dj = 1;
+  for (j = 1; j < l; )
+  {
+    gpmem_t av = avma;
+    long lv = lg(v);
+
+    for (k = 0; k < dj; k++) v[lv+k] = j+k;
+    setlg(v, lv + dj);
+    h2 = hnfall_i(vecextract_p(x, v), NULL, 1);
+    if (gegal(h, h2))
+    { /* these dj columns can be eliminated */
+      avma = av; setlg(v, lv);
+      j += dj;
+      if (j >= l) break; /* shouldn't occur */
+      dj <<= 1;
+      if (j + dj >= l) dj = (l - j) >> 1;
+    }
+    else if (dj > 1)
+    { /* at least one interesting column, try with first half of this set */
+      avma = av; setlg(v, lv);
+      dj >>= 1;
+    }
+    else
+    { /* this column should be kept */
+      if (gegal(h2, H)) break;
+      h = h2; j++;
+    }
+  }
+  return v;
 }
 
 GEN
 buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,
         long minsFB,long flun,long prec)
 {
-  ulong av = avma,av0,av1,limpile;
-  long N,R1,R2,RU,PRECREG,PRECLLL,KCCO,RELSUP,LIMC,LIMC2,lim;
-  long nlze,sreg,nrelsup,nreldep,phase,matmax,i,j,k,ss,cglob;
-  long first=1, sfb_increase=0, sfb_trials=0, precdouble=0, precadd=0;
-  double cbach,cbach2,drc,LOGD2;
-  GEN p1,vecT2,fu,zu,nf,LLLnf,D,xarch,W,reg,lfun,z,clh,vperm,subFB;
-  GEN B,C,c1,sublambda,pdep,parch,liste,invp,clg1,clg2, *mat;
-  GEN CHANGE=NULL, extramat=NULL, extraC=NULL, list_jideal=NULL;
+  gpmem_t av = avma, av0, av1, limpile;
+  long N,R1,R2,RU,PRECREG,PRECLLL,PRECLLLadd,KCCO,RELSUP,LIMC,LIMC2,lim;
+  long nlze,zc,nrelsup,nreldep,phase,matmax,i,j,k,seed,MAXRELSUP;
+  long sfb_increase, sfb_trials, precdouble=0, precadd=0;
+  long cglob; /* # of relations found so far */
+  double cbach, cbach2, drc, LOGD2;
+  GEN vecG,fu,zu,nf,D,A,W,R,Res,z,h,list_jideal;
+  GEN res,L,resc,B,C,lambda,pdep,liste,invp,clg1,clg2,Vbase;
+  GEN *mat;     /* raw relations found (as VECSMALL, not HNF-reduced) */
+  GEN first_nz; /* first_nz[i] = index of 1st non-0 entry in mat[i] */
+  GEN CHANGE=NULL, extramat=NULL, extraC=NULL;
   char *precpb = NULL;
+  FB_t F;
 
-  if (DEBUGLEVEL) timer2();
+  if (DEBUGLEVEL) (void)timer2();
 
-  if (typ(P)==t_POL) nf = NULL; else { nf = checknf(P); P = (GEN)nf[1]; }
+  P = get_nfpol(P, &nf);
   if (typ(gRELSUP) != t_INT) gRELSUP = gtrunc(gRELSUP);
   RELSUP = itos(gRELSUP);
   if (RELSUP<=0) err(talker,"not enough relations in bnfxxx");
 
   /* Initializations */
   fu = NULL; /* gcc -Wall */
-  N = degpol(P); PRECREG = max(BIGDEFAULTPREC,prec);
+  N = degpol(P);
+  PRECREG = max(BIGDEFAULTPREC,prec);
+  PRECLLLadd = DEFAULTPREC;
   if (!nf)
   {
-    nf = initalgall0(P, nf_REGULAR, PRECREG);
+    nf = initalg(P, PRECREG);
     /* P was non-monic and nfinit changed it ? */
     if (lg(nf)==3) { CHANGE = (GEN)nf[2]; nf = (GEN)nf[1]; }
   }
-  if (N<=1) return buchall_for_degree_one_pol(nf,CHANGE,flun);
+  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("initalg & rootsof1");
 
-  R1 = itos(gmael(nf,2,1)); R2 = (N-R1)>>1; RU = R1+R2;
+  nf_get_sign(nf, &R1, &R2); RU = R1+R2;
   D = (GEN)nf[3]; drc = fabs(gtodouble(D));
   LOGD2 = log(drc); LOGD2 = LOGD2*LOGD2;
   lim = (long) (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;
+  if (cbach == 0.) err(talker,"Bach constant = 0 in bnfxxx");
+  if (nbrelpid <= 0) gborne = gzero;
+
   cbach2 = gtodouble(gcbach2);
+  /* resc ~ sqrt(D) w / 2^r1 (2pi)^r2 = hR / Res(zeta_K, s=1) */
+  resc = gdiv(mulri(gsqrt(absi(D),DEFAULTPREC), (GEN)zu[1]),
+              gmul2n(gpowgs(Pi2n(1,DEFAULTPREC), R2), R1));
   if (DEBUGLEVEL)
-  {
-    fprintferr("N = %ld, R1 = %ld, R2 = %ld, RU = %ld\n",N,R1,R2,RU);
-    fprintferr("D = %Z\n",D);
-  }
-  av0 = avma; mat = NULL; FB = NULL;
+    fprintferr("N = %ld, R1 = %ld, R2 = %ld\nD = %Z\n",N,R1,R2,D);
+  av0 = avma; mat = NULL; first_nz = NULL;
 
 START:
-  avma = av0;
-  if (first) first = 0; else desallocate(&mat);
+  seed = getrand();
+  avma = av0; desallocate(&mat, &first_nz);
   if (precpb)
   {
     precdouble++;
@@ -2829,74 +2975,88 @@ START:
     cbach = check_bach(cbach,12.);
   precadd = 0;
 
-  /* T2-LLL-reduce nf.zk */
-  LLLnf = get_LLLnf(nf, PRECREG);
-  if (!LLLnf) { precpb = "LLLnf"; goto START; }
-
-  LIMC = (long)(cbach*LOGD2); if (LIMC < 20) LIMC = 20;
+  LIMC = (long)(cbach*LOGD2);
+  if (LIMC < 20) { LIMC = 20; cbach = LIMC / LOGD2; }
   LIMC2= max(3 * N, (long)(cbach2*LOGD2));
   if (LIMC2 < LIMC) LIMC2 = LIMC;
   if (DEBUGLEVEL) { fprintferr("LIMC = %ld, LIMC2 = %ld\n",LIMC,LIMC2); }
 
-  /* initialize FB, [sub]vperm */
-  lfun = FBgen(nf,LIMC2,LIMC);
-  if (!lfun) goto START;
-  vperm = cgetg(lg(vectbase), t_VECSMALL);
+  Res = FBgen(&F, nf, LIMC2, LIMC);
+  if (!Res || !subFBgen(&F, nf, min(lim,LIMC2) + 0.5, minsFB)) goto START;
+
+  if (DEBUGLEVEL)
+  {
+    if (DEBUGLEVEL>3)
+    {
+      fprintferr("\n***** IDEALS IN FACTORBASE *****\n\n");
+      for (i=1; i<=F.KC; i++) fprintferr("no %ld = %Z\n",i,F.LP[i]);
+      fprintferr("\n***** IDEALS IN SUB FACTORBASE *****\n\n");
+      outerr(vecextract_p(F.LP,F.subFB));
+      fprintferr("\n***** INITIAL PERMUTATION *****\n\n");
+      fprintferr("perm = %Z\n\n",F.perm);
+    }
+    msgtimer("sub factorbase (%ld elements)",lg(F.subFB)-1);
+  }
   sfb_trials = sfb_increase = 0;
-  subFB = subFBgen(N,min(lim,LIMC2), minsFB, vperm, &ss);
-  if (!subFB) goto START;
   nreldep = nrelsup = 0;
 
   /* PRECLLL = prec for LLL-reductions (idealred)
    * PRECREG = prec for archimedean components */
-  PRECLLL = DEFAULTPREC
-          + ((expi(D)*(lg(subFB)-2)+((N*N)>>2))>>TWOPOTBITS_IN_LONG);
+  PRECLLL = PRECLLLadd
+          + ((expi(D)*(lg(F.subFB)-2) + ((N*N)>>2)) >> TWOPOTBITS_IN_LONG);
   if (!precdouble) PRECREG = prec+1;
-  if (PRECREG < PRECLLL) PRECREG = PRECLLL;
-  KCCO = KC+RU-1 + max(ss,RELSUP); /* expected number of needed relations */
+  if (PRECREG < PRECLLL)
+  { /* very rare */
+    PRECREG = PRECLLL;
+    nf = nfnewprec(nf,PRECREG); av0 = avma;
+  }
+  KCCO = F.KC+RU-1 + RELSUP; /* expected # of needed relations */
   if (DEBUGLEVEL)
-    fprintferr("relsup = %ld, ss = %ld, KCZ = %ld, KC = %ld, KCCO = %ld\n",
-               RELSUP,ss,KCZ,KC,KCCO);
+    fprintferr("relsup = %ld, KCZ = %ld, KC = %ld, KCCO = %ld\n",
+               RELSUP, F.KCZ, F.KC, KCCO);
+  MAXRELSUP = min(50, 4*F.KC);
   matmax = 10*KCCO + MAXRELSUP; /* make room for lots of relations */
-  mat = (GEN*)gpmalloc(sizeof(GEN)*(matmax + 1));
+  reallocate(matmax, (GEN*)&mat, &first_nz);
   setlg(mat, KCCO+1);
   C = cgetg(KCCO+1,t_MAT);
   /* trivial relations (p) = prod P^e */
   cglob = 0; z = zerocol(RU);
-  for (i=1; i<=KCZ; i++)
+  for (i=1; i<=F.KCZ; i++)
   {
-    GEN P = idealbase[i];
-    if (isclone(P))
-    { /* all prime divisors in FB */
-      unsetisclone(P); cglob++;
-      C[cglob] = (long)z; /* 0 */
-      mat[cglob] = p1 = col_0(KC);
-      k = numideal[FB[i]];
-      p1[0] = k+1; /* for already_found_relation */
-      p1 += k;
-      for (j=lg(P)-1; j; j--) p1[j] = itos(gmael(P,j,3));
-    }
+    long p = F.FB[i];
+    GEN c, P = F.LV[p];
+    if (!isclone(P)) continue;
+
+    /* all prime divisors in FB */
+    cglob++;
+    C[cglob] = (long)z; /* 0 */
+    mat[cglob] = c = col_0(F.KC);
+    k = F.iLP[p];
+    first_nz[cglob] = k+1;
+    c += k;
+    for (j=lg(P)-1; j; j--) c[j] = itos(gmael(P,j,3));
   }
   if (DEBUGLEVEL) fprintferr("After trivial relations, cglob = %ld\n",cglob);
   /* initialize for other relations */
   for (i=cglob+1; i<=KCCO; i++)
   {
-    mat[i] = col_0(KC);
+    mat[i] = col_0(F.KC);
     C[i] = (long)cgetc_col(RU,PRECREG);
   }
-  av1 = avma; liste = cgetg(KC+1,t_VECSMALL);
-  for (i=1; i<=KC; i++) liste[i] = 0;
+  av1 = avma; liste = vecsmall_const(F.KC, 0);
   invp = relationrank(mat,cglob,liste);
 
   /* relations through elements of small norm */
-  cglob = small_norm_for_buchall(cglob,mat,C,(long)LIMC,PRECREG,
-                                 nf,gborne,nbrelpid,invp,liste,LLLnf);
+  if (gsigne(gborne) > 0)
+    cglob = small_norm_for_buchall(cglob,mat,first_nz,C,(long)LIMC,PRECREG,&F,
+                                   nf,nbrelpid,invp,liste);
   if (cglob < 0) { precpb = "small_norm"; goto START; }
   avma = av1; limpile = stack_lim(av1,1);
 
   phase = 0;
-  nlze = matmax = 0; /* for lint */
-  vecT2 = NULL;
+  nlze = 0; /* for lint */
+  vecG = NULL;
+  list_jideal = NULL;
 
   /* random relations */
   if (cglob == KCCO) /* enough rels, but init random_relations just in case */
@@ -2904,6 +3064,7 @@ START:
   else
   {
     GEN matarch;
+    long ss;
 
     if (DEBUGLEVEL) fprintferr("\n#### Looking for random relations\n");
 MORE:
@@ -2911,55 +3072,57 @@ MORE:
     {
       if (DEBUGLEVEL) fprintferr("*** Increasing sub factor base\n");
       sfb_increase = 0;
-      if (++sfb_trials >= SFB_MAX) goto START;
-      subFB = subFBgen(N,min(lim,LIMC2), lg(subFB)-1+SFB_STEP,NULL,&ss);
-      if (!subFB) goto START;
+      if (++sfb_trials > SFB_MAX || !subFBgen_increase(&F, nf, SFB_STEP))
+        goto START;
       nreldep = nrelsup = 0;
     }
 
     if (phase == 0) matarch = C;
     else
     { /* reduced the relation matrix at least once */
-      long lgex = max(nlze, MIN_EXTRA); /* number of new relations sought */
-      long slim; /* total number of relations (including lgex new ones) */
+      long lgex = max(nlze, MIN_EXTRA); /* # of new relations sought */
+      long slim; /* total # of relations (including lgex new ones) */
       setlg(extraC,  lgex+1);
       setlg(extramat,lgex+1); /* were allocated after hnfspec */
       slim = cglob+lgex;
       if (slim > matmax)
       {
-        mat = (long**)gprealloc(mat, (2*slim+1) * sizeof(long*),
-                                     (matmax+1) * sizeof(long*));
         matmax = 2 * slim;
+        reallocate(matmax, (GEN*)&mat, &first_nz);
       }
       setlg(mat, slim+1);
       if (DEBUGLEVEL)
 	fprintferr("\n(need %ld more relation%s)\n", lgex, (lgex==1)?"":"s");
-      for (j=cglob+1; j<=slim; j++) mat[j] = col_0(KC);
+      for (j=cglob+1; j<=slim; j++) mat[j] = col_0(F.KC);
       matarch = extraC - cglob; /* start at 0, the others at cglob */
     }
-    if (!vecT2)
+    if (!vecG)
     {
-      vecT2 = compute_vecT2(RU,nf);
+      vecG = compute_vecG(nf,PRECLLL);
       av1 = avma;
     }
-    if (!powsubFB)
+    if (!F.powsubFB)
     {
-      powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);
+      powsubFBgen(&F, nf, CBUCHG+1, PRECREG);
       av1 = avma;
     }
-    ss = random_relation(phase,cglob,(long)LIMC,PRECLLL,PRECREG,
-                         nf,subFB,vecT2,mat,matarch,list_jideal);
+    ss = random_relation(phase,cglob,(long)LIMC,PRECREG,MAXRELSUP,
+                         nf,vecG,mat,first_nz,matarch,list_jideal,&F);
     if (ss < 0)
     { /* could not find relations */
-      if (ss != -1) precpb = "random_relation"; /* precision pb */
+      if (ss != -1)
+      {
+        precpb = "random_relation"; /* precision pb */
+        PRECLLLadd = (PRECLLLadd<<1) - 2;
+      }
       goto START;
     }
-    if (DEBUGLEVEL > 2) dbg_outrel(phase,cglob,vperm,mat,matarch);
-    if (phase) 
+    if (DEBUGLEVEL > 2 && phase == 0) dbg_outrel(cglob,mat,matarch);
+    if (phase)
       for (j=lg(extramat)-1; j>0; j--)
       {
 	GEN c = mat[cglob+j], *g = (GEN*) extramat[j];
-	for (k=1; k<=KC; k++) g[k] = stoi(c[vperm[k]]);
+	for (k=1; k<=F.KC; k++) g[k] = stoi(c[F.perm[k]]);
       }
     cglob = ss;
   }
@@ -2968,21 +3131,17 @@ MORE:
   if (phase == 0)
   { /* never reduced before */
     long lgex;
-    W = hnfspec(mat,vperm,&pdep,&B,&C, lg(subFB)-1);
+    W = hnfspec(mat,F.perm,&pdep,&B,&C, lg(F.subFB)-1);
     phase = 1;
     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];
-    }
+    if (nlze) list_jideal = vecextract_i(F.perm, 1, nlze);
     lgex = max(nlze, MIN_EXTRA); /* set lgex > 0, in case regulator is 0 */
     /* allocate now, reduce dimension later (setlg) when lgex goes down */
     extramat= cgetg(lgex+1,t_MAT);
     extraC  = cgetg(lgex+1,t_MAT);
     for (j=1; j<=lgex; j++)
     {
-      extramat[j]=lgetg(KC+1,t_COL);
+      extramat[j]=lgetg(F.KC+1,t_COL);
       extraC[j] = (long)cgetc_col(RU,PRECREG);
     }
   }
@@ -2990,94 +3149,76 @@ MORE:
   {
     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);
+      gerepileall(av1,6, &W,&C,&B,&pdep,&extramat,&extraC);
     }
     list_jideal = NULL;
-    W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);
+    W = hnfadd(W,F.perm,&pdep,&B,&C, extramat,extraC);
     nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
     if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto MORE; }
   }
   if (nlze) goto MORE; /* dependent rows */
 
   /* first attempt at regulator for the check */
-  sreg = (lg(mat)-1) - (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,&sublambda);
-  if (!reg)
+  zc = (lg(mat)-1) - (lg(B)-1) - (lg(W)-1);
+  A = vecextract_i(C, 1, zc); /* cols corresponding to units */
+  R = compute_multiple_of_R(A,RU,N,&lambda);
+  if (!R)
   { /* not full rank for units */
     if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
     if (++nrelsup > MAXRELSUP) goto START;
     nlze = MIN_EXTRA; goto MORE;
   }
   /* anticipate precision problems */
-  if (!sublambda) { precpb = "bestappr"; goto START; }
+  if (!lambda) { precpb = "bestappr"; goto START; }
 
-  /* class number */
-  clh = dethnf_i(W);
-  if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", clh);
+  h = dethnf_i(W);
+  if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h);
 
-  /* 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);
-  /* precision problems? */
-  if (!c1 || gcmpgs(gmul2n(c1,1),3) < 0)
-  { /* has to be a precision problem unless we cheat on Bach constant */
-    if (!precdouble) precpb = "compute_check";
-    goto START;
-  }
-  if (gcmpgs(c1,3) > 0)
+  z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */
+  switch (compute_R(lambda, divir(h,z), &L, &R))
   {
-    if (++nrelsup <= MAXRELSUP)
-    {
-      if (DEBUGLEVEL) fprintferr("\n ***** check = %f\n",gtodouble(c1)/2);
-      nlze = MIN_EXTRA; goto MORE;
-    }
-    if (cbach < 11.99) { sfb_increase = 1; goto MORE; }
-    err(warner,"suspicious check. Try to increase extra relations");
+    case PRECI: /* precision problem unless we cheat on Bach constant */
+      if (!precdouble) precpb = "compute_R";
+      goto START;
+    case RELAT: /* not enough relations */
+      if (++nrelsup <= MAXRELSUP) nlze = MIN_EXTRA; else sfb_increase = 1;
+      goto MORE;
   }
+  /* DONE */
 
-  if (KCZ2 > KCZ)
-  { /* "be honest" */
-    if (!powsubFB) powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);
-    if (!be_honest(nf,subFB,PRECLLL)) goto START;
-  }
+  if (!be_honest(&F, nf, PRECLLL)) goto START;
 
   /* fundamental units */
-  if (flun < 0 || flun > 1)
+  if (flun & nf_INIT || flun & nf_UNITS)
   {
-    xarch = cleancol(gmul(xarch,parch),N,PRECREG);
-    if (DEBUGLEVEL) msgtimer("cleancol");
+    GEN v = extract_full_lattice(L); /* L may be very large */
+    if (v)
+    {
+      A = vecextract_p(A, v);
+      L = vecextract_p(L, v);
+    }
+    /* arch. components of fund. units */
+    A = cleanarch(gmul(A,lllint(L)), N, PRECREG);
+    if (DEBUGLEVEL) msgtimer("cleanarch");
   }
-  if (labs(flun) > 1)
+  if (flun & nf_UNITS)
   {
-    fu = getfu(nf,&xarch,reg,flun,&k,PRECREG);
-    if (k <= 0 && labs(flun) > 2)
+    fu = getfu(nf,&A,R,flun,&k,PRECREG);
+    if (k <= 0 && flun & nf_FORCE)
     {
-      if (k < 0)
-      {
-        precadd = (DEFAULTPREC-2) + ((-k) >> TWOPOTBITS_IN_LONG);
-        if (precadd <= 0) precadd = (DEFAULTPREC-2);
-      }
+      if (k < 0) precadd = (DEFAULTPREC-2) + ((-k) >> TWOPOTBITS_IN_LONG);
+      (void)setrand(seed);
       precpb = "getfu"; goto START;
     }
   }
+  desallocate(&mat, &first_nz);
 
   /* class group generators */
-  i = lg(C)-sreg; C += sreg; C[0] = evaltyp(t_MAT)|evallg(i);
-  C = cleancol(C,N,PRECREG);
-  class_group_gen(nf,W,C,vperm, &clg1, &clg2, PRECREG);
-
-  c1 = gdiv(gmul(reg,clh), z);
-  desallocate(&mat);
-
-  return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c1,zu,W,B,xarch,C,vectbase,vperm));
+  i = lg(C)-zc; C += zc; C[0] = evaltyp(t_MAT)|evallg(i);
+  C = cleanarch(C,N,PRECREG);
+  Vbase = vecextract_p(F.LP, F.perm);
+  class_group_gen(nf,W,C,Vbase,PRECREG,NULL, &clg1, &clg2);
+  res = get_clfu(clg1, R, gdiv(mpmul(R,h), z), zu, fu, k);
+  return gerepilecopy(av, buchall_end(nf,CHANGE,flun,res,clg2,W,B,A,C,Vbase));
 }