===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/buch3.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/buch3.c	2001/10/02 11:17:03	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/buch3.c	2002/09/11 07:26:50	1.2
@@ -1,4 +1,4 @@
-/* $Id: buch3.c,v 1.1 2001/10/02 11:17:03 noro Exp $
+/* $Id: buch3.c,v 1.2 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -21,33 +21,36 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"
 #include "parinf.h"
 
-extern GEN concatsp3(GEN x, GEN y, GEN z);
+extern void testprimes(GEN bnf, long bound);
+extern GEN Fp_PHlog(GEN a, GEN g, GEN p, GEN ord);
+extern GEN FqX_factor(GEN x, GEN T, GEN p);
+extern GEN anti_unif_mod_f(GEN nf, GEN pr, GEN sqf);
+extern GEN arch_mul(GEN x, GEN y);
 extern GEN check_and_build_cycgen(GEN bnf);
+extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
+extern GEN detcyc(GEN cyc);
+extern GEN famat_reduce(GEN fa);
+extern GEN famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id);
+extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);
 extern GEN gmul_mat_smallvec(GEN x, GEN y);
+extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
+extern GEN idealprodprime(GEN nf, GEN L);
 extern GEN ideleaddone_aux(GEN nf,GEN x,GEN ideal);
+extern GEN isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag);
 extern GEN logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid);
-extern GEN vconcat(GEN Q1, GEN Q2);
-extern void minim_alloc(long n, double ***q, GEN *x, double **y,  double **z, double **v);
+extern GEN make_integral(GEN nf, GEN L0, GEN f, GEN *listpr, GEN *ptd1);
+extern GEN sqred1_from_QR(GEN x, long prec);
+extern GEN subgroupcondlist(GEN cyc, GEN bound, GEN listKer);
+extern GEN to_Fp_simple(GEN nf, GEN x, GEN ffproj);
 extern GEN to_famat_all(GEN x, GEN y);
 extern GEN trivfact(void);
-extern GEN arch_mul(GEN x, GEN y);
-extern GEN isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag);
-extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
-extern GEN ideallllred_elt(GEN nf, GEN I);
+extern GEN unif_mod_f(GEN nf, GEN pr, GEN sqf);
+extern GEN vconcat(GEN Q1, GEN Q2);
+extern long FqX_is_squarefree(GEN P, GEN T, GEN p);
+extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
+extern void minim_alloc(long n, double ***q, GEN *x, double **y,  double **z, double **v);
+extern void rowselect_p(GEN A, GEN B, GEN p, long init);
 
-/* U W V = D, Ui = U^(-1) */
-GEN
-compute_class_number(GEN W, GEN *D,GEN *Ui,GEN *V)
-{
-  GEN S = smith2(W);
-
-  *Ui= ginv((GEN)S[1]);
-  if (V) *V = (GEN)S[2];
-  *D = (GEN)S[3];
-  if (DEBUGLEVEL>=4) msgtimer("smith/class group");
-  return dethnf_i(*D);
-}
-
 /* FIXME: obsolete, see zarchstar (which is much slower unfortunately). */
 static GEN
 get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen, long ngen, long rankmax)
@@ -63,7 +66,7 @@ get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen, 
     limr = (limr-1)>>1;
     for (k=rr;  k<=limr; k++)
     {
-      long av1=avma;
+      gpmem_t av1=avma;
       alpha = gzero;
       for (kk=k,i=1; i<=N; i++,kk/=rr)
       {
@@ -75,13 +78,12 @@ get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen, 
 	vecsign[i] = (gsigne(gsubst(alpha,va,(GEN)rac[i])) > 0)? (long)_0
                                                                : (long)_1;
       v1 = concatsp(v, vecsign);
-      if (rank(v1) == rankinit) avma=av1;
-      else
-      {
-	v=v1; rankinit++; ngen++;
-        gen[ngen] = (long) alpha;
-	if (rankinit == rankmax) return ginv(v); /* full rank */
-      }
+      if (rank(v1) == rankinit) { avma = av1; continue; }
+
+      v=v1; rankinit++; ngen++;
+      gen[ngen] = (long) alpha;
+      if (rankinit == rankmax) return ginv(v); /* full rank */
+      vecsign = cgetg(rankmax+1,t_COL);
     }
   }
 }
@@ -91,9 +93,9 @@ GEN
 buchnarrow(GEN bnf)
 {
   GEN nf,_0,_1,cyc,gen,v,matsign,arch,cycgen,logs;
-  GEN dataunit,p1,p2,h,clh,basecl,met,u1;
+  GEN dataunit,p1,p2,h,basecl,met,u1;
   long r1,R,i,j,ngen,sizeh,t,lo,c;
-  ulong av = avma;
+  gpmem_t av = avma;
 
   bnf = checkbnf(bnf);
   nf = checknf(bnf); r1 = nf_get_r1(nf);
@@ -134,12 +136,12 @@ buchnarrow(GEN bnf)
     vconcat(diagonal(cyc), logs),
     vconcat(zeromat(ngen, r1-t), gscalmat(gdeux,r1-t))
   );
-  clh = compute_class_number(h,&met,&u1,NULL);
-  lo = lg(met)-1;
-  for (c=1; c<=lo; c++)
-    if (gcmp1(gcoeff(met,c,c))) break;
+ 
+  lo = lg(h)-1;
+  met = smithrel(h,NULL,&u1);
+  c = lg(met);
+  if (DEBUGLEVEL>3) msgtimer("smith/class group");
 
-  u1 = reducemodmatrix(u1,h);
   basecl = cgetg(c,t_VEC);
   for (j=1; j<c; j++)
   {
@@ -152,32 +154,80 @@ buchnarrow(GEN bnf)
       if (signe(p1))
       {
 	p2 = idealmul(nf,p2, idealpow(nf,(GEN)gen[i],p1));
-        p1 = content(p2); if (!gcmp1(p1)) p2 = gdiv(p2,p1);
+        p2 = primpart(p2);
       }
     }
     basecl[j] = (long)p2;
   }
   v = cgetg(4,t_VEC);
-  v[1] = lcopy(clh); setlg(met,c);
-  v[2] = (long)mattodiagonal(met);
-  v[3] = lcopy(basecl); return gerepileupto(av, v);
+  v[1] = (long)dethnf_i(h);
+  v[2] = (long)met;
+  v[3] = (long)basecl; return gerepilecopy(av, v);
 }
 
-/* given two coprime ideals x (integral) and id, compute alpha in x,
- * alpha = 1 mod (id), with x/alpha nearly reduced.
- */
+/********************************************************************/
+/**                                                                **/
+/**                  REDUCTION MOD IDELE                           **/
+/**                                                                **/
+/********************************************************************/
+
 static GEN
-findalpha(GEN nf,GEN x,GEN id)
+compute_fact(GEN nf, GEN u1, GEN gen)
 {
-  GEN p1, y;
-  GEN alp = idealaddtoone_i(nf,x,id);
+  GEN G, basecl;
+  long prec,i,j, l = lg(u1), h = lg(u1[1]); /* l > 1 */
 
-  y = ideallllred_elt(nf, idealmullll(nf,x,id));
-  p1 = ground(element_div(nf,alp,y));
-  alp = gsub(alp, element_mul(nf,p1,y));
-  return gcmp0(alp)? y: alp;
+  basecl = cgetg(l,t_VEC);
+  prec = nfgetprec(nf);
+  G = cgetg(3,t_VEC);
+  G[2] = lgetg(1,t_MAT);
+
+  for (j=1; j<l; j++)
+  {
+    GEN g,e, z = NULL;
+    for (i=1; i<h; i++)
+    {
+      e = gcoeff(u1,i,j); if (!signe(e)) continue;
+
+      g = (GEN)gen[i];
+      if (typ(g) != t_MAT)
+      {
+        if (z) 
+          z[2] = (long)arch_mul((GEN)z[2], to_famat_all(g, e));
+        else
+        {
+          z = cgetg(3,t_VEC);
+          z[1] = 0;
+          z[2] = (long)to_famat_all(g, e);
+        }
+        continue;
+      }
+
+      G[1] = (long)g;
+      g = idealpowred(nf,G,e,prec);
+      z = z? idealmulred(nf,z,g,prec): g;
+    }
+    z[2] = (long)famat_reduce((GEN)z[2]);
+    basecl[j] = (long)z;
+  }
+  return basecl;
 }
 
+/* given two coprime integral ideals x and f (f HNF), compute "small"
+ * non-zero a in x, such that a = 1 mod (f). GTM 193: Algo 4.3.3 */
+static GEN
+redideal(GEN nf,GEN x,GEN f)
+{
+  GEN q, y, b;
+
+  if (gcmp1(gcoeff(f,1,1))) return idealred_elt(nf, x); /* f = 1 */
+
+  b = idealaddtoone_i(nf,x,f); /* a = b mod (x f) */
+  y = idealred_elt(nf, idealmullll(nf,x,f));
+  q = ground(element_div(nf,b,y));
+  return gsub(b, element_mul(nf,q,y)); /* != 0 since = 1 mod f */
+}
+
 static int
 too_big(GEN nf, GEN bet)
 {
@@ -191,76 +241,189 @@ too_big(GEN nf, GEN bet)
   return 0; /* not reached */
 }
 
+/* GTM 193: Algo 4.3.4. Reduce x mod idele */
 static GEN
-idealmodidele(GEN nf, GEN x, GEN ideal, GEN sarch, GEN arch)
+_idealmodidele(GEN nf, GEN x, GEN idele, GEN sarch)
 {
-  long av = avma,i,l;
-  GEN p1,p2,alp,bet,b;
+  gpmem_t av = avma;
+  GEN a,A,D,G, f = (GEN)idele[1];
 
-  nf=checknf(nf); alp=findalpha(nf,x,ideal);
-  p1=idealdiv(nf,alp,x);
-  bet = element_div(nf,findalpha(nf,p1,ideal),alp);
-  if (too_big(nf,bet) > 0) { avma=av; return x; }
-  p1=(GEN)sarch[2]; l=lg(p1);
-  if (l > 1)
-  {
-    b=bet; p2=lift_intern(gmul((GEN)sarch[3],zsigne(nf,bet,arch)));
-    for (i=1; i<l; i++)
-    if (signe(p2[i])) bet = element_mul(nf,bet,(GEN)p1[i]);
-    if (b != bet && too_big(nf,bet) > 0) { avma=av; return x; }
-  }
-  return idealmul(nf,bet,x);
+  G = redideal(nf, x, f);
+  D = redideal(nf, idealdiv(nf,G,x), f);
+  A = element_div(nf,D,G);
+  if (too_big(nf,A) > 0) { avma = av; return x; }
+  a = set_sign_mod_idele(nf, NULL, A, idele, sarch);
+  if (a != A && too_big(nf,A) > 0) { avma = av; return x; }
+  return idealmul(nf, a, x);
 }
 
