===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/buch1.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/buch1.c	2001/10/02 11:17:03	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/buch1.c	2002/09/11 07:26:50	1.2
@@ -1,4 +1,4 @@
-/* $Id: buch1.c,v 1.1 2001/10/02 11:17:03 noro Exp $
+/* $Id: buch1.c,v 1.2 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -23,30 +23,45 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
 const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */
 
+/* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
+ * 2 | index), hence the low order bit is not useful. So we hash
+ * HASHBITS bits starting at bit 1, not bit 0 */
+#define HASHBITS 10
+const int HASHT = 1 << HASHBITS;
+
+static long
+hash(long q) { return (q & ((1 << (HASHBITS+1)) - 1)) >> 1; }
+#undef HASHBITS
+
 /* See buch2.c:
  * precision en digits decimaux=2*(#digits decimaux de Disc)+50
- * on prendra les p decomposes tels que prod(p)>lim dans la subbase
+ * on prendra les p decomposes tels que prod(p)>lim dans la subFB
  * LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
  * LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
- * subbase contient les p decomposes tels que prod(p)>sqrt(Disc)
- * lgsub=subbase[0]=#subbase;
- * subfactorbase est la table des form[p] pour p dans subbase
- * nbram est le nombre de p divisant Disc elimines dans subbase
- * powsubfactorbase est la table des puissances des formes dans subfactorbase
+ * subFB contient les indices des p decomposes tels que prod(p) > sqrt(Disc)
+ * powsubFB est la table des puissances des formes dans subFB
  */
-#define HASHT 1024 /* power of 2 */
-static const long CBUCH = 15; /* of the form 2^k-1 */
-static const long randshift=BITS_IN_RANDOM-1 - 4; /* BITS_IN_RANDOM-1-k */
 
-static long KC,KC2,lgsub,limhash,RELSUP,PRECREG;
-static long *primfact,*primfact1, *exprimfact,*exprimfact1, *badprim;
-static long *factorbase,*numfactorbase, *subbase, *vectbase, **hashtab;
-static GEN  **powsubfactorbase,vperm,subfactorbase,Disc,sqrtD,isqrtD;
+#define RANDOM_BITS 4
+static const long CBUCH = (1<<RANDOM_BITS)-1;
+static const long randshift=BITS_IN_RANDOM-1 - RANDOM_BITS;
+#undef RANDOM_BITS
 
+static long KC,KC2,limhash,RELSUP,PRECREG;
+static long *primfact,*exprimfact,*badprim;
+static long *FB,*numFB, **hashtab;
+static GEN  powsubFB,vperm,subFB,Disc,sqrtD,isqrtD;
+
 GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec);
 extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
-extern GEN roots_to_pol(GEN L, long v);
-extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);
+extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
+extern GEN quadhilbertreal(GEN D, GEN flag, long prec);
+extern void comp_gen(GEN z,GEN x,GEN y);
+extern GEN codeform5(GEN x, long prec);
+extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD);
+extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
+extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
+extern GEN cgetalloc(GEN x, size_t l, long t);
 
 GEN
 quadclassunit0(GEN x, long flag, GEN data, long prec)
@@ -109,7 +124,7 @@ getallforms(GEN D, long *pth, GEN *ptz)
 
 #define MOD4(x) ((((GEN)x)[2])&3)
 #define MAXL 300
-/* find P and Q two non principal prime ideals (above p,q) such that 
+/* find P and Q two non principal prime ideals (above p,q) such that
  *   (pq, z) = 1  [coprime to all class group representatives]
  *   cl(P) = cl(Q) if P has order 2 in Cl(K)
  * Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */
@@ -143,7 +158,7 @@ get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq)
   ell = 3;
   while (l < 3 || ell<=MAXL)
   {
-    ell += *diffell++; if (!*diffell) err(primer1);
+    NEXT_PRIME_VIADIFF_CHECK(ell, diffell);
     if (smodis(z,ell) && kross(d,ell) > 0)
     {
       court[2]=ell; form = redimag(primeform(D,court,0));
@@ -210,11 +225,12 @@ gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, lo
 static GEN
 quadhilbertimag(GEN D, GEN flag)
 {
-  long av=avma,h,i,e,prec;
+  long h, i, e, prec;
+  gpmem_t av=avma;
   GEN z,L,P,p,q,qfp,qfq,up,uq,u;
   int raw = ((typ(flag)==t_INT && signe(flag)));
 
-  if (DEBUGLEVEL>=2) timer2();
+  if (DEBUGLEVEL>=2) (void)timer2();
   if (gcmpgs(D,-11) >= 0) return polx[0];
   L = getallforms(D,&h,&z);
   if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h);
@@ -237,9 +253,10 @@ quadhilbertimag(GEN D, GEN flag)
   prec = raw? DEFAULTPREC: 3;
   for(;;)
   {
-    long av0 = avma, ex, exmax = 0;
+    long ex, exmax = 0;
+    gpmem_t av0 = avma;
     GEN lead, sqd = gsqrt(negi(D),prec);
-    P = cgetg(h+1,t_VEC); 
+    P = cgetg(h+1,t_VEC);
     for (i=1; i<=h; i++)
     {
       GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec);
@@ -250,7 +267,7 @@ quadhilbertimag(GEN D, GEN flag)
         v[2] = (long)s;
       }
       else P[i] = (long)s;
-      ex = gexpo(s); if (ex > 0) exmax += ex; 
+      ex = gexpo(s); if (ex > 0) exmax += ex;
     }
     if (DEBUGLEVEL>=2) msgtimer("roots");
     if (raw) { P = gcopy(P); break; }
@@ -271,8 +288,6 @@ quadhilbertimag(GEN D, GEN flag)
   return gerepileupto(av,P);
 }
 
-GEN quadhilbertreal(GEN D, GEN flag, long prec);
-
 GEN
 quadhilbert(GEN D, GEN flag, long prec)
 {
@@ -310,11 +325,14 @@ get_om(GEN nf, GEN a)
  * Set list[j + 1] = g1^e1...gk^ek where j is the integer
  *   ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
 static GEN
-getallelts(GEN nf, GEN G)
+getallelts(GEN bnr)
 {
-  GEN C,c,g, *list, **pows, *gk;
+  GEN nf,G,C,c,g, *list, **pows, *gk;
   long lc,i,j,k,no;
 
+  nf = checknf(bnr);
+  G = (GEN)bnr[5];
+
   no = itos((GEN)G[1]);
   c = (GEN)G[2];
   g = (GEN)G[3]; lc = lg(c)-1;
@@ -330,11 +348,12 @@ getallelts(GEN nf, GEN G)
   {
     c[i] = k = itos((GEN)c[i]);
     gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i];
-    for (j=2; j<k; j++) gk[j] = idealmul(nf, gk[j-1], gk[1]);
+    for (j=2; j<k; j++)
+      gk[j] = idealmodidele(bnr, idealmul(nf, gk[j-1], gk[1]));
     pows[i] = gk; /* powers of g[i] */
   }
 
-  C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc]; 
+  C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
   for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
   /* C[i] = c(k-i+1) * ... * ck */
   /* j < C[i+1] <==> j only involves g(k-i)...gk */
@@ -345,8 +364,9 @@ getallelts(GEN nf, GEN G)
   {
     GEN p1,p2;
     if (j == C[i+1]) i++;
-    p2 = pows[lc-i][j/C[i]]; 
-    p1 = list[j%C[i] + 1]; if (p1) p2 = idealmul(nf,p2,p1);
+    p2 = pows[lc-i][j/C[i]];
+    p1 = list[j%C[i] + 1];
+    if (p1) p2 = idealmodidele(bnr, idealmul(nf,p2,p1));
     list[j + 1] = p2;
   }
   list[1] = idealhermite(nf,gun);
@@ -419,14 +439,6 @@ get_prec(GEN P, long prec)
 }
 
 /* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */
-GEN
-PiI2(long prec)
-{
-  GEN z = cgetg(3,t_COMPLEX);
-  GEN x = mppi(prec); setexpo(x,2);
-  z[1] = zero;
-  z[2] = (long)x; return z; /* = 2I Pi */
-}
 
 /* Compute data for ellphist */
 static GEN
@@ -436,9 +448,10 @@ ellphistinit(GEN om, long prec)
 
   if (gsigne(gimag(gdiv(om1,om2))) < 0)
   {
-    p1=om1; om1=om2; om2=p1;
-    p1=cgetg(3,t_VEC); p1[1]=(long)om1; p1[2]=(long)om2;
-    om = p1;
+    p1 = om1; om1 = om2; om2 = p1;
+    om = cgetg(3,t_VEC);
+    om[1] = (long)om1;
+    om[2] = (long)om2;
   }
   om1b = gconj(om1);
   om2b = gconj(om2); res = cgetg(4,t_VEC);