-static GEN
-idealmulmodidele(GEN nf,GEN x,GEN y, GEN ideal,GEN sarch,GEN arch)
+GEN
+idealmodidele(GEN bnr, GEN x)
 {
-  return idealmodidele(nf,idealmul(nf,x,y),ideal,sarch,arch);
+  GEN bid = (GEN)bnr[2], fa2 = (GEN)bid[4];
+  GEN idele = (GEN)bid[1];
+  GEN sarch = (GEN)fa2[lg(fa2)-1];
+  return _idealmodidele(checknf(bnr), x, idele, sarch);
 }
 
-/* assume n > 0 */
-/* FIXME: should compute x^n = a I using idealred, then reduce a mod idele */
+/* v_pr(L0 * cx). tau = pr[5] or (more efficient) mult. table for pr[5] */
+static long
+fast_val(GEN nf,GEN L0,GEN cx,GEN pr,GEN tau)
+{
+  GEN p = (GEN)pr[1];
+  long v = int_elt_val(nf,L0,p,tau,NULL);
+  if (cx)
+  {
+    long w = ggval(cx, p);
+    if (w) v += w * itos((GEN)pr[3]);
+  }
+  return v;
+}
+
 static GEN
-idealpowmodidele(GEN nf,GEN x,GEN n, GEN ideal,GEN sarch,GEN arch)
+compute_raygen(GEN nf, GEN u1, GEN gen, GEN bid)
 {
-  long i,m,av=avma;
-  GEN y;
-  ulong j;
+  GEN f, fZ, basecl, module, fa, fa2, *listpr, *listep, *vecinvpi, *vectau;
+  GEN pr, t, sqf, EX, sarch, cyc;
+  long i,j,l,lp;
 
-  if (cmpis(n, 16) < 0)
+  /* basecl = generators in factored form */
+  basecl = compute_fact(nf,u1,gen);
+
+  module = (GEN)bid[1];
+  cyc = gmael(bid,2,2); EX = (GEN)cyc[1]; /* exponent of (O/f)^* */
+  f   = (GEN)module[1]; fZ = gcoeff(f,1,1);
+  fa  = (GEN)bid[3];
+  fa2 = (GEN)bid[4]; sarch = (GEN)fa2[lg(fa2)-1];
+  listpr = (GEN*)fa[1];
+  listep = (GEN*)fa[2];
+
+  lp = lg(listpr);
+  /* sqf = squarefree kernel of f */
+  sqf = lp <= 2? NULL: idealprodprime(nf, (GEN)listpr);
+
+  vecinvpi = (GEN*)cgetg(lp, t_VEC);
+  vectau   = (GEN*)cgetg(lp, t_VEC);
+  for (i=1; i<lp; i++) 
   {
-    if (gcmp1(n)) return x;
-    x = idealpow(nf,x,n);
-    x = idealmodidele(nf,x,ideal,sarch,arch);
-    return gerepileupto(av,x);
+    pr = listpr[i];
+    vecinvpi[i] = NULL; /* to be computed if needed */
+    vectau[i] = eltmul_get_table(nf, (GEN)pr[5]);
   }
 
-  i = lgefint(n)-1; m=n[i]; j=HIGHBIT;
-  while ((m&j)==0) j>>=1;
-  y = x;
-  for (j>>=1; j; j>>=1)
+  l = lg(basecl);
+  for (i=1; i<l; i++)
   {
-    y = idealmul(nf,y,y);
-    if (m&j) y = idealmul(nf,y,x);
-    y = idealmodidele(nf,y,ideal,sarch,arch);
-  }
-  for (i--; i>=2; i--)
-    for (m=n[i],j=HIGHBIT; j; j>>=1)
+    GEN d, invpi, mulI, G, I, A, e, L, newL;
+    long la, v, k;
+    /* G = [I, A=famat(L,e)] is a generator, I integral */
+    G = (GEN)basecl[i];
+    I = (GEN)G[1];
+    A = (GEN)G[2];
+      L = (GEN)A[1];
+      e = (GEN)A[2];
+    /* if no reduction took place in compute_fact, everybody is still coprime
+     * to f + no denominators */
+    if (!I)
     {
-      y = idealmul(nf,y,y);
-      if (m&j) y = idealmul(nf,y,x);
-      y = idealmodidele(nf,y,ideal,sarch,arch);
+      basecl[i] = (long)famat_to_nf_modidele(nf, L, e, bid);
+      continue;
     }
-  return gerepileupto(av,y);
+    if (lg(A) == 1)
+    {
+      basecl[i] = (long)I;
+      continue;
+    }
+
+    /* compute mulI so that mulI * I coprime to f
+     * FIXME: use idealcoprime ??? (Less efficient. Fix idealcoprime!) */
+    mulI = NULL;
+    for (j=1; j<lp; j++)
+    {
+      invpi = vecinvpi[j];
+      pr    = listpr[j]; 
+      v = idealval(nf, I, pr);
+      if (v) {
+        if (!invpi) invpi = vecinvpi[j] = anti_unif_mod_f(nf, pr, sqf);
+        t = element_pow(nf,invpi,stoi(v));
+        mulI = mulI? element_mul(nf, mulI, t): t;
+      }
+    }
+
+    /* make all components of L coprime to f. 
+     * Assuming (L^e * I, f) = 1, then newL^e * mulI = L^e */
+    la = lg(e); newL = cgetg(la, t_VEC);
+    for (k=1; k<la; k++)
+    {
+      GEN L0, cx, LL = (GEN)L[k];
+      if (typ(LL) != t_COL) LL = algtobasis(nf, LL);
+
+      L0 = Q_primitive_part(LL, &cx); /* LL = L0*cx (faster element_val) */
+      for (j=1; j<lp; j++)
+      {
+        invpi = vecinvpi[j];
+        pr  = listpr[j];
+        v = fast_val(nf, L0,cx, pr,vectau[j]); /* = val_pr(LL) */
+        if (v) {
+          if (!invpi) invpi = vecinvpi[j] = anti_unif_mod_f(nf, pr, sqf);
+          t = element_pow(nf,invpi,stoi(v));
+          LL = element_mul(nf, LL, t);
+        }
+      }
+      LL = make_integral(nf, LL, f, listpr, NULL);
+      newL[k] = (long)FpV_red(LL, fZ);
+    }
+
+    /* G in nf, = L^e mod f */
+    G = famat_to_nf_modideal_coprime(nf, newL, gmod(e,EX), f);
+    d = NULL;
+    if (mulI)
+    {
+      GEN mod;
+
+      /* mulI integral, d | mulI * I,  sqf(dO_K) | f */
+      mulI = make_integral(nf,mulI,f,listpr, &d);
+
+      /* reduce mulI mod (f * (mulI,d)) */
+      if (!d) mod = f;
+      else
+      {
+        mod = hnfmodid(eltmul_get_table(nf,mulI), d); /* = (mulI,d) */
+        mod = idealmul(nf,f, mod);
+      }
+      mulI = colreducemodHNF(mulI, mod, NULL);
+      G = element_muli(nf, G, mulI);
+
+      /* L^e * I = (G/d * I) mod f */
+      G = colreducemodHNF(G, mod, NULL);
+    }
+
+    G = set_sign_mod_idele(nf,A,G,module,sarch);
+    I = idealmul(nf,I,G);
+    if (d) I = Q_div_to_int(I,d);
+    /* more or less useless, but cheap at this point */
+    I = _idealmodidele(nf,I,module,sarch);
+    basecl[i] = (long)I;
+  }
+  return basecl;
 }
 
+/********************************************************************/
+/**                                                                **/
+/**                   INIT RAY CLASS GROUP                         **/
+/**                                                                **/
+/********************************************************************/
+
 static GEN
 buchrayall(GEN bnf,GEN module,long flag)
 {
-  GEN nf,cyc,gen,genplus,fa2,sarch,hmatu,u,clg,logs;
-  GEN dataunit,p1,p2,h,clh,genray,met,u1,u2,u1old,cycgen;
+  GEN nf,cyc,gen,genplus,hmatu,u,clg,logs;
+  GEN dataunit,p1,h,met,u1,u2,U,cycgen;
   GEN racunit,bigres,bid,cycbid,genbid,x,y,funits,hmat,vecel;
-  long RU,Ri,i,j,ngen,lh,lo,c,av=avma;
+  long RU, Ri, j, ngen, lh;
+  const int add_gen = flag & nf_GEN;
+  const int do_init = flag & nf_INIT;
+  gpmem_t av=avma;
 
   bnf = checkbnf(bnf); nf = checknf(bnf);
   funits = check_units(bnf, "buchrayall"); RU = lg(funits);
@@ -275,7 +438,7 @@ buchrayall(GEN bnf,GEN module,long flag)
   Ri = lg(cycbid)-1; lh = ngen+Ri;
 
   x = idealhermite(nf,module);
-  if (Ri || flag & (nf_INIT|nf_GEN))
+  if (Ri || add_gen || do_init)
   {
     vecel = cgetg(ngen+1,t_VEC);
     for (j=1; j<=ngen; j++)
@@ -285,34 +448,32 @@ buchrayall(GEN bnf,GEN module,long flag)
       vecel[j]=(long)p1;
     }
   }