@@ -475,12 +488,13 @@ computeth2(GEN om, GEN la, long prec)
 static GEN
 computeP2(GEN bnr, GEN la, int raw, long prec)
 {
-  long av=avma, av2, clrayno,i, first = 1;
+  long clrayno, i, first = 1;
+  gpmem_t av=avma, av2;
   GEN listray,P0,P,f,lanum, nf = checknf(bnr);
-  
+
   f = gmael3(bnr,2,1,1);
   if (typ(la) != t_COL) la = algtobasis(nf,la);
-  listray = getallelts(nf,(GEN)bnr[5]);
+  listray = getallelts(bnr);
   clrayno = lg(listray)-1; av2 = avma;
 PRECPB:
   if (!first)
@@ -518,7 +532,7 @@ do_compo(GEN x, GEN y)
 {
   long a, ph = degpol(y);
   GEN z;
-  y = gmul(gpuigs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0])));
+  y = gmul(gpowgs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0])));
   for  (a = 0;; a = nexta(a))
   {
     if (a) x = gsubst(x, 0, gaddsg(a, polx[0]));
@@ -543,7 +557,8 @@ galoisapplypol(GEN nf, GEN s, GEN x)
 static GEN
 findquad(GEN a, GEN x, GEN p)
 {
-  long tu,tv, av = avma;
+  long tu, tv;
+  gpmem_t av = avma;
   GEN u,v;
   if (typ(x) == t_POLMOD) x = (GEN)x[2];
   if (typ(a) == t_POLMOD) a = (GEN)a[2];
@@ -645,14 +660,15 @@ treatspecialsigma(GEN nf, GEN gf, int raw, long prec)
     return NULL;
   }
 
-  p1 = gcoeff(gf,1,1);
+  p1 = gcoeff(gf,1,1); /* integer > 0 */
   p2 = gcoeff(gf,2,2);
   if (gcmp1(p2)) { fl = 0; tryf = p1; }
   else {
     if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL;
     fl = 1; tryf = shifti(p1,-1);
   }
-  if (cmpis(tryf, 3) <= 0 || !gcmp0(resii(D, tryf)) || !isprime(tryf))
+  /* tryf integer > 0 */
+  if (cmpis(tryf, 3) <= 0 || signe(resii(D, tryf)) || !isprime(tryf))
     return NULL;
 
   i = itos(tryf); if (fl) i *= 4;
@@ -701,7 +717,7 @@ quadrayimagsigma(GEN bnr, int raw, long prec)
   f2 = 2 * itos(gcoeff(f,1,1));
   u = getallrootsof1(bnf); lu = lg(u);
   for (i=1; i<lu; i++)
-    u[i] = (long)colreducemodmat((GEN)u[i], f, NULL); /* roots of 1, mod f */
+    u[i] = (long)colreducemodHNF((GEN)u[i], f, NULL); /* roots of 1, mod f */
   if (DEBUGLEVEL>1)
     fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
   for (a=0; a<f2; a++)
@@ -712,7 +728,7 @@ quadrayimagsigma(GEN bnr, int raw, long prec)
       if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b);
 
       labas = algtobasis(nf, la);
-      lamodf = colreducemodmat(labas, f, NULL);
+      lamodf = colreducemodHNF(labas, f, NULL);
       for (i=1; i<lu; i++)
         if (gegal(lamodf, (GEN)u[i])) break;
       if (i < lu) continue; /* la = unit mod f */
@@ -731,7 +747,8 @@ GEN
 quadray(GEN D, GEN f, GEN flag, long prec)
 {
   GEN bnr,y,p1,pol,bnf,lambda;
-  long av = avma, raw;
+  long raw;
+  gpmem_t av = avma;
 
   if (!flag) flag = gzero;
   if (typ(flag)==t_INT) lambda = NULL;
@@ -781,27 +798,39 @@ quadray(GEN D, GEN f, GEN flag, long prec)
 /*  Routines related to binary quadratic forms (for internal use)  */
 /*                                                                 */
 /*******************************************************************/
-extern void comp_gen(GEN z,GEN x,GEN y);
-extern GEN codeform5(GEN x, long prec);
-extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD);
-extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
-extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
-
 #define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD)
 #define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD)))
-#define comprealform(x,y) fix_signs(comprealform5(x,y,Disc,sqrtD,isqrtD))
 
+static GEN
+fix_signs(GEN x)
+{
+  GEN a = (GEN)x[1];
+  GEN c = (GEN)x[3];
+  if (signe(a) < 0)
+  {
+    if (narrow || absi_equal(a,c)) return rhorealform(x);
+    setsigne(a,1); setsigne(c,-1);
+  }
+  return x;
+}
+
+static GEN
+comprealform(GEN x,GEN y) {
+  return fix_signs( comprealform5(x,y,Disc,sqrtD,isqrtD) );
+}
+
 /* compute rho^n(x) */
-static GEN 
+static GEN
 rhoreal_pow(GEN x, long n)
 {
-  long i, av = avma, lim = stack_lim(av,1);
+  long i;
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   for (i=1; i<=n; i++)
   {
     x = rhorealform(x);
     if (low_stack(lim, stack_lim(av,1)))
     {
-      if(DEBUGMEM>1) err(warnmem,"rhoreal_pow");	
+      if(DEBUGMEM>1) err(warnmem,"rhoreal_pow");
       x = gerepilecopy(av, x);
     }
   }
@@ -809,54 +838,71 @@ rhoreal_pow(GEN x, long n)
 }
 
 static GEN
-fix_signs(GEN x)
+realpf5(GEN D, long p)
 {
-  GEN a = (GEN)x[1];
-  GEN c = (GEN)x[3];
-  if (signe(a) < 0)
-  {
-    if (narrow || absi_equal(a,c)) return rhorealform(x);
-    setsigne(a,1); setsigne(c,-1);
-  }
-  return x;
-}
-
-static GEN
-redrealprimeform5(GEN Disc, long p)
-{
-  long av = avma;
-  GEN y = primeform(Disc,stoi(p),PRECREG);
+  gpmem_t av = avma;
+  GEN y = primeform(D,stoi(p),PRECREG);
   y = codeform5(y,PRECREG);
   return gerepileupto(av, redrealform(y));
 }
 
 static GEN
-redrealprimeform(GEN Disc, long p)
+realpf(GEN D, long p)
 {
-  long av = avma;
-  GEN y = primeform(Disc,stoi(p),PRECREG);
+  gpmem_t av = avma;
+  GEN y = primeform(D,stoi(p),PRECREG);
   return gerepileupto(av, redrealform(y));
 }
 
 static GEN
+imagpf(GEN D, long p) { return primeform(D,stoi(p),0); }
+
+static GEN
 comprealform3(GEN x, GEN y)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN z = cgetg(4,t_VEC); comp_gen(z,x,y);
   return gerepileupto(av, redrealform(z));
 }
 
+/* Warning: ex[0] not set */
 static GEN
-initrealform5(long *ex)
+init_form(long *ex, GEN (*comp)(GEN,GEN))
 {
-  GEN form = powsubfactorbase[1][ex[1]];
-  long i;
+  long i, l = lg(powsubFB);
+  GEN F = NULL;
+  for (i=1; i<l; i++)
+    if (ex[i])
+    {
+      GEN t = gmael(powsubFB,i,ex[i]);
+      F = F? comp(F,t): t;
+    }
+  return F;
+}
+static GEN
+initrealform5(long *ex) { return init_form(ex, &comprealform); }
+static GEN
+initimagform(long *ex) { return init_form(ex, &compimag); }
 
-  for (i=2; i<=lgsub; i++)
-    form = comprealform(form, powsubfactorbase[i][ex[i]]);
-  return form;
+static GEN
+random_form(GEN ex, GEN (*comp)(GEN,GEN))
+{
+  long i, l = lg(ex);
+  gpmem_t av = avma;
+  GEN F;
+  for(;;)
+  {
+    for (i=1; i<l; i++) ex[i] = mymyrand()>>randshift;
+    if ((F = init_form(ex, comp))) return F;
+    avma = av;
+  }
 }
 
+static GEN
+real_random_form(GEN ex) { return random_form(ex, &comprealform3); }
+static GEN
+imag_random_form(GEN ex) { return random_form(ex, &compimag); }
+
 /*******************************************************************/
 /*                                                                 */
 /*                     Common subroutines                          */
@@ -865,11 +911,9 @@ initrealform5(long *ex)
 static void
 buch_init(void)
 {
-  if (DEBUGLEVEL) timer2();
+  if (DEBUGLEVEL) (void)timer2();
   primfact  = new_chunk(100);
-  primfact1 = new_chunk(100);
   exprimfact  = new_chunk(100);
-  exprimfact1 = new_chunk(100);
   badprim = new_chunk(100);
   hashtab = (long**) new_chunk(HASHT);
 }
@@ -885,9 +929,10 @@ check_bach(double cbach, double B)
 }
 
 static long
-factorisequad(GEN f, long kcz, long limp)
+factorquad(GEN f, long kcz, long limp)
 {
-  long i,p,k,av,lo;
+  long i, p, k, lo;
+  gpmem_t av;
   GEN q,r, x = (GEN)f[1];
 
   if (is_pm1(x)) { primfact[0]=0; return 1; }
@@ -895,24 +940,23 @@ factorisequad(GEN f, long kcz, long limp)
   if (signe(x) < 0) x = absi(x);
   for (i=1; ; i++)
   {
-    p=factorbase[i]; q=dvmdis(x,p,&r);
+    p=FB[i]; q=dvmdis(x,p,&r);
     if (!signe(r))
     {
-      k=0; while (!signe(r)) { x=q; k++; q=dvmdis(x,p,&r); }
-      lo++; primfact[lo]=p; exprimfact[lo]=k;
+      for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
+      primfact[++lo]=p; exprimfact[lo]=k;
     }
     if (cmpis(q,p)<=0) break;
     if (i==kcz) { avma=av; return 0; }
   }
-  p = x[2]; avma=av;
-  /* p = itos(x) if lgefint(x)=3 */
-  if (lgefint(x)!=3 || p > limhash) return 0;
+  avma=av;
+  if (lgefint(x)!=3 || (p=x[2]) > limhash) return 0;
 
   if (p != 1 && p <= limp)
   {
     for (i=1; i<=badprim[0]; i++)
       if (p % badprim[i] == 0) return 0;
-    lo++; primfact[lo]=p; exprimfact[lo]=1;
+    primfact[++lo]=p; exprimfact[lo]=1;
     p = 1;
   }
   primfact[0]=lo; return p;
@@ -922,219 +966,207 @@ factorisequad(GEN f, long kcz, long limp)
 static long *
 largeprime(long q, long *ex, long np, long nrho)
 {
-  const long hashv = ((q & (2 * HASHT - 1)) - 1) >> 1;
-  long *pt, i;
+  const long hashv = hash(q);
+  long *pt, i, l = lg(subFB);
 
- /* If q = 0 (2048), very slim chance of getting a relation.
-  * And hashtab[-1] is undefined anyway */
-  if (hashv < 0) return NULL;
   for (pt = hashtab[hashv]; ; pt = (long*) pt[0])
   {
     if (!pt)
     {
-      pt = (long*) gpmalloc((lgsub+4)<<TWOPOTBYTES_IN_LONG);
+      pt = (long*) gpmalloc((l+3)<<TWOPOTBYTES_IN_LONG);
       *pt++ = nrho; /* nrho = pt[-3] */
       *pt++ = np;   /* np   = pt[-2] */
       *pt++ = q;    /* q    = pt[-1] */
       pt[0] = (long)hashtab[hashv];
-      for (i=1; i<=lgsub; i++) pt[i]=ex[i];
+      for (i=1; i<l; i++) pt[i]=ex[i];
       hashtab[hashv]=pt; return NULL;
     }
     if (pt[-1] == q) break;
   }
-  for(i=1; i<=lgsub; i++)
+  for(i=1; i<l; i++)
     if (pt[i] != ex[i]) return pt;
   return (pt[-2]==np)? (GEN)NULL: pt;
 }
 
-static long
-badmod8(GEN x)
+/* p | conductor of order of disc D ? */
+static int
+is_bad(GEN D, ulong p)
 {
-  long r = mod8(x);
-  if (!r) return 1;
-  if (signe(Disc) < 0) r = 8-r;
-  return (r < 4);
+  gpmem_t av = avma;
+  int r;
+  if (p == 2)
+  {
+    r = mod16(D) >> 1;
+    if (r && signe(D) < 0) r = 8-r;
+    return (r < 4);
+  }
+  r = (resii(D, muluu(p,p)) == gzero); /* p^2 | D ? */
+  avma = av; return r;
 }
 
-/* cree factorbase, numfactorbase, vectbase; affecte badprim */
-static void
-factorbasequad(GEN Disc, long n2, long n)
+/* create FB, numFB; fill badprim. Return L(kro_D, 1) */
+static GEN
+FBquad(GEN Disc, long n2, long n)
 {
-  long i,p,bad, av = avma;
-  byteptr d=diffptr;
+  GEN Res = realun(DEFAULTPREC);
+  long i, p, bad, s;
+  gpmem_t av;
+  byteptr d = diffptr;
 
-  numfactorbase = (long*) gpmalloc(sizeof(long)*(n2+1));
-  factorbase    = (long*) gpmalloc(sizeof(long)*(n2+1));
-  KC=0; bad=0; i=0; p = *d++;
-  while (p<=n2)
+  numFB = cgetg(n2+1, t_VECSMALL);
+  FB    = cgetg(n2+1, t_VECSMALL);
+  av = avma;
+  KC = 0; bad = 0; i = 0;
+  if (maxprime() < n2) err(primer1);
+  for (p = 0;;) /* p <= n2 */
   {
-    switch(krogs(Disc,p))
+    NEXT_PRIME_VIADIFF(p, d);
+    if (!KC && p > n) KC = i;
+    if (p > n2) break;
+    s = krogs(Disc,p);
+    switch (s)
     {
       case -1: break; /* inert */
       case  0: /* ramified */
-      {
-        GEN p1 = divis(Disc,p);
-	if (smodis(p1,p) == 0)
-          if (p!=2 || badmod8(p1)) { badprim[++bad]=p; break; }
-        i++; numfactorbase[p]=i; factorbase[i] = -p; break;
-      }
+        if (is_bad(Disc, (ulong)p)) { badprim[++bad]=p; break; }
+        /* fall through */
       default:  /* split */
-        i++; numfactorbase[p]=i; factorbase[i] = p;
+        i++; numFB[p]=i; FB[i] = p; break;
     }
-    p += *d++; if (!*d) err(primer1);
-    if (KC == 0 && p>n) KC = i;
+    Res = mulsr(p, divrs(Res, p - s));
   }
-  if (!KC) { free(factorbase); free(numfactorbase); return; }
+  if (!KC) return NULL;
+  Res = gerepileupto(av, Res);
   KC2 = i;
-  vectbase = (long*) gpmalloc(sizeof(long)*(KC2+1));
-  for (i=1; i<=KC2; i++)
-  {
-    p = factorbase[i];
-    vectbase[i]=p; factorbase[i]=labs(p);
-  }
+  setlg(FB, KC2+1);
   if (DEBUGLEVEL)
   {
     msgtimer("factor base");
     if (DEBUGLEVEL>7)
     {
-      fprintferr("factorbase:\n");
-      for (i=1; i<=KC; i++) fprintferr("%ld ",factorbase[i]);
-      fprintferr("\n"); flusherr();
+      fprintferr("FB:\n");
+      for (i=1; i<=KC; i++) fprintferr("%ld ", FB[i]);
+      fprintferr("\n");
     }
   }
-  avma=av; badprim[0] = bad;
+  badprim[0] = bad; return Res;
 }
 
-/* cree vectbase and subfactorbase. Affecte lgsub */
-static long
-subfactorbasequad(double ll, long KC)
+/* create vperm, return subFB */
+static GEN
+subFBquad(GEN D, double PROD, long KC)
 {
-  long i,j,k,nbidp,p,pro[100], ss;
-  GEN y;
-  double prod;
+  long i, j, lgsub, ino, lv = KC+1;
+  double prod = 1.;
+  gpmem_t av;
+  GEN no;
 
-  i=0; ss=0; prod=1;
-  for (j=1; j<=KC && prod<=ll; j++)
+  vperm = cgetg(lv, t_VECSMALL);
+  av = avma;
+  no    = cgetg(lv, t_VECSMALL); ino = 1;
+  for (i=j=1; j < lv; j++)
   {
-    p = vectbase[j];
-    if (p>0) { pro[++i]=p; prod*=p; vperm[i]=j; } else ss++;
+    long p = FB[j];
+    if (smodis(D, p) == 0) no[ino++] = j; /* ramified */
+    else
+    {
+      vperm[i] = j; i++;
+      prod *= p;
+      if (prod > PROD) break;
+    }
   }
-  if (prod<=ll) return -1;
-  nbidp=i;
-  for (k=1; k<j; k++)
-    if (vectbase[k]<=0) vperm[++i]=k;
-
-  y=cgetg(nbidp+1,t_COL);
-  if (PRECREG) /* real */
-    for (j=1; j<=nbidp; j++)
-      y[j] = (long)redrealprimeform5(Disc, pro[j]);
-  else
-    for (j=1; j<=nbidp; j++) /* imaginary */
-      y[j] = (long)primeform(Disc,stoi(pro[j]),0);
-  subbase = (long*) gpmalloc(sizeof(long)*(nbidp+1));
-  lgsub = nbidp; for (j=1; j<=lgsub; j++) subbase[j]=pro[j];
-  if (DEBUGLEVEL>7)
-  {
-    fprintferr("subfactorbase: ");
-    for (i=1; i<=lgsub; i++)
-      { fprintferr("%ld: ",subbase[i]); outerr((GEN)y[i]); }
-    fprintferr("\n"); flusherr();
-  }
-  subfactorbase = y; return ss;
+  if (j == lv) return NULL;
+  lgsub = i;
+  for (j = 1; j <ino; i++,j++) vperm[i] = no[j];
+  for (     ; i < lv; i++)     vperm[i] = i;
+  avma = av;
+  if (DEBUGLEVEL) msgtimer("subFBquad (%ld elt.)", lgsub-1);
+  return vecextract_i(vperm, 1, lgsub-1);
 }
 
-static void
-powsubfact(long n, long a)
+/* assume n >= 1, x[i][j] = subFB[i]^j, for j = 1..n */
+static GEN
+powsubFBquad(long n)
 {
-  GEN unform, *y, **x = (GEN**) gpmalloc(sizeof(GEN*)*(n+1));
-  long i,j;
+  long i,j, l = lg(subFB);
+  GEN F, y, x = cgetg(l, t_VEC);
 
-  for (i=1; i<=n; i++)
-    x[i] = (GEN*) gpmalloc(sizeof(GEN)*(a+1));
   if (PRECREG) /* real */
   {
-    unform=cgetg(6,t_VEC);
-    unform[1]=un;
-    unform[2]=(mod2(Disc) == mod2(isqrtD))? (long)isqrtD: laddsi(-1,isqrtD);
-    unform[3]=lshifti(subii(sqri((GEN)unform[2]),Disc),-2);
-    unform[4]=zero;
-    unform[5]=(long)realun(PRECREG);
-    for (i=1; i<=n; i++)
+    for (i=1; i<l; i++)
     {
-      y = x[i]; y[0] = unform;
-      for (j=1; j<=a; j++)
-	y[j] = comprealform(y[j-1],(GEN)subfactorbase[i]);
+      F = realpf5(Disc, FB[subFB[i]]);
+      y = cgetg(n+1, t_VEC); x[i] = (long)y;
+      y[1] = (long)F;
+      for (j=2; j<=n; j++) y[j] = (long)comprealform((GEN)y[j-1], F);
     }
   }
   else /* imaginary */
   {
-    unform = primeform(Disc, gun, 0);
-    for (i=1; i<=n; i++)
+    for (i=1; i<l; i++)
     {
-      y = x[i]; y[0] = unform;
-      for (j=1; j<=a; j++)
-	y[j] = compimag(y[j-1],(GEN)subfactorbase[i]);
+      F = imagpf(Disc, FB[subFB[i]]);
+      y = cgetg(n+1, t_VEC); x[i] = (long)y;
+      y[1] = (long)F;
+      for (j=2; j<=n; j++) y[j] = (long)compimag((GEN)y[j-1], F);
     }
   }
-  if (DEBUGLEVEL) msgtimer("powsubfact");
-  powsubfactorbase = x;
+  if (DEBUGLEVEL) msgtimer("powsubFBquad");
+  return x;
 }
 
 static void
 desalloc(long **mat)
 {
   long i,*p,*q;
+  for (i=1; i<lg(mat); i++) free(mat[i]);
+  free(mat);
+  for (i=1; i<HASHT; i++)
+    for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }
+}
 