-  if (flag & nf_GEN)
+  if (add_gen)
   {
     genplus = cgetg(lh+1,t_VEC);
     for (j=1; j<=ngen; j++)
-      genplus[j] = (long) idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);
+      genplus[j] = (long)idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);
     for (  ; j<=lh; j++)
       genplus[j] = genbid[j-ngen];
   }
   if (!Ri)
   {
-    if (!(flag & nf_GEN)) clg = cgetg(3,t_VEC);
-    else
-      { clg = cgetg(4,t_VEC); clg[3] = (long)genplus; }
+    clg = cgetg(add_gen? 4: 3,t_VEC);
+    if (add_gen) clg[3] = (long)genplus;
     clg[1] = mael(bigres,1,1);
     clg[2] = (long)cyc;
-    if (!(flag & nf_INIT)) return gerepilecopy(av,clg);
+    if (!do_init) return gerepilecopy(av,clg);
     y = cgetg(7,t_VEC);
-    y[1] = lcopy(bnf);
-    y[2] = lcopy(bid);
-    y[3] = lcopy(vecel);
+    y[1] = (long)bnf;
+    y[2] = (long)bid;
+    y[3] = (long)vecel;
     y[4] = (long)idmat(ngen);
-    y[5] = lcopy(clg); u = cgetg(3,t_VEC);
+    y[5] = (long)clg; u = cgetg(3,t_VEC);
     y[6] = (long)u;
       u[1] = lgetg(1,t_MAT);
       u[2] = (long)idmat(RU);
-    return gerepileupto(av,y);
+    return gerepilecopy(av,y);
   }
-  fa2 = (GEN)bid[4]; sarch = (GEN)fa2[lg(fa2)-1];
 
   cycgen = check_and_build_cycgen(bnf);
   dataunit = cgetg(RU+1,t_MAT); racunit = gmael(bigres,4,2);
@@ -340,61 +501,29 @@ buchrayall(GEN bnf,GEN module,long flag)
     vconcat(diagonal(cyc), gneg_i(logs)),
     vconcat(zeromat(ngen, Ri), hmat)
   );
-  clh = compute_class_number(h,&met,&u1,NULL);
-  u1old = u1; lo = lg(met)-1;
-  for (c=1; c<=lo; c++)
-    if (gcmp1(gcoeff(met,c,c))) break;
+  met = smithrel(hnf(h), &U, add_gen? &u1: NULL);
+  clg = cgetg(add_gen? 4: 3, t_VEC);
+  clg[1] = (long)detcyc(met);
+  clg[2] = (long)met;
+  if (add_gen) clg[3] = (long)compute_raygen(nf,u1,genplus,bid);
+  if (!do_init) return gerepilecopy(av, clg);
 
-  if (flag & nf_GEN)
-  {
-    GEN Id = idmat(degpol(nf[1])), arch = (GEN)module[2];
-    u1 = reducemodmatrix(u1,h);
-    genray = cgetg(c,t_VEC);
-    for (j=1; j<c; j++)
-    {
-      GEN *op, minus = Id, plus = Id;
-      long av1 = avma, s;
-      for (i=1; i<=lo; i++)
-      {
-	p1 = gcoeff(u1,i,j);
-        if (!(s = signe(p1))) continue;
-
-        if (s > 0) op = &plus; else { op = &minus; p1 = negi(p1); }
-        p1 = idealpowmodidele(nf,(GEN)genplus[i],p1,x,sarch,arch);
-        *op = *op==Id? p1
-                     : idealmulmodidele(nf,*op,p1,x,sarch,arch);
-      }
-      if (minus == Id) p1 = plus;
-      else
-      {
-        p2 = ideleaddone_aux(nf,minus,module);
-        p1 = idealdivexact(nf,idealmul(nf,p2,plus),minus);
-        p1 = idealmodidele(nf,p1,x,sarch,arch);
-      }
-      genray[j]=lpileupto(av1,p1);
-    }
-    clg = cgetg(4,t_VEC); clg[3] = lcopy(genray);
-  } else clg = cgetg(3,t_VEC);
-  clg[1] = licopy(clh); setlg(met,c);
-  clg[2] = (long)mattodiagonal(met);
-  if (!(flag & nf_INIT)) return gerepileupto(av,clg);
-
   u2 = cgetg(Ri+1,t_MAT);
   u1 = cgetg(RU+1,t_MAT); u = (GEN)hmatu[2];
   for (j=1; j<=RU; j++) { u1[j]=u[j]; setlg(u[j],RU+1); }
   u += RU;
   for (j=1; j<=Ri; j++) { u2[j]=u[j]; setlg(u[j],RU+1); }
-  p1 = lllint(u1); p2 = ginv(hmat);
+ 
   y = cgetg(7,t_VEC);
-  y[1] = lcopy(bnf);
-  y[2] = lcopy(bid);
-  y[3] = lcopy(vecel);
-  y[4] = linv(u1old);
-  y[5] = lcopy(clg); u = cgetg(3,t_VEC);
+  y[1] = (long)bnf;
+  y[2] = (long)bid;
+  y[3] = (long)vecel;
+  y[4] = (long)U;
+  y[5] = (long)clg; u = cgetg(3,t_VEC);
   y[6] = (long)u;
-    u[1] = lmul(u2,p2);
-    u[2] = lmul(u1,p1);
-  return gerepileupto(av,y);
+    u[1] = lmul(u2, ginv(hmat));
+    u[2] = (long)lllint_ip(u1,100);
+  return gerepilecopy(av,y);
 }
 
 GEN
@@ -445,7 +574,7 @@ rayclassno(GEN bnf,GEN ideal)
 {
   GEN nf,h,dataunit,racunit,bigres,bid,cycbid,funits,H;
   long RU,i;
-  ulong av = avma;
+  gpmem_t av = avma;
 
   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
   bigres = (GEN)bnf[8]; h = gmael(bigres,1,1); /* class number */
@@ -470,7 +599,7 @@ quick_isprincipalgen(GEN bnf, GEN x)
   GEN z = cgetg(3, t_VEC), gen = gmael3(bnf,8,1,3);
   GEN idep, ep = isprincipal(bnf,x);
   /* x \prod g[i]^(-ep[i]) = factorisation of principal ideal */
-  idep = isprincipalfact(bnf, gen, gneg(ep), x, nf_GENMAT | nf_GEN);
+  idep = isprincipalfact(bnf, gen, gneg(ep), x, nf_GENMAT|nf_FORCE);
   z[1] = (long)ep;
   z[2] = idep[2]; return z;
 }
@@ -478,8 +607,9 @@ quick_isprincipalgen(GEN bnf, GEN x)
 GEN
 isprincipalrayall(GEN bnr, GEN x, long flag)
 {
-  long av=avma,i,j,c;
-  GEN bnf,nf,bid,matu,vecel,ep,p1,beta,idep,y,rayclass;
+  long i, j, c;
+  gpmem_t av=avma;
+  GEN bnf,nf,bid,matu,vecel,ep,p1,beta,idep,ex,rayclass;
   GEN divray,genray,alpha,alphaall,racunit,res,funit;
 
   checkbnr(bnr);
@@ -488,6 +618,9 @@ isprincipalrayall(GEN bnr, GEN x, long flag)
   vecel=(GEN)bnr[3];
   matu =(GEN)bnr[4];
   rayclass=(GEN)bnr[5];
+  divray = (GEN)rayclass[2]; c = lg(divray)-1;
+  ex = cgetg(c+1,t_COL);
+  if (c == 0 && !(flag & nf_GEN)) return ex;
 
   if (typ(x) == t_VEC && lg(x) == 3)
   { idep = (GEN)x[2]; x = (GEN)x[1]; }  /* precomputed */
@@ -495,16 +628,14 @@ isprincipalrayall(GEN bnr, GEN x, long flag)
     idep = quick_isprincipalgen(bnf, x);
   ep  = (GEN)idep[1];
   beta= (GEN)idep[2];
-  c = lg(ep);
-  for (i=1; i<c; i++) /* modify beta as if gen -> vecel.gen (coprime to bid) */
+  j = lg(ep);
+  for (i=1; i<j; i++) /* modify beta as if gen -> vecel.gen (coprime to bid) */
     if (typ(vecel[i]) != t_INT && signe(ep[i])) /* <==> != 1 */
       beta = arch_mul(to_famat_all((GEN)vecel[i], negi((GEN)ep[i])), beta);
   p1 = gmul(matu, concatsp(ep, zideallog(nf,beta,bid)));
-  divray = (GEN)rayclass[2]; c = lg(divray);
-  y = cgetg(c,t_COL);
-  for (i=1; i<c; i++)
-    y[i] = lmodii((GEN)p1[i],(GEN)divray[i]);
-  if (!(flag & nf_GEN)) return gerepileupto(av, y);
+  for (i=1; i<=c; i++)
+    ex[i] = lmodii((GEN)p1[i],(GEN)divray[i]);
+  if (!(flag & nf_GEN)) return gerepileupto(av, ex);
 
   /* compute generator */
   if (lg(rayclass)<=3)
@@ -512,7 +643,7 @@ isprincipalrayall(GEN bnr, GEN x, long flag)
 
   genray = (GEN)rayclass[3];
   /* TODO: should be using nf_GENMAT and function should return a famat */
-  alphaall = isprincipalfact(bnf, genray, gneg(y), x, nf_GEN | nf_FORCE);
+  alphaall = isprincipalfact(bnf, genray, gneg(ex), x, nf_GEN | nf_FORCE);
   if (!gcmp0((GEN)alphaall[1])) err(bugparier,"isprincipalray (bug1)");
 
   res = (GEN)bnf[8];
@@ -533,7 +664,7 @@ isprincipalrayall(GEN bnr, GEN x, long flag)
     alpha = gdiv(alpha,p1);
   }
   p1 = cgetg(4,t_VEC);
-  p1[1] = lcopy(y);
+  p1[1] = lcopy(ex);
   p1[2] = (long)algtobasis(nf,alpha);
   p1[3] = lcopy((GEN)alphaall[3]);
   return gerepileupto(av, p1);
@@ -542,7 +673,7 @@ isprincipalrayall(GEN bnr, GEN x, long flag)
 GEN
 isprincipalray(GEN bnr, GEN x)
 {
-  return isprincipalrayall(bnr,x,nf_REGULAR);
+  return isprincipalrayall(bnr,x,0);
 }
 
 GEN
@@ -551,10 +682,11 @@ isprincipalraygen(GEN bnr, GEN x)
   return isprincipalrayall(bnr,x,nf_GEN);
 }
 
+/* N! / N^N * (4/pi)^r2 * sqrt(|D|) */
 GEN
 minkowski_bound(GEN D, long N, long r2, long prec)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN p1;
   p1 = gdiv(mpfactr(N,prec), gpowgs(stoi(N),N));
   p1 = gmul(p1, gpowgs(gdivsg(4,mppi(prec)), r2));