-  free(vectbase); free(factorbase); free(numfactorbase);
-  if (mat)
+static void
+sub_fact(GEN col, GEN F)
+{
+  GEN b = (GEN)F[2];
+  long i, p, e;
+  for (i=1; i<=primfact[0]; i++)
   {
-    free(subbase);
-    for (i=1; i<lg(subfactorbase); i++) free(powsubfactorbase[i]);
-    for (i=1; i<lg(mat); i++) free(mat[i]);
-    free(mat); free(powsubfactorbase);
-    for (i=1; i<HASHT; i++)
-      for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }
+    p = primfact[i]; e = exprimfact[i];
+    if (smodis(b, p<<1) > p) e = -e;
+    col[numFB[p]] -= e;
   }
 }
-
-/* L-function */
-static GEN
-lfunc(GEN Disc)
+static void
+add_fact(GEN col, GEN F)
 {
-  long av=avma, p;
-  GEN y=realun(DEFAULTPREC);
-  byteptr d=diffptr;
-
-  for(p = *d++; p<=30000; p += *d++)
+  GEN b = (GEN)F[2];
+  long i, p, e;
+  for (i=1; i<=primfact[0]; i++)
   {
-    if (!*d) err(primer1);
-    y = mulsr(p, divrs(y, p-krogs(Disc,p)));
+    p = primfact[i]; e = exprimfact[i];
+    if (smodis(b, p<<1) > p) e = -e;
+    col[numFB[p]] += e;
   }
-  return gerepileupto(av,y);
 }
 
 #define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y
 static GEN
 get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
 {
-  GEN p1,p2,res,*init, u1u2=smith2(W), u1=(GEN)u1u2[1], met=(GEN)u1u2[3];
-  long c,i,j, l = lg(met);
+  GEN res,*init, u1, met = smithrel(W,NULL,&u1);
+  long c,i,j, l = lg(W);
 
-  u1 = reducemodmatrix(ginv(u1), W);
-  for (c=1; c<l; c++)
-    if (gcmp1(gcoeff(met,c,c))) break;
+  c = lg(met);
   if (DEBUGLEVEL) msgtimer("smith/class group");
   res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC);
   for (i=1; i<l; i++)
-    init[i] = primeform(Disc,stoi(labs(vectbase[vperm[i]])),prec);
+    init[i] = primeform(Disc,stoi(FB[vperm[i]]),prec);
   for (j=1; j<c; j++)
   {
-    p1 = NULL;
+    GEN p1 = NULL;
     for (i=1; i<l; i++)
-    {
-      p2 = gpui(init[i], gcoeff(u1,i,j), prec);
-      p1 = comp(p1,p2);
-    }
+      p1 = comp(p1, powgi(init[i], gcoeff(u1,i,j)));
     res[j] = (long)p1;
   }
   if (DEBUGLEVEL) msgtimer("generators");
@@ -1142,47 +1174,46 @@ get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
 }
 
 static GEN
-extra_relations(long LIMC, long *ex, long nlze, GEN extramatc)
+extra_relations(long LIMC, long nlze, GEN *ptextraC)
 {
-  long av,fpc,p,ep,i,j,k,nlze2, *col, *colg, s = 0, extrarel = nlze+2;
+  long fpc, i, k, nlze2;
+  long s = 0, extrarel = nlze+2, lgsub = lg(subFB);
+  gpmem_t av;
   long MAXRELSUP = min(50,4*KC);
-  GEN p1,form, extramat = cgetg(extrarel+1,t_MAT);
+  GEN p1, form, ex, extramat, extraC, col;
 
-  if (DEBUGLEVEL)
-  {
-    fprintferr("looking for %ld extra relations\n",extrarel);
-    flusherr();
-  }
-  for (j=1; j<=extrarel; j++) extramat[j]=lgetg(KC+1,t_COL);
-  nlze2 = PRECREG? max(nlze,lgsub): min(nlze+1,KC);
+  extramat = cgetg(extrarel+1, t_MAT);
+  extraC   = cgetg(extrarel+1, t_VEC);
+  for (i=1; i<=extrarel; i++) extramat[i] = lgetg(KC+1,t_VECSMALL);
+  if (!PRECREG)
+    for (i=1; i<=extrarel; i++) extraC[i] = lgetg(1,t_COL);
+
+  if (DEBUGLEVEL) fprintferr("looking for %ld extra relations\n",extrarel);
+  nlze2 = PRECREG? max(nlze,lgsub-1): min(nlze+1,KC);
   if (nlze2 < 3 && KC > 2) nlze2 = 3;
+  ex = cgetg(nlze2+1, t_VECSMALL);
   av = avma;
   while (s<extrarel)
   {
     form = NULL;
     for (i=1; i<=nlze2; i++)
     {
-      ex[i]=mymyrand()>>randshift;
+      ex[i] = mymyrand()>>randshift;
       if (ex[i])
       {
-        p1 = primeform(Disc,stoi(factorbase[vperm[i]]),PRECREG);
-        p1 = gpuigs(p1,ex[i]); form = comp(form,p1);
+        p1 = primeform(Disc,stoi(FB[vperm[i]]),PRECREG);
+        p1 = gpowgs(p1,ex[i]); form = comp(form,p1);
       }
     }
     if (!form) continue;
 
-    fpc = factorisequad(form,KC,LIMC);
+    fpc = factorquad(form,KC,LIMC);
     if (fpc==1)
     {
       s++; col = (GEN)extramat[s];
       for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i];
-      for (   ; i<=KC; i++) col[vperm[i]]= 0;
-      for (j=1; j<=primfact[0]; j++)
-      {
-        p=primfact[j]; ep=exprimfact[j];
-        if (smodis((GEN)form[2], p<<1) > p) ep = -ep;
-        col[numfactorbase[p]] += ep;
-      }
+      for (   ; i<=KC;    i++) col[vperm[i]] = 0;
+      add_fact(col, form);
       for (i=1; i<=KC; i++)
         if (col[i]) break;
       if (i>KC)
@@ -1190,165 +1221,165 @@ extra_relations(long LIMC, long *ex, long nlze, GEN ex
         s--; avma = av;
         if (--MAXRELSUP == 0) return NULL;
       }
-      else { av = avma; if (PRECREG) coeff(extramatc,1,s) = form[4]; }
+      else { av = avma; if (PRECREG) extraC[s] = form[4]; }
     }
     else avma = av;
     if (DEBUGLEVEL)
     {
       if (fpc == 1) fprintferr(" %ld",s);
       else if (DEBUGLEVEL>1) fprintferr(".");
-      flusherr();
     }
   }
-  for (j=1; j<=extrarel; j++)
+
+  for (i=1; i<=extrarel; i++)
   {
-    colg = cgetg(KC+1,t_COL);
-    col = (GEN)extramat[j]; extramat[j] = (long) colg;
-    for (k=1; k<=KC; k++)
-      colg[k] = lstoi(col[vperm[k]]);
+    GEN colg = cgetg(KC+1,t_COL);
+    col = (GEN)extramat[i];
+    extramat[i] = (long)colg;
+    for (k=1; k<=KC; k++) colg[k] = lstoi(col[vperm[k]]);
   }
-  if (DEBUGLEVEL)
+  if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("extra relations"); }
+  *ptextraC = extraC; return extramat;
+}
+#undef comp
+
+static long
+trivial_relations(long **mat, GEN C)
+{
+  long i, s = 0;
+  GEN col;
+  for (i=1; i<=KC; i++) 
   {
-    fprintferr("\n");
-    msgtimer("extra relations");
+    if (smodis(Disc, FB[i])) continue;
+    /* ramified prime ==> trivial relation */
+    col = mat[++s];
+    col[i] = 2;
+    if (C) affsr(0, (GEN)C[s]);
   }
-  return extramat;
+  return s;
 }
-#undef comp
 
+static void
+dbg_all(long phase, long s, long n)
+{
+  if (DEBUGLEVEL>1) fprintferr("\n");
+  msgtimer("%s rel [#rel/#test = %ld/%ld]", phase? "random": "initial", s, n);
+}
+
+void
+wr_rel(GEN col)
+{
+  long i, l = lg(col);
+  fprintferr("\nrel = ");
+  for (i=1; i<l; i++)
+    if (col[i]) fprintferr("%ld^%ld ",i,col[i]);
+  fprintferr("\n");
+}
+
+void
+dbg_rel(long s, GEN col)
+{
+  if (DEBUGLEVEL == 1) fprintferr("%4ld",s);
+  else { fprintferr("cglob = %ld. ", s); wr_rel(col); }
+  flusherr(); 
+}
 /*******************************************************************/
 /*                                                                 */
 /*                    Imaginary Quadratic fields                   */
 /*                                                                 */
 /*******************************************************************/
 
+/* Strategy 1 until lim relation found, then Strategy 2 until LIM. */
 static GEN
-imag_random_form(long current, long *ex)
+imag_relations(long LIM, long lim, long LIMC, long **mat)
 {
-  long av = avma,i;
-  GEN form,pc;
+  long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0;
+  int first = 1;
+  gpmem_t av;
+  GEN C, col, form, ex = cgetg(lgsub, t_VECSMALL);
+  GEN dummy = cgetg(1,t_COL);
 
+  C = cgetg(LIM+1, t_VEC);
+  for (i=1; i<=LIM; i++) C[i] = (long)dummy;
+  av = avma;
+  s = trivial_relations(mat, NULL);
+  lim += s; if (lim > LIM) lim = LIM;
   for(;;)
   {
-    form = pc = primeform(Disc,stoi(factorbase[current]),PRECREG);
-    for (i=1; i<=lgsub; i++)
-    {
-      ex[i] = mymyrand()>>randshift;
-      if (ex[i])
-        form = compimag(form,powsubfactorbase[i][ex[i]]);
+    if (s >= lim) {
+      if (s >= LIM) break;
+      lim = LIM; first = 0;
+      if (DEBUGLEVEL) dbg_all(0, s, nbtest);
     }
-    if (form != pc) return form;
-    avma = av; /* ex = 0, try again */
-  }
-}
-
-static void
-imag_relations(long lim, long s, long LIMC, long *ex, long **mat)
-{
-  static long nbtest;
-  long av = avma, i,j,pp,fpc,b1,b2,ep,current, first = (s==0);
-  long *col,*fpd,*oldfact,*oldexp;
-  GEN pc,form,form1;
-
-  if (first) nbtest = 0 ;
-  while (s<lim)
-  {
-    avma=av; nbtest++; current = first? 1+(s%KC): 1+s-RELSUP;
-    form = imag_random_form(current,ex);
-    fpc = factorisequad(form,KC,LIMC);
+    avma = av; current = first? 1+(s%KC): 1+s-RELSUP;
+    form = imagpf(Disc, FB[current]);
+    form = compimag(form, imag_random_form(ex));
+    nbtest++; fpc = factorquad(form,KC,LIMC);
     if (!fpc)
     {
-      if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
+      if (DEBUGLEVEL>1) fprintferr(".");
       continue;
     }
     if (fpc > 1)
     {
-      fpd = largeprime(fpc,ex,current,0);
+      long *fpd = largeprime(fpc,ex,current,0);
+      long b1, b2, p;
+      GEN form2;
       if (!fpd)
       {
-        if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
+        if (DEBUGLEVEL>1) fprintferr(".");
         continue;
       }
-      form1 = powsubfactorbase[1][fpd[1]];
-      for (i=2; i<=lgsub; i++)
-        form1 = compimag(form1,powsubfactorbase[i][fpd[i]]);
-      pc=primeform(Disc,stoi(factorbase[fpd[-2]]),0);
-      form1=compimag(form1,pc);
-      pp = fpc << 1;
-      b1=smodis((GEN)form1[2], pp);
-      b2=smodis((GEN)form[2],  pp);
-      if (b1 != b2 && b1+b2 != pp) continue;
+      form2 = initimagform(fpd);
+      form2 = compimag(form2, imagpf(Disc, FB[fpd[-2]]));
+      p = fpc << 1;
+      b1 = smodis((GEN)form2[2], p);
+      b2 = smodis((GEN)form[2],  p);
+      if (b1 != b2 && b1+b2 != p) continue;
 
-      s++; col = mat[s];
-      if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }
-      oldfact = primfact; oldexp = exprimfact;
-      primfact = primfact1; exprimfact = exprimfact1;
-      factorisequad(form1,KC,LIMC);
-
+      col = mat[++s];
+      add_fact(col, form);
+      (void)factorquad(form2,KC,LIMC);
       if (b1==b2)
       {
-        for (i=1; i<=lgsub; i++)
-          col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];
-        col[fpd[-2]]++;
-        for (j=1; j<=primfact[0]; j++)
-        {
-          pp=primfact[j]; ep=exprimfact[j];
-          if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;
-          col[numfactorbase[pp]] -= ep;
-        }
+        for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
+        sub_fact(col, form2); col[fpd[-2]]++;
       }
       else
       {
-        for (i=1; i<=lgsub; i++)
-          col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];
-        col[fpd[-2]]--;
-        for (j=1; j<=primfact[0]; j++)
-        {
-          pp=primfact[j]; ep=exprimfact[j];
-          if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;
-          col[numfactorbase[pp]] += ep;
-        }
+        for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
+        add_fact(col, form2); col[fpd[-2]]--;
       }
-      primfact = oldfact; exprimfact = oldexp;
-    }	
+    }
     else
     {
-      s++; col = mat[s];
-      if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }
-      for (i=1; i<=lgsub; i++)
-        col[numfactorbase[subbase[i]]] = -ex[i];
+      col = mat[++s];
+      for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
+      add_fact(col, form);
     }
-    for (j=1; j<=primfact[0]; j++)
-    {
-      pp=primfact[j]; ep=exprimfact[j];
-      if (smodis((GEN)form[2], pp<<1) > pp) ep = -ep;
-      col[numfactorbase[pp]] += ep;
-    }
     col[current]--;
+    if (DEBUGLEVEL) dbg_rel(s, col);
     if (!first && fpc == 1 && col[current] == 0)
     {
       s--; for (i=1; i<=KC; i++) col[i]=0;
     }
   }
-  if (DEBUGLEVEL)
-  {
-    fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);
-    msgtimer("%s relations", first? "initial": "random");
-  }
+  if (DEBUGLEVEL) dbg_all(1, s, nbtest);
+  return C;
 }
 
 static int