@@ -566,13 +698,13 @@ minkowski_bound(GEN D, long N, long r2, long prec)
 static long
 zimmertbound(long N,long R2,GEN DK)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN w;
 
   if (N < 2) return 1;
   if (N < 21)
   {
-    static double c[21][11] = {
+    static double c[19][11] = {
 {/*2*/  0.6931,     0.45158},
 {/*3*/  1.71733859, 1.37420604},
 {/*4*/  2.91799837, 2.50091538, 2.11943331},
@@ -619,68 +751,12 @@ zimmertbound(long N,long R2,GEN DK)
   avma = av; return itos(w);
 }
 
-/* all primes up to Minkowski bound factor on factorbase ? */
-static void
-testprime(GEN bnf, long minkowski)
-{
-  ulong av = avma;
-  long pp,i,nbideal,k,pmax;
-  GEN f,p1,vectpp,fb,dK, nf=checknf(bnf);
-  byteptr delta = diffptr;
-
-  if (DEBUGLEVEL>1)
-    fprintferr("PHASE 1: check primes to Zimmert bound = %ld\n\n",minkowski);
-  f=(GEN)nf[4]; dK=(GEN)nf[3];
-  if (!gcmp1(f))
-  {
-    GEN different = gmael(nf,5,5);
-    if (DEBUGLEVEL>1)
-      fprintferr("**** Testing Different = %Z\n",different);
-    p1 = isprincipalall(bnf,different,nf_FORCE);
-    if (DEBUGLEVEL>1) fprintferr("     is %Z\n",p1);
-  }
-  fb=(GEN)bnf[5];
-  p1 = gmael(fb, lg(fb)-1, 1); /* largest p in factorbase */
-  pp = 0; pmax = is_bigint(p1)? VERYBIGINT: itos(p1);
-  if ((ulong)minkowski > maxprime()) err(primer1);
-  while (pp < minkowski)
-  {
-    pp += *delta++;
-    if (DEBUGLEVEL>1) fprintferr("*** p = %ld\n",pp);
-    vectpp=primedec(bnf,stoi(pp)); nbideal=lg(vectpp)-1;
-    /* loop through all P | p if ramified, all but one otherwise */
-    if (!smodis(dK,pp)) nbideal++;
-    for (i=1; i<nbideal; i++)
-    {
-      GEN P = (GEN)vectpp[i];
-      if (DEBUGLEVEL>1)
-        fprintferr("  Testing P = %Z\n",P);
-      if (cmpis(idealnorm(bnf,P), minkowski) < 1)
-      {
-	if (pp <= pmax && (k = tablesearch(fb, P, cmp_prime_ideal)))
-	{
-	  if (DEBUGLEVEL>1) fprintferr("    #%ld in factor base\n",k);
-	}
-	else
-	{
-	  p1 = isprincipal(bnf,P);
-	  if (DEBUGLEVEL>1) fprintferr("    is %Z\n",p1);
-	}
-      }
-      else if (DEBUGLEVEL>1)
-        fprintferr("    Norm(P) > Zimmert bound\n");
-    }
-    avma = av;
-  }
-  if (DEBUGLEVEL>1) { fprintferr("End of PHASE 1.\n\n"); flusherr(); }
-}
-
 /* return \gamma_n^n if known, an upper bound otherwise */
 static GEN
 hermiteconstant(long n)
 {
   GEN h,h1;
-  long av;
+  gpmem_t av;
 
   switch(n)
   {
@@ -694,7 +770,7 @@ hermiteconstant(long n)
     case 8: return stoi(256);
   }
   av = avma;
-  h  = gpuigs(divsr(2,mppi(DEFAULTPREC)), n);
+  h  = gpowgs(divsr(2,mppi(DEFAULTPREC)), n);
   h1 = gsqr(ggamma(gdivgs(stoi(n+4),2),DEFAULTPREC));
   return gerepileupto(av, gmul(h,h1));
 }
@@ -731,33 +807,31 @@ isprimitive(GEN nf)
 }
 
 static GEN
+dft_bound()
+{
+  if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");
+  return dbltor(0.2);
+}
+
+static GEN
 regulatorbound(GEN bnf)
 {
-  long N,R1,R2,R;
-  GEN nf,dKa,bound,p1,c1;
+  long N, R1, R2, R;
+  GEN nf, dKa, p1, c1;
 
   nf = (GEN)bnf[7]; N = degpol(nf[1]);
-  bound = dbltor(0.2);
-  if (!isprimitive(nf))
-  {
-    if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");
-    return bound;
-  }
+  if (!isprimitive(nf)) return dft_bound();
+
   dKa = absi((GEN)nf[3]);
-  R1 = itos(gmael(nf,2,1));
-  R2 = itos(gmael(nf,2,2)); R = R1+R2-1;
-  if (!R2 && N<12) c1 = gpuigs(stoi(4),N>>1); else c1 = gpuigs(stoi(N),N);
-  if (cmpii(dKa,c1) <= 0)
-  {
-    if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");
-    return bound;
-  }
+  nf_get_sign(nf, &R1, &R2); R = R1+R2-1;
+  if (!R2 && N<12) c1 = gpowgs(stoi(4),N>>1); else c1 = gpowgs(stoi(N),N);
+  if (cmpii(dKa,c1) <= 0) return dft_bound();
+
   p1 = gsqr(glog(gdiv(dKa,c1),DEFAULTPREC));
-  p1 = gdivgs(gmul2n(gpuigs(gdivgs(gmulgs(p1,3),N*(N*N-1)-6*R2),R),R2),N);
-  p1 = gsqrt(gdiv(p1, hermiteconstant(R)), DEFAULTPREC);
-  if (gcmp(p1,bound) > 0) bound = p1;
+  p1 = divrs(gmul2n(gpowgs(divrs(mulrs(p1,3),N*(N*N-1)-6*R2),R),R2), N);
+  p1 = mpsqrt(gdiv(p1, hermiteconstant(R)));
   if (DEBUGLEVEL>1) fprintferr("Mahler bound for regulator: %Z\n",p1);
-  return bound;
+  return gmax(p1, dbltor(0.2));
 }
 
 /* x given by its embeddings */
@@ -780,23 +854,23 @@ norm_by_embed(long r1, GEN x)
 static int
 is_unit(GEN M, long r1, GEN x)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN Nx = ground( norm_by_embed(r1, gmul_mat_smallvec(M,x)) );
   int ok = is_pm1(Nx);
   avma = av; return ok;
 }
 
-#define NBMAX 5000
 /* FIXME: should use smallvectors */
 static GEN
-minimforunits(GEN nf, long BORNE, long stockmax)
+minimforunits(GEN nf, long BORNE, GEN w)
 {
   const long prec = MEDDEFAULTPREC;
-  long av = avma,n,i,j,k,s,norme,normax,*x,cmpt,r1;
-  GEN u,r,S,a,M,p1;
-  double p;
+  long n, i, j, k, s, *x, r1, cnt = 0;
+  gpmem_t av = avma;
+  GEN u,r,a,M;
+  double p, norme, normin, normax;
   double **q,*v,*y,*z;
-  double eps=0.000001, BOUND = BORNE + eps;
+  double eps=0.000001, BOUND = BORNE * 1.00001;
 
   if (DEBUGLEVEL>=2)
   {
@@ -804,23 +878,18 @@ minimforunits(GEN nf, long BORNE, long stockmax)
     if (DEBUGLEVEL>2) fprintferr("   BOUND = %ld\n",BOUND);
     flusherr();
   }
-  r1 = nf_get_r1(nf);
-  a = gmael(nf,5,3); n = lg(a);
-  minim_alloc(n, &q, &x, &y, &z, &v);
-  n--;
-  u = lllgram(a,prec);
-  M = gmul(gmael(nf,5,1), u); /* embeddings of T2-reduced basis */
-  M = gprec_w(M, prec);
-  a = gmul(qf_base_change(a,u,1), realun(prec));
-  r = sqred1(a);
+  r1 = nf_get_r1(nf); n = degpol(nf[1]);
+  minim_alloc(n+1, &q, &x, &y, &z, &v);
+  M = gprec_w(gmael(nf,5,1), prec);
+  a = gmul(gmael(nf,5,2), realun(prec));
+  r = sqred1_from_QR(a, prec);
   for (j=1; j<=n; j++)
   {
     v[j] = rtodbl(gcoeff(r,j,j));
     for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
   }
-  normax=0;
-  S = cgetg(stockmax+1,t_MAT);
-  s=0; k=n; cmpt=0; y[n]=z[n]=0;
+  normax = 0.; normin = (double)BOUND;
+  s=0; k=n; y[n]=z[n]=0;
   x[n]=(long)(sqrt(BOUND/v[n]));
 
   for(;;)
@@ -846,31 +915,27 @@ minimforunits(GEN nf, long BORNE, long stockmax)
     }
     while (k>1);
     if (!x[1] && y[1]<=eps) break;
-    if (++cmpt > NBMAX) { av=avma; return NULL; }
 
     if (DEBUGLEVEL>8){ fprintferr("."); flusherr(); }
-    p = x[1]+z[1]; norme = (long)(y[1] + p*p*v[1] + eps);
+    if (++cnt == 5000) return NULL; /* too expensive */
+
+    p = x[1]+z[1]; norme = y[1] + p*p*v[1] + eps;
     if (norme > normax) normax = norme;
-    if (is_unit(M,r1, x))
+    if (is_unit(M,r1, x)
+    && (norme > 2*n  /* exclude roots of unity */
+        || !isnfscalar(element_pow(nf, small_to_col(x), w))))
     {
+      if (norme < normin) normin = norme;
       if (DEBUGLEVEL>=2) { fprintferr("*"); flusherr(); }
-      cmpt = 0;
-      if (++s <= stockmax)
-      {
-	p1 = cgetg(n+1,t_COL);
-	for (i=1; i<=n; i++) p1[i]=lstoi(x[i]);
-	S[s] = lmul(u,p1);
-      }
     }
     x[k]--;
   }
   if (DEBUGLEVEL>=2){ fprintferr("\n"); flusherr(); }
-  k = (s<stockmax)? s:stockmax; setlg(S,k+1);
-  S = gerepilecopy(av, S);
-  u = cgetg(4,t_VEC);
+  avma = av; u = cgetg(4,t_VEC);
   u[1] = lstoi(s<<1);
-  u[2] = lstoi(normax);
-  u[3] = (long)S; return u;
+  u[2] = (long)dbltor(normax);
+  u[3] = (long)dbltor(normin);
+  return u;
 }
 
 #undef NBMAX