-imag_be_honest(long *ex)
+imag_be_honest()
 {
-  long p,fpc, s = KC, nbtest = 0, av = avma;
-  GEN form;
+  long p, fpc, s = KC, nbtest = 0;
+  gpmem_t av = avma;
+  GEN F, ex = cgetg(lg(subFB), t_VECSMALL);
 
   while (s<KC2)
   {
-    p = factorbase[s+1];
-    if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }
-    form = imag_random_form(s+1,ex);
-    fpc = factorisequad(form,s,p-1);
+    p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
+    F = compimag(imagpf(Disc, p), imag_random_form(ex));
+    fpc = factorquad(F,s,p-1);
     if (fpc == 1) { nbtest=0; s++; }
     else
       if (++nbtest > 20) return 0;
@@ -1363,233 +1394,187 @@ imag_be_honest(long *ex)
 /*                                                                 */
 /*******************************************************************/
 
-static GEN
-real_random_form(long *ex)
+/* (1/2) log (d * 2^{e * EXP220}) */
+GEN
+get_dist(GEN e, GEN d, long prec)
 {
-  long av = avma, i;
-  GEN p1,form = NULL;
-
-  for(;;)
-  {
-    for (i=1; i<=lgsub; i++)
-    {
-      ex[i] = mymyrand()>>randshift;
-/*    if (ex[i]) KB: less efficient if I put this in. Why ??? */
-      {
-        p1 = powsubfactorbase[i][ex[i]];
-        form = form? comprealform3(form,p1): p1;
-      }
-    }
-    if (form) return form;
-    avma = av;
-  }
+  GEN t = mplog(absr(d));
+  if (signe(e)) t = addrr(t, mulir(mulsi(EXP220,e), mplog2(prec)));
+  return shiftr(t, -1);
 }
 
-static void
-real_relations(long lim, long s, long LIMC, long *ex, long **mat, GEN glog2,
-               GEN vecexpo)
+static GEN
+real_relations(long LIM, long lim, long LIMC, long **mat)
 {
-  static long nbtest;
-  long av = avma,av1,i,j,p,fpc,b1,b2,ep,current, first = (s==0);
-  long *col,*fpd,*oldfact,*oldexp,limstack;
-  long findecycle,nbrhocumule,nbrho;
-  GEN p1,p2,form,form0,form1,form2;
+  long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0, endcycle, rhoacc, rho;
+  int first = 1;
+  gpmem_t av, av1, limstack;
+  GEN C, d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
 
-  limstack=stack_lim(av,1);
-  if (first) nbtest = 0;
+  C = cgetg(LIM+1, t_VEC);
+  for (i=1; i<=LIM; i++) C[i] = lgetr(PRECREG);
+
   current = 0;
-  p1 = NULL; /* gcc -Wall */
-  while (s<lim)
+  av = avma; limstack = stack_lim(av,1);
+  s = trivial_relations(mat, C);
+  lim += s; if (lim > LIM) lim = LIM;
+NEW:
+  for(;;)
   {
+    if (s >= lim) {
+      if (lim == LIM) break;
+      lim = LIM; first = 0;
+      if (DEBUGLEVEL) dbg_all(0, s, nbtest);
+    }
+    avma = av;
     form = real_random_form(ex);
     if (!first)
     {
       current = 1+s-RELSUP;
-      p1 = redrealprimeform(Disc, factorbase[current]);
-      form = comprealform3(form,p1);
+      form = comprealform3(form, realpf(Disc, FB[current]));
     }
+    av1 = avma;
     form0 = form; form1 = NULL;
-    findecycle = nbrhocumule = 0;
-    nbrho = -1; av1 = avma;
-    while (s<lim)
+    endcycle = rhoacc = 0;
+    rho = -1;
+
+CYCLE:
+    if (endcycle) goto NEW;
+    if (low_stack(limstack, stack_lim(av,1)))
     {
-      if (low_stack(limstack, stack_lim(av,1)))
-      {
-        GEN *gptr[2];
-        int c = 1;
-	if(DEBUGMEM>1) err(warnmem,"real_relations");	
-        gptr[0]=&form; if (form1) gptr[c++]=&form1;
-        gerepilemany(av1,gptr,c);
-      }
-      if (nbrho < 0) nbrho = 0; /* first time in */
+      if(DEBUGMEM>1) err(warnmem,"real_relations");
+      gerepileall(av1, form1? 2: 1, &form, &form1);
+    }
+    if (rho < 0) rho = 0; /* first time in */
+    else
+    {
+      form = rhorealform(form); rho++;
+      rhoacc++;
+      if (first)
+        endcycle = (absi_equal((GEN)form[1],(GEN)form0[1])
+             && egalii((GEN)form[2],(GEN)form0[2])
+             && (!narrow || signe(form0[1])==signe(form[1])));
       else
       {
-        if (findecycle) break;
-        form = rhorealform(form);
-        nbrho++; nbrhocumule++;
-        if (first)
+        if (narrow)
+          { form = rhorealform(form); rho++; }
+        else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */
         {
-          if (absi_equal((GEN)form[1],(GEN)form0[1])
-               && egalii((GEN)form[2],(GEN)form0[2])
-               && (!narrow || signe(form0[1])==signe(form[1]))) findecycle=1;
+          if (absi_equal((GEN)form[1],(GEN)form0[1]) &&
+                  egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
+          form = rhorealform(form); rho++;
         }
         else
-        {
-          if (narrow)
-            { form=rhorealform(form); nbrho++; }
-          else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */
-          {
-            if (absi_equal((GEN)form[1],(GEN)form0[1]) &&
-                    egalii((GEN)form[2],(GEN)form0[2])) break;
-            form=rhorealform(form); nbrho++;
-          }
-          else
-            { setsigne(form[1],1); setsigne(form[3],-1); }
-          if (egalii((GEN)form[1],(GEN)form0[1]) &&
-              egalii((GEN)form[2],(GEN)form0[2])) break;
-        }
+          { setsigne(form[1],1); setsigne(form[3],-1); }
+        if (egalii((GEN)form[1],(GEN)form0[1]) &&
+            egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
       }
-      nbtest++; fpc = factorisequad(form,KC,LIMC);
-      if (!fpc)
+    }
+    nbtest++; fpc = factorquad(form,KC,LIMC);
+    if (!fpc)
+    {
+      if (DEBUGLEVEL>1) fprintferr(".");
+      goto CYCLE;
+    }
+    if (fpc > 1)
+    { /* look for Large Prime relation */
+      long *fpd = largeprime(fpc,ex,current,rhoacc);
+      long b1, b2, p;
+      GEN form2;
+      if (!fpd)
       {
-        if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
-        continue;
+        if (DEBUGLEVEL>1) fprintferr(".");
+        goto CYCLE;
       }
-      if (fpc > 1)
+      if (!form1) form1 = initrealform5(ex);
+      if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
+      form1 = rhoreal_pow(form1, rho);
+      rho = 0;
+
+      form2 = initrealform5(fpd);
+      if (fpd[-2]) form2 = comprealform(form2, realpf5(Disc, FB[fpd[-2]]));
+      form2 = rhoreal_pow(form2, fpd[-3]);
+      if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3]))
       {
-	fpd = largeprime(fpc,ex,current,nbrhocumule);
-        if (!fpd)
-        {
-          if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
-          continue;
-        }
-        if (!form1) form1 = initrealform5(ex);
-        if (!first)
-        {
-          p1 = redrealprimeform5(Disc, factorbase[current]);
-          form1=comprealform(form1,p1);
-        }
-        form1 = rhoreal_pow(form1, nbrho); nbrho = 0;
-        form2 = initrealform5(fpd);
-        if (fpd[-2])
-        {
-          p1 = redrealprimeform5(Disc, factorbase[fpd[-2]]);
-          form2=comprealform(form2,p1);
-        }
-        form2 = rhoreal_pow(form2, fpd[-3]);
-        if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3]))
-        {
-          setsigne(form2[1],1);
-          setsigne(form2[3],-1);
-        }
-        p = fpc << 1;
-        b1=smodis((GEN)form2[2], p);
-        b2=smodis((GEN)form1[2], p);
-        if (b1 != b2 && b1+b2 != p) continue;
+        setsigne(form2[1],1);
+        setsigne(form2[3],-1);
+      }
+      p = fpc << 1;
+      b1 = smodis((GEN)form2[2], p);
+      b2 = smodis((GEN)form1[2], p);
+      if (b1 != b2 && b1+b2 != p) goto CYCLE;
 
-        s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);
-        oldfact = primfact; oldexp = exprimfact;
-        primfact = primfact1; exprimfact = exprimfact1;
-        factorisequad(form2,KC,LIMC);
-        if (b1==b2)
-        {
-          for (i=1; i<=lgsub; i++)
-            col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];
-          for (j=1; j<=primfact[0]; j++)
-          {
-            p=primfact[j]; ep=exprimfact[j];
-            if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;
-            col[numfactorbase[p]] -= ep;
-          }
-          if (fpd[-2]) col[fpd[-2]]++; /* implies !first */
-          p1 = subii((GEN)form1[4],(GEN)form2[4]);
-          p2 = divrr((GEN)form1[5],(GEN)form2[5]);
-        }
-        else
-        {
-          for (i=1; i<=lgsub; i++)
-            col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];
-          for (j=1; j<=primfact[0]; j++)
-          {
-            p=primfact[j]; ep=exprimfact[j];
-            if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;
-            col[numfactorbase[p]] += ep;
-          }
-          if (fpd[-2]) col[fpd[-2]]--;
-          p1 = addii((GEN)form1[4],(GEN)form2[4]);
-          p2 = mulrr((GEN)form1[5],(GEN)form2[5]);
-        }
-        primfact = oldfact; exprimfact = oldexp;
+      col = mat[++s];
+      add_fact(col, form1);
+      (void)factorquad(form2,KC,LIMC);
+      if (b1==b2)
+      {
+        for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
+        sub_fact(col, form2);
+        if (fpd[-2]) col[fpd[-2]]++; /* implies !first */
+        d = get_dist(subii((GEN)form1[4],(GEN)form2[4]),
+                     divrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
       }
       else
       {
-	if (!form1) form1 = initrealform5(ex);
-        if (!first)
-        {
-          p1 = redrealprimeform5(Disc, factorbase[current]);
-          form1=comprealform(form1,p1);
-        }
-        form1 = rhoreal_pow(form1,nbrho); nbrho = 0;
-
-	s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);
-	for (i=1; i<=lgsub; i++)
-          col[numfactorbase[subbase[i]]] = -ex[i];
-        p1 = (GEN) form1[4];
-        p2 = (GEN) form1[5];
+        for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
+        add_fact(col, form2);
+        if (fpd[-2]) col[fpd[-2]]--;
+        d = get_dist(addii((GEN)form1[4],(GEN)form2[4]),
+                     mulrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
       }
-      for (j=1; j<=primfact[0]; j++)
+    }
+    else
+    { /* standard relation */
+      if (!form1) form1 = initrealform5(ex);
+      if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
+      form1 = rhoreal_pow(form1,rho);
+      rho = 0;
+
+      col = mat[++s];
+      for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
+      add_fact(col, form1);
+      d = get_dist((GEN)form1[4], (GEN)form1[5], PRECREG);
+    }
+    if (DEBUGLEVEL) fprintferr(" %ld",s);
+    affrr(d, (GEN)C[s]);
+    if (first)
+    {
+      if (s >= lim) goto NEW;
+      goto CYCLE;
+    }
+    else
+    {
+      col[current]--;
+      if (fpc == 1 && col[current] == 0)
       {
-        p=primfact[j]; ep=exprimfact[j];
-        if (smodis((GEN)form1[2], p<<1) > p) ep = -ep;
-        col[numfactorbase[p]] += ep;
+        s--; for (i=1; i<=KC; i++) col[i]=0;
       }
-      p1 = mpadd(mulir(mulsi(EXP220,p1), glog2), mplog(absr(p2)));
-      affrr(shiftr(p1,-1), (GEN)vecexpo[s]);
-      if (!first)
-      {
-        col[current]--;
-        if (fpc == 1 && col[current] == 0)
-          { s--; for (i=1; i<=KC; i++) col[i]=0; }
-        break;
-      }
     }
-    avma = av;
   }
-  if (DEBUGLEVEL)
-  {
-    fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);
-    msgtimer("%s relations", first? "initial": "random");
-  }
+  if (DEBUGLEVEL) dbg_all(1, s, nbtest);
+  return C;
 }
 
 static int
-real_be_honest(long *ex)
+real_be_honest()
 {
-  long p,fpc,s = KC,nbtest = 0,av = avma;
-  GEN p1,form,form0;
+  long p, fpc, s = KC, nbtest = 0;
+  gpmem_t av = avma;
+  GEN F,F0, ex = cgetg(lg(subFB), t_VECSMALL);
 
   while (s<KC2)
   {
-    p = factorbase[s+1];
-    if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }
-    form = real_random_form(ex);
-    p1 = redrealprimeform(Disc, p);
-    form0 = form = comprealform3(form,p1);
-    for(;;)
+    p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
+    F = comprealform3(realpf(Disc, p), real_random_form(ex));
+    for (F0 = F;;)
     {
-      fpc = factorisequad(form,s,p-1);
+      fpc = factorquad(F,s,p-1);
       if (fpc == 1) { nbtest=0; s++; break; }
-      form = rhorealform(form);
       if (++nbtest > 20) return 0;
-      if (narrow || absi_equal((GEN)form[1],(GEN)form[3]))
-	form = rhorealform(form);
-      else
-      {
-	setsigne(form[1],1);
-	setsigne(form[3],-1);
-      }
-      if (egalii((GEN)form[1],(GEN)form0[1])
-       && egalii((GEN)form[2],(GEN)form0[2])) break;
+      F = fix_signs( rhorealform(F) );
+      if (egalii((GEN)F[1],(GEN)F0[1])
+       && egalii((GEN)F[2],(GEN)F0[2])) break;
     }
     avma = av;
   }
@@ -1597,69 +1582,99 @@ real_be_honest(long *ex)
 }
 
 static GEN
-gcdrealnoer(GEN a,GEN b,long *pte)
+gcdreal(GEN a,GEN b)
 {
   long e;
   GEN k1,r;
 
+  if (!signe(a)) return mpabs(b);
+  if (!signe(b)) return mpabs(a);
+
   if (typ(a)==t_INT)
   {
     if (typ(b)==t_INT) return mppgcd(a,b);
-    k1=cgetr(lg(b)); affir(a,k1); a=k1;
+    a = itor(a, lg(b));
   }
   else if (typ(b)==t_INT)
-    { k1=cgetr(lg(a)); affir(b,k1); b=k1; }
+  {
+    b = itor(b, lg(a));
+  }
   if (expo(a)<-5) return absr(b);
   if (expo(b)<-5) return absr(a);
   a=absr(a); b=absr(b);
   while (expo(b) >= -5  && signe(b))
   {
-    k1=gcvtoi(divrr(a,b),&e);
+    k1 = gcvtoi(divrr(a,b),&e);
     if (e > 0) return NULL;
     r=subrr(a,mulir(k1,b)); a=b; b=r;
   }
-  *pte=expo(b); return absr(a);
+  return absr(a);
 }
 
-static GEN
-get_reg(GEN matc, long sreg)
+enum { RELAT, LARGE, PRECI };
+
+static int
+get_R(GEN C, long sreg, GEN z, GEN *ptR)
 {
-  long i,e,maxe;
-  GEN reg = mpabs(gcoeff(matc,1,1));
+  GEN R = gun;
+  double c;
+  long i;
 
-  e = maxe = 0;
-  for (i=2; i<=sreg; i++)
+  if (PRECREG)
   {
-    reg = gcdrealnoer(gcoeff(matc,1,i),reg,&e);
-    if (!reg) return NULL;
-    maxe = maxe? max(maxe,e): e;
+    R = mpabs((GEN)C[1]);
+    for (i=2; i<=sreg; i++)
+    {
+      R = gcdreal((GEN)C[i], R);
+      if (!R) return PRECI;
+    }
+    if (DEBUGLEVEL)
+    {
+      if (DEBUGLEVEL>7) fprintferr("R = %Z",R);
+      msgtimer("regulator");
+    }
+    if (gexpo(R) <= -3)
+    {
+      if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
+      return RELAT;
+    }
   }
+  c = gtodouble(gmul(z, R));
+  if (c < 0.75 || c > 1.5) return RELAT;
+  *ptR = R; return LARGE;
+}
+
+static int
+quad_be_honest()
+{
+  int r;
+  if (KC2 <= KC) return 1;
   if (DEBUGLEVEL)
-  {
-    if (DEBUGLEVEL>7) { fprintferr("reg = "); outerr(reg); }
-    msgtimer("regulator");
-  }
-  return reg;
+    fprintferr("be honest for primes from %ld to %ld\n", FB[KC+1],FB[KC2]);
+  r = PRECREG? real_be_honest(): imag_be_honest();
+  if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("be honest"); }
+  return r;
 }
 
 GEN
 buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec)
 {
-  long av0 = avma,av,tetpil,KCCO,KCCOPRO,i,j,s, *ex,**mat;
-  long extrarel,nrelsup,nreldep,LIMC,LIMC2,cp,nbram,nlze;
-  GEN p1,h,W,met,res,basecl,dr,c_1,pdep,C,B,extramat,extraC;
-  GEN reg,vecexpo,glog2,cst;
-  double drc,lim,LOGD;
+  gpmem_t av0 = avma, av;
+  long KCCO, i, j, s, **mat;
+  long nrelsup, nreldep, LIMC, LIMC2, cp, nlze;
+  GEN h, W, cyc, res, gen, dep, C, B, extramat, extraC;
+  GEN R, resc, Res, z;
+  double drc, lim, LOGD, LOGD2;
+  const long MAXRELSUP = 7;
 
   Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad");
-  s=mod4(Disc);
-  glog2 = vecexpo = NULL; /* gcc -Wall */
+  s = mod4(Disc);
   switch(signe(Disc))
   {
     case -1:
       if (cmpis(Disc,-4) >= 0)
       {
-        p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;
+        GEN p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;
         p1[2]=p1[3]=lgetg(1,t_VEC); return p1;
       }
       if (s==2 || s==1) err(funder2,"buchquad");
@@ -1677,156 +1692,100 @@ buchquad(GEN D, double cbach, double cbach2, long RELS
   if (!isfundamental(Disc))
     err(warner,"not a fundamental discriminant in quadclassunit");
   buch_init(); RELSUP = RELSUP0;
-  dr=cgetr(3); affir(Disc,dr); drc=fabs(rtodbl(dr)); LOGD=log(drc);
-  lim=sqrt(drc); cst = mulrr(lfunc(Disc), dbltor(lim));
+  drc = fabs(gtodouble(Disc));
+  LOGD = log(drc);
+  LOGD2 = LOGD * LOGD;
+
+  lim = sqrt(drc);
+  /* resc = sqrt(D) w / 2^r1 (2pi)^r2 ~ hR / L(chi,1) */
+  if (PRECREG) resc = dbltor(lim / 2.);
+  else         resc = dbltor(lim / PI);
   if (!PRECREG) lim /= sqrt(3.);
+  R = gun;
   cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0));
   if (cp < 13) cp = 13;
   av = avma; cbach /= 2;
+  mat = NULL;
 
-INCREASE: avma = av; cbach = check_bach(cbach,6.);
+START: avma = av; cbach = check_bach(cbach,6.);
+  if (mat) { desalloc(mat); mat = NULL; }
   nreldep = nrelsup = 0;
-  LIMC = (long)(cbach*LOGD*LOGD);
-  if (LIMC < cp) LIMC=cp;
-  LIMC2 = max(20, (long)(max(cbach,cbach2)*LOGD*LOGD));
+  LIMC = (long)(cbach*LOGD2);
+  if (LIMC < cp) { LIMC = cp; cbach = LIMC / LOGD2; }
+  LIMC2 = (long)(max(cbach,cbach2)*LOGD2);
   if (LIMC2 < LIMC) LIMC2 = LIMC;
   if (PRECREG)
   {
     PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG));