@@ -901,7 +966,7 @@ compute_M0(GEN M_star,long N)
     m2 = (N-n1)>>1;
     for (n2=n1; n2<=m2; n2++)
     {
-      long av = avma; n3=N-n1-n2; prec=PREC;
+      gpmem_t av = avma; n3=N-n1-n2; prec=PREC;
       if (n1==n2 && n1==n3) /* n1 = n2 = n3 = m1 = N/3 */
       {
 	p1=gdivgs(M_star,m1);
@@ -924,8 +989,8 @@ compute_M0(GEN M_star,long N)
       { /* n3 > N/3 >= n1 */
 	k = n2; kk = N-2*k;
 	p2=gsub(M_star,gmulgs(X,k));
-	p3=gmul(gpuigs(stoi(kk),kk),gpuigs(gsubgs(gmul(M_star,p2),kk*kk),k));
-	pol=gsub(p3,gmul(gmul(gpuigs(stoi(k),k),gpuigs(X,k)),gpuigs(p2,N-k)));
+	p3=gmul(gpowgs(stoi(kk),kk),gpowgs(gsubgs(gmul(M_star,p2),kk*kk),k));
+	pol=gsub(p3,gmul(gmul(gpowgs(stoi(k),k),gpowgs(X,k)),gpowgs(p2,N-k)));
 	prec=gprecision(pol); if (!prec) prec = MEDDEFAULTPREC;
 	r=roots(pol,prec); lr = lg(r);
 	for (i=1; i<lr; i++)
@@ -943,7 +1008,7 @@ compute_M0(GEN M_star,long N)
           if (gsigne(v) <= 0) continue;
 
           u=gmul2n(gadd(S,p6),-1);
-          w=gpui(P,gdivgs(stoi(-k),kk),prec);
+          w=gpow(P,gdivgs(stoi(-k),kk),prec);
           p6=gmulsg(k,gadd(gsqr(glog(u,prec)),gsqr(glog(v,prec))));
           M0_pro=gmul2n(gadd(p6,gmulsg(kk,gsqr(glog(w,prec)))),-2);
           if (DEBUGLEVEL>2)
@@ -961,7 +1026,7 @@ compute_M0(GEN M_star,long N)
 	f2 = gadd(f2,gmulsg(n2,gmul(X,Z)));
 	f2 = gadd(f2,gmulsg(n3,gmul(X,Y)));
 	f2 = gsub(f2,gmul(M,gmul(X,gmul(Y,Z))));
-	f3 = gsub(gmul(gpuigs(X,n1),gmul(gpuigs(Y,n2),gpuigs(Z,n3))), gun);
+	f3 = gsub(gmul(gpowgs(X,n1),gmul(gpowgs(Y,n2),gpowgs(Z,n3))), gun);
         /* f1 = n1 X + n2 Y + n3 Z - M */
         /* f2 = n1 YZ + n2 XZ + n3 XY */
         /* f3 = X^n1 Y^n2 Z^n3 - 1*/
@@ -1019,85 +1084,70 @@ compute_M0(GEN M_star,long N)
       if (!M0) avma = av; else M0 = gerepilecopy(av, M0);
     }
   }
-  for (i=1;i<=4;i++) delete_var();
+  for (i=1;i<=4;i++) (void)delete_var();
   return M0? M0: gzero;
 }
 
 static GEN
 lowerboundforregulator_i(GEN bnf)
 {
-  long N,R1,R2,RU,i,nbrootsofone,nbmin;
-  GEN rootsofone,nf,M0,M,m,col,T2,bound,minunit,newminunit;
-  GEN vecminim,colalg,p1,pol,y;
+  long N,R1,R2,RU,i;
+  GEN nf,M0,M,G,bound,minunit,newminunit;
+  GEN vecminim,p1,pol,y;
   GEN units = check_units(bnf,"bnfcertify");
 
-  rootsofone=gmael(bnf,8,4); nbrootsofone=itos((GEN)rootsofone[1]);
-  nf=(GEN)bnf[7]; T2=gmael(nf,5,3); N=degpol(nf[1]);
-  R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); RU=R1+R2-1;
-  if (RU==0) return gun;
+  nf = (GEN)bnf[7]; N = degpol(nf[1]);
+  nf_get_sign(nf, &R1, &R2); RU = R1+R2-1;
+  if (!RU) return gun;
 
-  units=algtobasis(bnf,units); minunit=qfeval(T2,(GEN)units[1]);
+  G = gmael(nf,5,2);
+  units = algtobasis(bnf,units);
+  minunit = gnorml2(gmul(G, (GEN)units[1])); /* T2(units[1]) */
   for (i=2; i<=RU; i++)
   {
-    newminunit=qfeval(T2,(GEN)units[i]);
-    if (gcmp(newminunit,minunit)<0) minunit=newminunit;
+    newminunit = gnorml2(gmul(G, (GEN)units[i]));
+    if (gcmp(newminunit,minunit) < 0) minunit = newminunit;
   }
-  if (gcmpgs(minunit,1000000000)>0) return NULL;
+  if (gexpo(minunit) > 30) return NULL;
 
-  vecminim = minimforunits(nf,itos(gceil(minunit)),10000);
+  vecminim = minimforunits(nf, itos(gceil(minunit)), gmael3(bnf,8,4,1));
   if (!vecminim) return NULL;
-  m=(GEN)vecminim[3]; nbmin=lg(m)-1;
-  if (nbmin==10000) return NULL;
-  bound=gaddgs(minunit,1);
-  for (i=1; i<=nbmin; i++)
-  {
-    col=(GEN)m[i]; colalg=basistoalg(nf,col);
-    if (!gcmp1(lift_intern(gpuigs(colalg,nbrootsofone))))
-    {
-      newminunit=qfeval(T2,col);
-      if (gcmp(newminunit,bound)<0) bound=newminunit;
-    }
-  }
-  if (gcmp(bound,minunit)>0) err(talker,"bug in lowerboundforregulator");
+  bound = (GEN)vecminim[3];
+  i = gexpo(gsub(bound,minunit));
+  if (i > -5) err(bugparier,"lowerboundforregulator");
+
   if (DEBUGLEVEL>1)
   {
-    fprintferr("M* = %Z\n",gprec_w(bound,3));
+    fprintferr("M* = %Z\n", bound);
     if (DEBUGLEVEL>2)
     {
-      p1=polx[0]; pol=gaddgs(gsub(gpuigs(p1,N),gmul(bound,p1)),N-1);
+      p1=polx[0]; pol=gaddgs(gsub(gpowgs(p1,N),gmul(bound,p1)),N-1);
       p1 = roots(pol,DEFAULTPREC);
       if (N&1) y=greal((GEN)p1[3]); else y=greal((GEN)p1[2]);
       M0 = gmul2n(gmulsg(N*(N-1),gsqr(glog(y,DEFAULTPREC))),-2);
       fprintferr("pol = %Z\n",pol);
       fprintferr("old method: y = %Z, M0 = %Z\n",y,gprec_w(M0,3));
     }
-    flusherr();
   }
   M0 = compute_M0(bound,N);
   if (DEBUGLEVEL>1) { fprintferr("M0 = %Z\n",gprec_w(M0,3)); flusherr(); }
-  M = gmul2n(gdivgs(gdiv(gpuigs(M0,RU),hermiteconstant(RU)),N),R2);
-  if (gcmp(M,dbltor(0.04))<0) return NULL;
+  M = gmul2n(gdivgs(gdiv(gpowgs(M0,RU),hermiteconstant(RU)),N),R2);
+  if (gcmp(M, dbltor(0.04)) < 0) return NULL;
   M = gsqrt(M,DEFAULTPREC);
   if (DEBUGLEVEL>1)
-  {
     fprintferr("(lower bound for regulator) M = %Z\n",gprec_w(M,3));
-    flusherr();
-  }
   return M;
 }
 
 static GEN
 lowerboundforregulator(GEN bnf)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN x = lowerboundforregulator_i(bnf);
   if (!x) { avma = av; x = regulatorbound(bnf); }
   return x;
 }
 
-extern GEN to_Fp_simple(GEN x, GEN prh);
-extern GEN Fp_PHlog(GEN a, GEN g, GEN p, GEN ord);
-
 /* Compute a square matrix of rank length(beta) associated to a family
  * (P_i), 1<=i<=length(beta), of primes s.t. N(P_i) = 1 mod pp, and
  * (P_i,beta[j]) = 1 for all i,j */
@@ -1105,7 +1155,7 @@ static void
 primecertify(GEN bnf,GEN beta,long pp,GEN big)
 {
   long i,j,qq,nbcol,lb,nbqq,ra,N;
-  GEN nf,mat,mat1,qgen,decqq,newcol,Qh,Q,g,ord;
+  GEN nf,mat,mat1,qgen,decqq,newcol,Q,g,ord,modpr;
 
   ord = NULL; /* gcc -Wall */
   nbcol = 0; nf = (GEN)bnf[7]; N = degpol(nf[1]);
@@ -1126,11 +1176,11 @@ primecertify(GEN bnf,GEN beta,long pp,GEN big)
         g = lift_intern(gener(qgen)); /* primitive root */
         ord = decomp(stoi(qq-1));
       }
-      Qh = prime_to_ideal(nf,Q);
+      modpr = zkmodprinit(nf, Q);
       newcol = cgetg(lb+1,t_COL);
       for (j=1; j<=lb; j++)
       {
-        GEN t = to_Fp_simple((GEN)beta[j], Qh);
+        GEN t = to_Fp_simple(nf, (GEN)beta[j], modpr);
         newcol[j] = (long)Fp_PHlog(t,g,qgen,ord);
       }
       if (DEBUGLEVEL>3)
@@ -1143,7 +1193,7 @@ primecertify(GEN bnf,GEN beta,long pp,GEN big)
       mat1 = concatsp(mat,newcol); ra = rank(mat1);
       if (ra==nbcol) continue;
 
-      if (DEBUGLEVEL>2) fprintferr("       new rank: %ld\n\n",ra);
+      if (DEBUGLEVEL>2) fprintferr("       new rank: %ld\n",ra);
       if (++nbcol == lb) return;
       mat = mat1;
     }