-    glog2  = glog(gdeux,PRECREG);
     sqrtD  = gsqrt(Disc,PRECREG);
     isqrtD = gfloor(sqrtD);
   }
-  factorbasequad(Disc,LIMC2,LIMC);
-  if (!KC) goto INCREASE;
 
-  vperm = cgetg(KC+1,t_VECSMALL); for (i=1; i<=KC; i++) vperm[i]=i;
-  nbram = subfactorbasequad(lim,KC);
-  if (nbram == -1) { desalloc(NULL); goto INCREASE; }
-  KCCO = KC + RELSUP;
-  if (DEBUGLEVEL) { fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO); flusherr(); }
-  powsubfact(lgsub,CBUCH+7);
-
-  mat = (long**) gpmalloc((KCCO+1)*sizeof(long*));
-  setlg(mat, KCCO+1);
-  for (i=1; i<=KCCO; i++)
-  {
-    mat[i] = (long*) gpmalloc((KC+1)*sizeof(long));
-    for (j=1; j<=KC; j++) mat[i][j]=0;
-  }
-  ex = new_chunk(lgsub+1);
+  Res = FBquad(Disc,LIMC2,LIMC);
+  if (!Res) goto START;
+  subFB = subFBquad(Disc, lim + 0.5, KC);
+  if (!subFB) goto START;
+  powsubFB = powsubFBquad(CBUCH+1);
   limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1;
   for (i=0; i<HASHT; i++) hashtab[i]=NULL;
 
-  s = lgsub+nbram+RELSUP;
-  if (PRECREG)
+  KCCO = KC + RELSUP;
+  if (DEBUGLEVEL) fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO);
+  mat = (long**)cgetalloc(NULL, KCCO+1, t_VEC);
+  for (i=1; i<=KCCO; i++)
   {
-    vecexpo=cgetg(KCCO+1,t_VEC);
-    for (i=1; i<=KCCO; i++) vecexpo[i]=lgetr(PRECREG);
-    real_relations(s,0,LIMC,ex,mat,glog2,vecexpo);
-    real_relations(KCCO,s,LIMC,ex,mat,glog2,vecexpo);
+    GEN t = cgetalloc(NULL, KC+1, t_VECSMALL);
+    for (j=1; j<=KC; j++) t[j]=0;
+    mat[i] = t;
   }
-  else
-  {
-    imag_relations(s,0,LIMC,ex,mat);
-    imag_relations(KCCO,s,LIMC,ex,mat);
-  }
-  if (KC2 > KC)
-  {
-    if (DEBUGLEVEL)
-      fprintferr("be honest for primes from %ld to %ld\n",
-                  factorbase[KC+1],factorbase[KC2]);
-    s = PRECREG? real_be_honest(ex): imag_be_honest(ex);
-    if (DEBUGLEVEL)
-    {
-      fprintferr("\n");
-      msgtimer("be honest");
-    }
-    if (!s) { desalloc(mat); goto INCREASE; }
-  }
-  C=cgetg(KCCO+1,t_MAT);
-  if (PRECREG)
-  {
-    for (i=1; i<=KCCO; i++)
-    {
-      C[i]=lgetg(2,t_COL); coeff(C,1,i)=vecexpo[i];
-    }
-    if (DEBUGLEVEL>7) { fprintferr("C: "); outerr(C); flusherr(); }
-  }
-  else
-    for (i=1; i<=KCCO; i++) C[i]=lgetg(1,t_COL);
-  W = hnfspec(mat,vperm,&pdep,&B,&C,lgsub);
-  nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
 
-  KCCOPRO=KCCO;
+  s = lg(subFB)-1 + RELSUP;
+  C = PRECREG? real_relations(KCCO,s,LIMC,mat)
+             : imag_relations(KCCO,s,LIMC,mat);
+  W = hnfspec(mat,vperm,&dep,&B,&C,lg(subFB)-1);
+  nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
+
   if (nlze)
   {
-EXTRAREL:
-    s = PRECREG? 2: 1; extrarel=nlze+2;
-    extraC=cgetg(extrarel+1,t_MAT);
-    for (i=1; i<=extrarel; i++) extraC[i]=lgetg(s,t_COL);
-    extramat = extra_relations(LIMC,ex,nlze,extraC);
-    if (!extramat) { desalloc(mat); goto INCREASE; }
-    W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);
-    nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
-    KCCOPRO += extrarel;
+MORE:
+    extramat = extra_relations(LIMC,nlze, &extraC);
+    if (!extramat) { goto START; }
+    W = hnfadd(W,vperm,&dep,&B,&C, extramat,extraC);
+    nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
+    KCCO += lg(extramat)-1;
     if (nlze)
     {
-      if (++nreldep > 5) { desalloc(mat); goto INCREASE; }
-      goto EXTRAREL;
+      if (++nreldep > 5) goto START;
+      goto MORE;
     }
   }
-  /* tentative class number */
-  h=gun; for (i=1; i<lg(W); i++) h=mulii(h,gcoeff(W,i,i));
-  if (PRECREG)
-  {
-    /* tentative regulator */
-    reg = get_reg(C, KCCOPRO - (lg(B)-1) - (lg(W)-1));
-    if (!reg)
-    {
-      desalloc(mat);
-      prec = (PRECREG<<1)-2; goto INCREASE;
-    }
-    if (gexpo(reg)<=-3)
-    {
-      if (++nrelsup <= 7)
-      {
-        if (DEBUGLEVEL) { fprintferr("regulateur nul\n"); flusherr(); }
-        nlze=min(KC,nrelsup); goto EXTRAREL;
-      }
-      desalloc(mat); goto INCREASE;
-    }
-    c_1 = divrr(gmul2n(gmul(h,reg),1), cst);
-  }
-  else
-  {
-    reg = gun;
-    c_1 = divrr(gmul(h,mppi(DEFAULTPREC)), cst);
-  }
 
-  if (gcmpgs(gmul2n(c_1,2),3)<0) { c_1=stoi(10); nrelsup=7; }
-  if (gcmpgs(gmul2n(c_1,1),3)>0)
+  h = dethnf_i(W);
+  if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h);
+
+  z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */
+  switch(get_R(C, KCCO - (lg(B)-1) - (lg(W)-1), divir(h,z), &R))
   {
-    nrelsup++;
-    if (nrelsup<=7)
-    {
-      if (DEBUGLEVEL)
-        { fprintferr("***** check = %f\n\n",gtodouble(c_1)); flusherr(); }
-      nlze=min(KC,nrelsup); goto EXTRAREL;
-    }
-    if (cbach < 5.99) { desalloc(mat); goto INCREASE; }
-    err(warner,"suspicious check. Suggest increasing extra relations.");
+    case PRECI:
+      prec = (PRECREG<<1)-2;
+      goto START;
+
+    case RELAT:
+      if (++nrelsup <= MAXRELSUP) { nlze = min(KC, nrelsup); goto MORE; }
+      goto START;
   }
-  basecl = get_clgp(Disc,W,&met,PRECREG);
-  s = lg(basecl); desalloc(mat); tetpil=avma;
+  /* DONE */
+  if (!quad_be_honest()) goto START;
 
+  gen = get_clgp(Disc,W,&cyc,PRECREG);
+  desalloc(mat);
+
   res=cgetg(6,t_VEC);
-  res[1]=lcopy(h); p1=cgetg(s,t_VEC);
-  for (i=1; i<s; i++) p1[i] = (long)icopy(gcoeff(met,i,i));
-  res[2]=(long)p1;
-  res[3]=lcopy(basecl);
-  res[4]=lcopy(reg);
-  res[5]=lcopy(c_1); return gerepile(av0,tetpil,res);
+  res[1]=(long)h;
+  res[2]=(long)cyc;
+  res[3]=(long)gen;
+  res[4]=(long)R;
+  res[5]=ldiv(mpmul(R,h), z); return gerepilecopy(av0,res);
 }
 
 GEN