@@ -1153,7 +1203,7 @@ primecertify(GEN bnf,GEN beta,long pp,GEN big)
 static void
 check_prime(long p, GEN bnf, GEN cyc, GEN cycgen, GEN fu, GEN mu, GEN big)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,b, lc = lg(cyc), w = itos((GEN)mu[1]), lf = lg(fu);
   GEN beta = cgetg(lf+lc, t_VEC);
 
@@ -1178,16 +1228,16 @@ check_prime(long p, GEN bnf, GEN cyc, GEN cycgen, GEN 
 long
 certifybuchall(GEN bnf)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long nbgen,i,j,p,N,R1,R2,R,bound;
   GEN big,nf,reg,rootsofone,funits,gen,p1,gbound,cycgen,cyc;
   byteptr delta = diffptr;
 
   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
   N=degpol(nf[1]); if (N==1) return 1;
-  R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); R=R1+R2-1;
+  nf_get_sign(nf, &R1, &R2); R = R1+R2-1;
   funits = check_units(bnf,"bnfcertify");
-  testprime(bnf, zimmertbound(N,R2,absi((GEN)nf[3])));
+  testprimes(bnf, zimmertbound(N,R2,absi((GEN)nf[3])));
   reg = gmael(bnf,8,2);
   cyc = gmael3(bnf,8,1,2); nbgen = lg(cyc)-1;
   gen = gmael3(bnf,8,1,3); rootsofone = gmael(bnf,8,4);
@@ -1224,8 +1274,10 @@ certifybuchall(GEN bnf)
   rootsofone = dummycopy(rootsofone);
   rootsofone[2] = (long)algtobasis(nf, (GEN)rootsofone[2]);
 
-  for (p = *delta++; p <= bound; p += *delta++)
+  for (p = *delta++; p <= bound; ) {  
     check_prime(p,bnf,cyc,cycgen,funits,rootsofone,big);
+    NEXT_PRIME_VIADIFF(p, delta);
+  }
 
   if (nbgen)
   {
@@ -1266,7 +1318,7 @@ imageofgroup0(GEN gen,GEN bnr,GEN H)
 static GEN
 imageofgroup(GEN gen,GEN bnr,GEN H)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   return gerepileupto(av,imageofgroup0(gen,bnr,H));
 }
 
@@ -1327,18 +1379,22 @@ bnrisconductor(GEN arg0,GEN arg1,GEN arg2)
   return itos(conductor(bnr,sub,-1));
 }
 
+GEN
+isprincipalray_init(GEN bnf, GEN x)
+{
+  GEN z = cgetg(3,t_VEC);
+  z[2] = (long)quick_isprincipalgen(bnf, x);
+  z[1] = (long)x; return z;
+}
+
 /* special for isprincipalrayall */
 static GEN
 getgen(GEN bnf, GEN gen)
 {
   long i,l = lg(gen);
-  GEN p1, g = cgetg(l, t_VEC);
+  GEN g = cgetg(l, t_VEC);
   for (i=1; i<l; i++)
-  {
-    p1 = cgetg(3,t_VEC); g[i] = (long)p1;
-    p1[1] = (long)gen[i];
-    p1[2] = (long)quick_isprincipalgen(bnf, (GEN)gen[i]);
-  }
+    g[i] = (long)isprincipalray_init(bnf, (GEN)gen[i]);
   return g;
 }
 
@@ -1352,7 +1408,7 @@ getgen(GEN bnf, GEN gen)
 GEN
 conductor(GEN bnr, GEN H, long all)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long r1,j,k,ep;
   GEN bnf,nf,gen,bid,ideal,arch,p1,clhray,clhss,fa,arch2,bnr2,P,ex,mod;
 
@@ -1363,8 +1419,8 @@ conductor(GEN bnr, GEN H, long all)
   nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);
   ideal= gmael(bid,1,1);
   arch = gmael(bid,1,2);
-  if (gcmp0(H)) H = NULL;
-  else
+  if (H && gcmp0(H)) H = NULL;
+  if (H)
   {
     p1 = gauss(H, diagonal(gmael(bnr,5,2)));
     if (!gcmp1(denom(p1))) err(talker,"incorrect subgroup in conductor");
@@ -1384,7 +1440,7 @@ conductor(GEN bnr, GEN H, long all)
     ep = (all>=0)? itos((GEN)ex[k]): 1;
     for (j=1; j<=ep; j++)
     {
-      mod[1] = (long)idealdivexact(nf,ideal,pr);
+      mod[1] = (long)idealdivpowprime(nf,ideal,pr,gun);
       clhss = orderofquotient(bnf,mod,H,gen);
       if (!egalii(clhss,clhray)) break;
       if (all < 0) { avma = av; return gzero; }
@@ -1414,22 +1470,21 @@ conductor(GEN bnr, GEN H, long all)
 
 /* etant donne un bnr et un polynome relatif, trouve le groupe des normes
    correspondant a l'extension relative en supposant qu'elle est abelienne
-   et que le module donnant bnr est multiple du conducteur. Verifie que
-   l'extension est bien abelienne (sous GRH) si rnf != NULL, dans ce cas
-   rnf est le rnf de l'extension relative. */
-static GEN
-rnfnormgroup0(GEN bnr, GEN polrel, GEN rnf)
+   et que le module donnant bnr est multiple du conducteur. */
+GEN
+rnfnormgroup(GEN bnr, GEN polrel)
 {
-  long av=avma,i,j,reldeg,sizemat,p,pmax,nfac,k;
-  GEN bnf,polreldisc,discnf,nf,raycl,group,detgroup,fa,greldeg;
-  GEN contreld,primreld,reldisc,famo,ep,fac,col,p1,bd,upnf;
-  byteptr d = diffptr + 1; /* start at p = 2 */
+  long i, j, reldeg, sizemat, p, nfac, k;
+  gpmem_t av = avma;
+  GEN bnf,index,discnf,nf,raycl,group,detgroup,fa,greldeg;
+  GEN famo,ep,fac,col;
+  byteptr d = diffptr;
 
   checkbnr(bnr); bnf=(GEN)bnr[1]; raycl=(GEN)bnr[5];
   nf=(GEN)bnf[7];
   polrel = fix_relative_pol(nf,polrel,1);
   if (typ(polrel)!=t_POL) err(typeer,"rnfnormgroup");
-  reldeg=degpol(polrel);
+  reldeg = degpol(polrel);
   /* reldeg-th powers are in norm group */
   greldeg = stoi(reldeg);
   group = diagonal(gmod((GEN)raycl[2], greldeg));
@@ -1437,124 +1492,95 @@ rnfnormgroup0(GEN bnr, GEN polrel, GEN rnf)
     if (!signe(gcoeff(group,i,i))) coeff(group,i,i) = (long)greldeg;
   detgroup = dethnf_i(group);
   k = cmpis(detgroup,reldeg);
-  if (k<0)
-  {
-    if (rnf) return NULL;
+  if (k < 0)
     err(talker,"not an Abelian extension in rnfnormgroup?");
-  }
-  if (!rnf && !k) return gerepilecopy(av, group);
+  if (!k) return gerepilecopy(av, group);
 
-  polreldisc=discsr(polrel);
-
-  if (rnf)
-  {
-    reldisc=gmael(rnf,3,1);
-    upnf=nfinit0(gmael(rnf,11,1),1,DEFAULTPREC);
-  }
-  else
-  {
-    reldisc = idealhermite(nf,polreldisc);
-    upnf = NULL;
-  }
-
-  reldisc = idealmul(nf, reldisc, gmael3(bnr,2,1,1));
-  contreld= content(reldisc);
-  primreld= gcmp1(contreld)? reldisc: gdiv(reldisc, contreld);
-
   discnf = (GEN)nf[3];
-  k = degpol(nf[1]);
-  bd = gmulsg(k, glog(absi(discnf), DEFAULTPREC));
-  bd = gadd(bd,glog(mpabs(det(reldisc)),DEFAULTPREC));
-  p1 = dbltor(reldeg * k * 2.5 + 5);
-  bd = gfloor(gsqr(gadd(gmulsg(4,bd),p1)));
-
-  pmax = is_bigint(bd)? 0: itos(bd);
-  if (rnf)
-  {
-    if (DEBUGLEVEL) fprintferr("rnfnormgroup: bound for primes = %Z\n", bd);
-    if (!pmax) err(warner,"rnfnormgroup: prime bound too large, can't certify");
-  }
+  index  = (GEN)nf[4];
   sizemat=lg(group)-1;
-  for (p=2; !pmax || p < pmax; p += *d++)
+  for (p=0 ;;)
   {
     long oldf = -1, lfa;
     /* If all pr are unramified and have the same residue degree, p =prod pr
      * and including last pr^f or p^f is the same, but the last isprincipal
      * is much easier! oldf is used to track this */
 
-    if (!*d) err(primer1);
-    if (!smodis(contreld,p)) continue; /* all pr|p ramified */
+    NEXT_PRIME_VIADIFF_CHECK(p,d);
+    if (!smodis(index, p)) continue; /* can't be treated efficiently */
 
-    fa = primedec(nf,stoi(p)); lfa = lg(fa)-1;
-
+    fa = primedec(nf, stoi(p)); lfa = lg(fa)-1;
     for (i=1; i<=lfa; i++)
     {
-      GEN pr = (GEN)fa[i];
+      GEN pr = (GEN)fa[i], pp, T, polr;
+      GEN modpr = nf_to_ff_init(nf, &pr, &T, &pp);
       long f;
+
+      polr = modprX(polrel, nf, modpr);
+      /* if pr (probably) ramified, we have to use all (non-ram) P | pr */
+      if (!FqX_is_squarefree(polr, T,pp)) { oldf = 0; continue; }
+
+      famo = FqX_factor(polr, T, pp);
+      fac = (GEN)famo[1]; f = degpol((GEN)fac[1]);
+      ep  = (GEN)famo[2]; nfac = lg(ep)-1;
       /* check decomposition of pr has Galois type */
-      if (element_val(nf,polreldisc,pr) != 0)
+      for (j=1; j<=nfac; j++)
       {
-        /* if pr ramified, we will have to use all (non-ram) P | pr */
-        if (idealval(nf,primreld,pr)!=0) { oldf = 0; continue; }
-
-        famo=idealfactor(upnf,rnfidealup(rnf,pr));	
-        ep=(GEN)famo[2];
-        fac=(GEN)famo[1];
-        nfac=lg(ep)-1;
-        f = itos(gmael(fac,1,4));
-        for (j=1; j<=nfac; j++)
-        {
-          if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");
-          if (itos(gmael(fac,j,4)) != f)
-          {
-            if (rnf) return NULL;
-            err(talker,"non Galois extension in rnfnormgroup");
-          }
-        }
+        if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");
+        if (degpol(fac[j]) != f)
+          err(talker,"non Galois extension in rnfnormgroup");
       }
-      else
-      {
-        famo=nffactormod(nf,polrel,pr);
-        ep=(GEN)famo[2];
-        fac=(GEN)famo[1];
-        nfac=lg(ep)-1;
-        f = degpol((GEN)fac[1]);
-        for (j=1; j<=nfac; j++)
-        {
-          if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");
-          if (degpol(fac[j]) != f)
-          {
-            if (rnf) return NULL;
-            err(talker,"non Galois extension in rnfnormgroup");
-          }
-        }
-      }
       if (oldf < 0) oldf = f; else if (oldf != f) oldf = 0;
       if (f == reldeg) continue; /* reldeg-th powers already included */
 
       if (oldf && i == lfa && !smodis(discnf, p)) pr = stoi(p);
 
       /* pr^f = N P, P | pr, hence is in norm group */
-      col = gmulsg(f, isprincipalrayall(bnr,pr,nf_REGULAR));
+      col = gmulsg(f, isprincipalrayall(bnr,pr,0));
       group = hnf(concatsp(group, col));
       detgroup = dethnf_i(group);
       k = cmpis(detgroup,reldeg);
       if (k < 0)
-      {
-        if (rnf) return NULL;
         err(talker,"not an Abelian extension in rnfnormgroup");
-      }
-      if (!rnf && !k) { cgiv(detgroup); return gerepileupto(av,group); }
+      if (!k) { cgiv(detgroup); return gerepileupto(av,group); }
     }
   }
   if (k>0) err(bugparier,"rnfnormgroup");
   cgiv(detgroup); return gerepileupto(av,group);
 }
 
-GEN
-rnfnormgroup(GEN bnr, GEN polrel)
+static int
+rnf_is_abelian(GEN nf, GEN pol)
 {
-  return rnfnormgroup0(bnr,polrel,NULL);
+  GEN ro, rores, nfL, L, mod, d;
+  long i,j, l;
+
+  nf = checknf(nf);
+  L = rnfequation(nf,pol);
+  mod = dummycopy(L); setvarn(mod, varn(nf[1]));
+  nfL = _initalg(mod, nf_PARTIALFACT, DEFAULTPREC);
+  ro = nfroots(nfL, L);
+  l = lg(ro)-1;
+  if (l != degpol(pol)) return 0;
+  if (isprime(stoi(l))) return 1;
+  ro = Q_remove_denom(ro, &d);
+  if (!d) rores = ro;
+  else
+  {
+    rores = cgetg(l+1, t_VEC);
+    for (i=1; i<=l; i++)
+      rores[i] = (long)rescale_pol((GEN)ro[i], d);
+  }
+  /* assume roots are sorted by increasing degree */
+  for (i=1; i<l; i++)
+    for (j=1; j<i; j++)
+    {
+      GEN a = RX_RXQ_compo((GEN)rores[j], (GEN)ro[i], mod);
+      GEN b = RX_RXQ_compo((GEN)rores[i], (GEN)ro[j], mod);
+      if (d) a = gmul(a, gpowgs(d, degpol(ro[i]) - degpol(ro[j])));
+      if (!gegal(a, b)) return 0;
+    }
+  return 1;
 }
 
 /* Etant donne un bnf et un polynome relatif polrel definissant une extension
@@ -1562,38 +1588,29 @@ rnfnormgroup(GEN bnr, GEN polrel)
    a l'extension definie par polrel sous la forme
    [[ideal,arch],[hm,cyc,gen],group] ou [ideal,arch] est le conducteur, et
    [hm,cyc,gen] est le groupe de classes de rayon correspondant. Verifie
-   (sous GRH) que polrel definit bien une extension abelienne si flag != 0 */
+   que polrel definit bien une extension abelienne si flag != 0 */
 GEN
 rnfconductor(GEN bnf, GEN polrel, long flag)
 {
-  long av=avma,tetpil,R1,i,v;
-  GEN nf,module,rnf,arch,bnr,group,p1,pol2;
+  long R1, i;
+  gpmem_t av = avma;
+  GEN nf,module,arch,bnr,group,p1,pol2;
 
-  bnf = checkbnf(bnf); nf=(GEN)bnf[7];
+  bnf = checkbnf(bnf); nf = (GEN)bnf[7];
   if (typ(polrel)!=t_POL) err(typeer,"rnfconductor");
-  module=cgetg(3,t_VEC); R1=itos(gmael(nf,2,1));
-  v=varn(polrel);
-  p1=unifpol((GEN)bnf[7],polrel,0);
+  p1=unifpol(nf, polrel, 0);
   p1=denom(gtovec(p1));
-  pol2=gsubst(polrel,v,gdiv(polx[v],p1));
-  pol2=gmul(pol2,gpuigs(p1,degpol(pol2)));
-  if (flag)
-  {
-    rnf=rnfinitalg(bnf,pol2,DEFAULTPREC);
-    module[1]=mael(rnf,3,1);
-  }
-  else
-  {
-    rnf=NULL;
-    module[1]=rnfdiscf(nf,pol2)[1];
-  }
-  arch=cgetg(R1+1,t_VEC);
-  module[2]=(long)arch; for (i=1; i<=R1; i++) arch[i]=un;
-  bnr=buchrayall(bnf,module,nf_INIT | nf_GEN);
-  group=rnfnormgroup0(bnr,pol2,rnf);
+  pol2=rescale_pol(polrel, p1);
+  if (flag && !rnf_is_abelian(bnf, pol2)) { avma = av; return gzero; }
+
+  module = cgetg(3,t_VEC);
+  module[1] = rnfdiscf(nf,pol2)[1];
+  R1 = nf_get_r1(nf); arch = cgetg(R1+1,t_VEC);
+  module[2] = (long)arch; for (i=1; i<=R1; i++) arch[i]=un;
+  bnr   = buchrayall(bnf,module,nf_INIT | nf_GEN);
+  group = rnfnormgroup(bnr,pol2);
   if (!group) { avma = av; return gzero; }
-  tetpil=avma;
-  return gerepile(av,tetpil,conductor(bnr,group,1));
+  return gerepileupto(av, conductor(bnr,group,1));
 }
 
 /* Given a number field bnf=bnr[1], a ray class group structure
@@ -1604,7 +1621,7 @@ rnfconductor(GEN bnf, GEN polrel, long flag)
 static GEN
 discrayrelall(GEN bnr, GEN H, long flag)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long r1,j,k,ep, nz, flrel = flag & nf_REL, flcond = flag & nf_COND;
   GEN bnf,nf,gen,bid,ideal,arch,p1,clhray,clhss,fa,arch2,idealrel,P,ex,mod,dlk;
 
@@ -1639,13 +1656,13 @@ discrayrelall(GEN bnr, GEN H, long flag)
     mod[1] = (long)ideal;
     for (j=1; j<=ep; j++)
     {
-      mod[1] = (long)idealdivexact(nf,(GEN)mod[1],pr);
+      mod[1] = (long)idealdivpowprime(nf,(GEN)mod[1],pr,gun);
       clhss = orderofquotient(bnf,mod,H,gen);
       if (flcond && j==1 && egalii(clhss,clhray)) { avma = av; return gzero; }
       if (is_pm1(clhss)) { S = addis(S, ep-j+1); break; }
       S = addii(S, clhss);
     }
-    idealrel = flrel? idealmul(nf,idealrel, idealpow(nf,pr, S))
+    idealrel = flrel? idealmulpowprime(nf,idealrel, pr, S)
                     : mulii(idealrel, powgi(idealnorm(nf,pr),S));
   }
   dlk = flrel? idealdivexact(nf,idealpow(nf,ideal,clhray), idealrel)
@@ -1673,7 +1690,7 @@ discrayrelall(GEN bnr, GEN H, long flag)
 static GEN
 discrayabsall(GEN bnr, GEN subgroup,long flag)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long degk,clhray,n,R1;
   GEN z,p1,D,dk,nf,dkabs,bnf;
 
@@ -1717,7 +1734,7 @@ discrayrelcond(GEN bnr, GEN H)
 GEN
 discrayabs(GEN bnr, GEN H)
 {
-  return discrayabsall(bnr,H,nf_REGULAR);
+  return discrayabsall(bnr,H,0);
 }
 
 GEN
@@ -1732,9 +1749,9 @@ discrayabscond(GEN bnr, GEN H)
 GEN
 bnrconductorofchar(GEN bnr, GEN chi)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long nbgen,i;
-  GEN p1,m,U,d1,cyc;
+  GEN m,U,d1,cyc;
 
   checkbnrgen(bnr);
   cyc = gmael(bnr,5,2); nbgen = lg(cyc)-1;
@@ -1746,11 +1763,10 @@ bnrconductorofchar(GEN bnr, GEN chi)
   for (i=1; i<=nbgen; i++)
   {
     if (typ(chi[i]) != t_INT) err(typeer,"conductorofchar");
-    p1 = cgetg(2,t_COL); m[i] = (long)p1;
-    p1[1] = lmulii((GEN)chi[i], divii(d1, (GEN)cyc[i]));
+    m[i] = (long)_col(mulii((GEN)chi[i], divii(d1, (GEN)cyc[i])));
   }
-  p1 = cgetg(2,t_COL); m[i] = (long)p1;
-  p1[1] = (long)d1; U = (GEN)hnfall(m)[2];
+  m[i] = (long)_col(d1);
+  U = (GEN)hnfall(m)[2];
   setlg(U,nbgen+1);
   for (i=1; i<=nbgen; i++) setlg(U[i],nbgen+1); /* U = Ker chi */
   return gerepileupto(av, conductor(bnr,U,0));
@@ -1761,7 +1777,7 @@ bnrconductorofchar(GEN bnr, GEN chi)
 GEN
 rayclassnolist(GEN bnf,GEN lists)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,j,lx,ly;
   GEN blist,ulist,Llist,h,b,u,L,m;
 
@@ -1875,7 +1891,7 @@ factorpow(GEN fa,long n)
 GEN
 discrayabslist(GEN bnf,GEN lists)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long ii,jj,i,j,k,clhss,ep,clhray,lx,ly,r1,degk,nz;
   long n,R1,lP;
   GEN hlist,blist,dlist,nf,dkabs,b,h,d;
@@ -2016,7 +2032,8 @@ zsimp(GEN bid, GEN matunit)
 static GEN
 zsimpjoin(GEN bidsimp, GEN bidp, GEN dummyfa, GEN matunit)
 {
-  long i,l1,l2,nbgen,c, av = avma;
+  long i, l1, l2, nbgen, c;
+  gpmem_t av = avma;
   GEN U,U1,U2,cyc1,cyc2,cyc,u1u2,met, y = cgetg(5,t_VEC);
 
   y[1] = (long)vconcat((GEN)bidsimp[1],dummyfa);
@@ -2068,8 +2085,6 @@ rayclassnointern(GEN blist, GEN h)
   return L;
 }
 
-void rowselect_p(GEN A, GEN B, GEN p, long init);
-
 static GEN
 rayclassnointernarch(GEN blist, GEN h, GEN matU)
 {
@@ -2113,7 +2128,8 @@ rayclassnointernarch(GEN blist, GEN h, GEN matU)
 GEN
 decodemodule(GEN nf, GEN fa)
 {
-  long n,k,j,fauxpr,av=avma;
+  long n, k, j, fauxpr;
+  gpmem_t av=avma;
   GEN g,e,id,pr;
 
   nf = checknf(nf);
@@ -2127,7 +2143,7 @@ decodemodule(GEN nf, GEN fa)
     fauxpr = itos((GEN)g[k]);
     j = (fauxpr%n)+1; fauxpr /= n*n;
     pr = (GEN)primedec(nf,stoi(fauxpr))[j];
-    id = idealmul(nf,id, idealpow(nf,pr,(GEN)e[k]));
+    id = idealmulpowprime(nf,id, pr,(GEN)e[k]);
   }
   return gerepileupto(av,id);
 }
@@ -2157,7 +2173,8 @@ discrayabslistarchintern(GEN bnf, GEN arch, long bound
   long degk,i,j,k,p2s,lfa,lp1,sqbou,cex, allarch;
   long ffs,fs,resp,flbou,nba, k2,karch,kka,nbarch,jjj,jj,square;
   long ii2,ii,ly,clhray,lP,ep,S,clhss,normps,normi,nz,r1,R1,n,c;
-  ulong q, av0 = avma, av,av1,lim;
+  ulong q;
+  gpmem_t av0 = avma, av, av1, lim;
   GEN nf,p,z,p1,p2,p3,fa,pr,normp,ideal,bidp,z2,matarchunit;
   GEN funits,racunit,embunit,sous,clh,sousray,raylist;
   GEN clhrayall,discall,faall,Id,idealrel,idealrelinit;
@@ -2170,7 +2187,7 @@ discrayabslistarchintern(GEN bnf, GEN arch, long bound
   mod = Id = dlk = ideal = clhrayall = discall = faall = NULL;
 
   /* ce qui suit recopie d'assez pres ideallistzstarall */
-  if (DEBUGLEVEL>2) timer2();
+  if (DEBUGLEVEL>2) (void)timer2();
   bnf = checkbnf(bnf); flbou=0;
   nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);
   degk = degpol(nf[1]);
@@ -2219,7 +2236,7 @@ discrayabslistarchintern(GEN bnf, GEN arch, long bound
   if (bound > (long)maxprime()) err(primer1);
   for (p[2]=0; p[2]<=bound; )
   {
-    p[2] += *ptdif++;
+    NEXT_PRIME_VIADIFF(p[2], ptdif);
     if (!flbou && p[2]>sqbou)
     {
       if (DEBUGLEVEL>1) fprintferr("\nStarting rayclassno computations\n");
@@ -2387,10 +2404,10 @@ discrayabslistarchintern(GEN bnf, GEN arch, long bound
           if (!allarch && nba)
           {
             p1 = (GEN)primedec(nf,gprime)[ffs%degk+1];
-            ideal = idealmul(nf,ideal,idealpow(nf,p1,(GEN)ex[k]));
+            ideal = idealmulpowprime(nf,ideal,p1,(GEN)ex[k]);
           }
           S=0; clhss=0;
-          normi = ii; normps= itos(gpuigs(gprime,resp));
+          normi = ii; normps= itos(gpowgs(gprime,resp));
           for (j=1; j<=ep; j++)
           {
             GEN fad, fad1, fad2;
@@ -2485,78 +2502,82 @@ GEN
 discrayabslistarchsquare(GEN bnf, GEN arch, long bound)
 { return discrayabslistarchintern(bnf,arch,bound, -1); }
 
-static GEN
-computehuv(GEN bnr,GEN id, GEN arch)
+/* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the
+   matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */
+GEN
+bnrGetSurj(GEN bnr1, GEN bnr2)
 {
-  long i,nbgen, av = avma;
-  GEN bnf,bnrnew,listgen,P,U,DC;
-  GEN newmod=cgetg(3,t_VEC);
-  newmod[1]=(long)id;
-  newmod[2]=(long)arch;
+  long nbg, i;
+  GEN gen, M;
 
-  bnf=(GEN)bnr[1];
-  bnrnew=buchrayall(bnf,newmod,nf_INIT);
-  listgen=gmael(bnr,5,3); nbgen=lg(listgen);
-  P=cgetg(nbgen,t_MAT);
-  for (i=1; i<nbgen; i++)
-    P[i] = (long)isprincipalray(bnrnew,(GEN)listgen[i]);
-  DC=diagonal(gmael(bnrnew,5,2));
-  U=(GEN)hnfall(concatsp(P,DC))[2];
-  setlg(U,nbgen); for (i=1; i<nbgen; i++) setlg(U[i], nbgen);
-  return gerepileupto(av, hnf(concatsp(U,diagonal(gmael(bnr,5,2)))));
+  gen = gmael(bnr1, 5, 3);
+  nbg = lg(gen) - 1;
+
+  M = cgetg(nbg + 1, t_MAT);
+  for (i = 1; i <= nbg; i++)
+    M[i] = (long)isprincipalray(bnr2, (GEN)gen[i]);
+
+  return M;
 }
 
-/* 0 if hinv*list[i] has a denominator for all i, 1 otherwise */
-static int
-hnflistdivise(GEN list,GEN h)
+/* Kernel of the above map */
+static GEN
+bnrGetKer(GEN bnr, GEN mod2)
 {
-  long av = avma, i, I = lg(list);
-  GEN hinv = ginv(h);
+  long i, n;
+  gpmem_t av = avma;
+  GEN P,U, bnr2 = buchrayall(bnr,mod2,nf_INIT);
 
-  for (i=1; i<I; i++)
-    if (gcmp1(denom(gmul(hinv,(GEN)list[i])))) break;
-  avma = av; return i < I;
+  P = bnrGetSurj(bnr, bnr2); n = lg(P);
+  P = concatsp(P, diagonal(gmael(bnr2,5,2)));
+  U = (GEN)hnfall(P)[2]; setlg(U,n);
+  for (i=1; i<n; i++) setlg(U[i], n);
+  U = concatsp(U, diagonal(gmael(bnr, 5,2)));
+  return gerepileupto(av, hnf(U));
 }
 
 static GEN
-subgroupcond(GEN bnr, long indexbound)
+subgroupcond(GEN bnr, GEN indexbound)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,j,lgi,lp;
-  GEN li,p1,lidet,perm,nf,bid,ideal,arch,primelist,listkernels;
+  GEN li,p1,lidet,perm,nf,bid,ideal,arch,arch2,primelist,listKer;
+  GEN mod = cgetg(3, t_VEC);
 
   checkbnrgen(bnr); bid=(GEN)bnr[2];
   ideal=gmael(bid,1,1);
   arch =gmael(bid,1,2); nf=gmael(bnr,1,7);
   primelist=gmael(bid,3,1); lp=lg(primelist)-1;
-  listkernels=cgetg(lp+lg(arch),t_VEC);
+  mod[2] = (long)arch;
+  listKer=cgetg(lp+lg(arch),t_VEC);
   for (i=1; i<=lp; )
   {
-    p1=idealdiv(nf,ideal,(GEN)primelist[i]);
-    listkernels[i++]=(long)computehuv(bnr,p1,arch);
+    mod[1] = (long)idealdivpowprime(nf,ideal,(GEN)primelist[i],gun);
+    listKer[i++] = (long)bnrGetKer(bnr,mod);
   }
+  mod[1] = (long)ideal; arch2 = dummycopy(arch);
+  mod[2] = (long)arch2;
   for (j=1; j<lg(arch); j++)
     if (signe((GEN)arch[j]))
     {
-      p1=dummycopy(arch); p1[j]=zero;
-      listkernels[i++]=(long)computehuv(bnr,ideal,p1);
+      arch2[j] = zero; 
+      listKer[i++] = (long)bnrGetKer(bnr,mod);
+      arch2[j] = un;
     }
-  setlg(listkernels,i);
+  setlg(listKer,i);
 
-  li=subgrouplist(gmael(bnr,5,2),indexbound);
-  lgi=lg(li);
-  for (i=1,j=1; j<lgi; j++)
-    if (!hnflistdivise(listkernels, (GEN)li[j])) li[i++] = li[j];
+  li = subgroupcondlist(gmael(bnr,5,2),indexbound,listKer);
+  lgi = lg(li);
   /* sort by increasing index */
-  lgi = i; setlg(li,i); lidet=cgetg(lgi,t_VEC);
-  for (i=1; i<lgi; i++) lidet[i]=(long)dethnf_i((GEN)li[i]);
-  perm = sindexsort(lidet); p1=li; li=cgetg(lgi,t_VEC);
+  lidet = cgetg(lgi,t_VEC);
+  for (i=1; i<lgi; i++) lidet[i] = (long)dethnf_i((GEN)li[i]);
+  perm = sindexsort(lidet); p1 = li; li = cgetg(lgi,t_VEC);
   for (i=1; i<lgi; i++) li[i] = p1[perm[lgi-i]];
   return gerepilecopy(av,li);
 }
 
 GEN
-subgrouplist0(GEN bnr, long indexbound, long all)
+subgrouplist0(GEN bnr, GEN indexbound, long all)
 {
   if (typ(bnr)!=t_VEC) err(typeer,"subgrouplist");
   if (lg(bnr)!=1 && typ(bnr[1])!=t_INT)