===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/base3.c,v
retrieving revision 1.1.1.1
retrieving revision 1.2
diff -u -p -r1.1.1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/base3.c	2001/10/02 11:17:02	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/base3.c	2002/09/11 07:26:49	1.2
@@ -1,4 +1,4 @@
-/* $Id: base3.c,v 1.1.1.1 2001/10/02 11:17:02 noro Exp $
+/* $Id: base3.c,v 1.2 2002/09/11 07:26:49 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -19,11 +19,10 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /*                                                                 */
 /*******************************************************************/
 #include "pari.h"
-extern int ker_trivial_mod_p(GEN x, GEN p);
-extern long ideal_is_zk(GEN ideal,long N);
+extern GEN u_FpM_deplin(GEN x, ulong p);
+extern GEN u_FpM_inv(GEN a, ulong p);
 extern GEN famat_ideallog(GEN nf, GEN g, GEN e, GEN bid);
 extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);
-extern GEN hnfall_i(GEN A, GEN *ptB, long remove);
 extern GEN vconcat(GEN A, GEN B);
 
 /*******************************************************************/
@@ -39,6 +38,7 @@ isnfscalar(GEN x)
 {
   long lx=lg(x),i;
 
+  if (typ(x) != t_COL) return 0;
   for (i=2; i<lx; i++)
     if (!gcmp0((GEN)x[i])) return 0;
   return 1;
@@ -47,7 +47,7 @@ isnfscalar(GEN x)
 static GEN
 scal_mul(GEN nf, GEN x, GEN y, long ty)
 {
-  long av=avma, tetpil;
+  gpmem_t av=avma, tetpil;
   GEN p1;
 
   if (!is_extscalar_t(ty))
@@ -66,24 +66,14 @@ get_tab(GEN nf, long *N)
   *N = lg(tab[1])-1; return tab;
 }
 
-/* product of x and y in nf */
 GEN
-element_mul(GEN nf, GEN x, GEN y)
+mul_by_tab(GEN tab, GEN x, GEN y)
 {
-  long av,i,j,k,N,tx,ty;
-  GEN s,v,c,p1,tab;
+  gpmem_t av;
+  long i,j,k,N;
+  GEN s,v,c,p1;
 
-  if (x == y) return element_sqr(nf,x);
-
-  tx=typ(x); ty=typ(y);
-  nf=checknf(nf); tab = get_tab(nf, &N);
-  if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_mul");
-  if (ty==t_POLMOD) y=checknfelt_mod(nf,y,"element_mul");
-  if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
-  if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
-  if (isnfscalar(x)) return gmul((GEN)x[1],y);
-  if (isnfscalar(y)) return gmul((GEN)y[1],x);
-
+  N = lg(x)-1;
   v=cgetg(N+1,t_COL); av=avma;
   for (k=1; k<=N; k++)
   {
@@ -118,17 +108,41 @@ element_mul(GEN nf, GEN x, GEN y)
   return v;
 }
 
+/* product of x and y in nf */
+GEN
+element_mul(GEN nf, GEN x, GEN y)
+{
+  long N,tx,ty;
+  GEN tab;
+
+  if (x == y) return element_sqr(nf,x);
+
+  tx=typ(x); ty=typ(y);
+  nf=checknf(nf);
+  if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_mul");
+  if (ty==t_POLMOD) y=checknfelt_mod(nf,y,"element_mul");
+  if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
+  if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
+  if (tx != t_COL || ty != t_COL) err(typeer,"element_mul");
+  if (isnfscalar(x)) return gmul((GEN)x[1],y);
+  if (isnfscalar(y)) return gmul((GEN)y[1],x);
+
+  tab = get_tab(nf, &N);
+  return mul_by_tab(tab,x,y);
+}
+
 /* inverse of x in nf */
 GEN
 element_inv(GEN nf, GEN x)
 {
-  long av=avma,i,N,tx=typ(x);
+  gpmem_t av=avma;
+  long i,N,tx=typ(x);
   GEN p1,p;
 
   nf=checknf(nf); N=degpol(nf[1]);
   if (is_extscalar_t(tx))
   {
-    if (tx==t_POLMOD) checknfelt_mod(nf,x,"element_inv");
+    if (tx==t_POLMOD) (void)checknfelt_mod(nf,x,"element_inv");
     else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
     return gerepileupto(av, algtobasis(nf, ginv(x)));
   }
@@ -138,6 +152,7 @@ element_inv(GEN nf, GEN x)
     for (i=2; i<=N; i++) p1[i]=lcopy((GEN)x[i]);
     return p1;
   }
+  if (tx != t_COL) err(typeer,"element_inv");
   p = NULL;
   for (i=1; i<=N; i++)
     if (typ(x[i])==t_INTMOD)
@@ -145,8 +160,8 @@ element_inv(GEN nf, GEN x)
       p = gmael(x,i,1);
       x = lift(x); break;
     }
-  p1 = ginvmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
-  p1 = algtobasis_intern(nf,p1);
+  p1 = QX_invmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
+  p1 = algtobasis_i(nf,p1);
 
   if (p) p1 = FpV(p1, p);
   return gerepileupto(av,p1);
@@ -156,14 +171,15 @@ element_inv(GEN nf, GEN x)
 GEN
 element_div(GEN nf, GEN x, GEN y)
 {
-  long av=avma,i,N,tx=typ(x),ty=typ(y);
+  gpmem_t av=avma;
+  long i,N,tx=typ(x),ty=typ(y);
   GEN p1,p;
 
   nf=checknf(nf); N=degpol(nf[1]);
-  if (tx==t_POLMOD) checknfelt_mod(nf,x,"element_div");
+  if (tx==t_POLMOD) (void)checknfelt_mod(nf,x,"element_div");
   else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
 
-  if (ty==t_POLMOD) checknfelt_mod(nf,y,"element_div");
+  if (ty==t_POLMOD) (void)checknfelt_mod(nf,y,"element_div");
   else if (ty==t_POL) y=gmodulcp(y,(GEN)nf[1]);
 
   if (is_extscalar_t(tx))
@@ -182,6 +198,7 @@ element_div(GEN nf, GEN x, GEN y)
     p1=gdiv(gmodulcp(gmul((GEN)nf[7],x),(GEN)nf[1]),y);
     return gerepileupto(av, algtobasis(nf,p1));
   }
+  if (tx != t_COL || ty != t_COL) err(typeer,"element_div");
 
   if (isnfscalar(y)) return gdiv(x,(GEN)y[1]);
   if (isnfscalar(x))
@@ -206,8 +223,8 @@ element_div(GEN nf, GEN x, GEN y)
       y = lift(y); break;
     }
 
-  p1 = gmul(gmul((GEN)nf[7],x), ginvmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
-  p1 = algtobasis_intern(nf, gres(p1, (GEN)nf[1]));
+  p1 = gmul(gmul((GEN)nf[7],x), QX_invmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
+  p1 = algtobasis_i(nf, gres(p1, (GEN)nf[1]));
   if (p) p1 = FpV(p1,p);
   return gerepileupto(av,p1);
 }
@@ -216,10 +233,13 @@ element_div(GEN nf, GEN x, GEN y)
 GEN
 element_muli(GEN nf, GEN x, GEN y)
 {
-  long av,i,j,k,N;
+  gpmem_t av;
+  long i,j,k,N;
   GEN p1,s,v,c,tab;
 
   tab = get_tab(nf, &N);
+  if (typ(x) != t_COL || lg(x) != N+1
+   || typ(y) != t_COL || lg(y) != N+1) err(typeer,"element_muli");
   v=cgetg(N+1,t_COL); av=avma;
   for (k=1; k<=N; k++)
   {
@@ -249,7 +269,7 @@ element_muli(GEN nf, GEN x, GEN y)
 	}
       }
     }
-    v[k]=(long) gerepileuptoint(av,s); av=avma;
+    v[k] = lpileuptoint(av,s); av = avma;
   }
   return v;
 }
@@ -259,7 +279,8 @@ GEN
 element_sqri(GEN nf, GEN x)
 {
   GEN p1,s,v,c,tab;
-  long av,i,j,k,N;
+  gpmem_t av;
+  long i,j,k,N;
 
   tab = get_tab(nf, &N);
 
@@ -290,23 +311,19 @@ element_sqri(GEN nf, GEN x)
 	}
       }
     }
-    v[k]=(long) gerepileuptoint(av,s); av=avma;
+    v[k] = lpileuptoint(av,s); av = avma;
   }
   return v;
 }
 
-/* square of x in nf */
 GEN
-element_sqr(GEN nf, GEN x)
+sqr_by_tab(GEN tab, GEN x)
 {
-  long av = avma,i,j,k,N, tx = typ(x);
-  GEN p1,s,v,c,tab;
+  gpmem_t av = avma;
+  long i,j,k,N;
+  GEN p1,s,v,c;
 
-  nf = checknf(nf);
-  tab = get_tab(nf, &N);
-  if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_sqr");
-  if (is_extscalar_t(tx))
-    return gerepileupto(av, algtobasis(nf, gsqr(x)));
+  N = lg(x)-1;
   if (isnfscalar(x))
   {
     s=cgetg(N+1,t_COL); s[1]=lsqr((GEN)x[1]);
@@ -335,7 +352,8 @@ element_sqr(GEN nf, GEN x)
 	if (signe(c))
 	{
           p1 = gmul((GEN)x[i],(GEN)x[j]);
-	  p1 = gcmp1(c)? gmul2n(p1,1): gmul(p1,shifti(c,1));
+	  p1 = gmul2n(p1,1);
+          if (!gcmp1(c)) p1 = gmul(p1, c);
 	  s = gadd(s,p1);
 	}
       }
@@ -345,113 +363,115 @@ element_sqr(GEN nf, GEN x)
   return v;
 }
 
+/* square of x in nf */
+GEN
+element_sqr(GEN nf, GEN x)
+{
+  long N, tx = typ(x);
+  GEN tab;
+
+  nf = checknf(nf);
+  if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_sqr");
+  if (is_extscalar_t(tx))
+  {
+    gpmem_t av = avma;
+    return gerepileupto(av, algtobasis(nf, gsqr(x)));
+  }
+  if (tx != t_COL) err(typeer,"element_sqr");
+
+  tab = get_tab(nf, &N);
+  return sqr_by_tab(tab,x);
+}
+
+static GEN
+_mul(void *data, GEN x, GEN y) { return element_mul((GEN)data,x,y); }
+static GEN
+_sqr(void *data, GEN x) { return element_sqr((GEN)data,x); }
+
 /* Compute x^n in nf, left-shift binary powering */
 GEN
 element_pow(GEN nf, GEN x, GEN n)
 {
-  ulong av = avma;
-  long s,N,i,j,m;
-  GEN y,p1,cx;
+  gpmem_t av = avma;
+  long s,N;
+  GEN y, cx;
 
   if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
   nf=checknf(nf); N=degpol(nf[1]);
   s=signe(n); if (!s) return gscalcol_i(gun,N);
-  if (typ(x)!=t_COL) x=algtobasis(nf,x);
+  if (typ(x) != t_COL)
+  {
+    x = algtobasis(nf,x);
+    if (typ(x) != t_COL) err(typeer,"element_pow");
+  }
 
   if (isnfscalar(x))
   {
     y = gscalcol_i(gun,N);
     y[1] = (long)powgi((GEN)x[1],n); return y;
   }
-  cx = content(x);
-  if (gcmp1(cx)) cx = NULL; else x = gdiv(x,cx);
-  p1 = n+2; m = *p1;
-  y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; j; m<<=1,j--)
-    {
-      y = element_sqr(nf, y);
-      if (m<0) y=element_mul(nf, y, x);
-    }
-    if (--i == 0) break;
-    m = *++p1, j = BITS_IN_LONG;
-  }
+  x = primitive_part(x, &cx);
+  y = leftright_pow(x, n, (void*)nf, _sqr, _mul);
   if (s<0) y = element_inv(nf, y);
   if (cx) y = gmul(y, powgi(cx, n));
   return av==avma? gcopy(y): gerepileupto(av,y);
 }
 
-/* x in Z^n, compute lift(x^n mod p) */
-GEN
-element_pow_mod_p(GEN nf, GEN x, GEN n, GEN p)
-{
-  ulong av = avma;
-  long s,N,i,j,m;
-  GEN y,p1;
+typedef struct {
+  GEN nf, p;
+  long I;
+} eltmod_muldata;
 
-  if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
-  nf=checknf(nf); N=degpol(nf[1]);
-  s=signe(n); if (!s) return gscalcol_i(gun,N);
-  if (typ(x)!=t_COL) x=algtobasis(nf,x);
-
-  if (isnfscalar(x))
-  {
-    y = gscalcol_i(gun,N);
-    y[1] = (long)powmodulo((GEN)x[1],n,p); return y;
-  }
-  p1 = n+2; m = *p1;
-  y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; j; m<<=1,j--)
-    {
-      y = element_sqri(nf, y);
-      if (m<0) y=element_muli(nf, y, x);
-      y = FpV_red(y, p);
-    }
-    if (--i == 0) break;
-    m = *++p1, j = BITS_IN_LONG;
-  }
-  if (s<0)  y = FpV_red(element_inv(nf,y), p);
-  return av==avma? gcopy(y): gerepileupto(av,y);
+static GEN
+_mulidmod(void *data, GEN x, GEN y)
+{
+  eltmod_muldata *D = (eltmod_muldata*)data;
+  (void)y; /* ignored */
+  return FpV_red(element_mulid(D->nf, x, D->I), D->p);
 }
+static GEN
+_sqrmod(void *data, GEN x)
+{
+  eltmod_muldata *D = (eltmod_muldata*)data;
+  return FpV_red(element_sqri(D->nf, x), D->p);
+}
 
 /* x = I-th vector of the Z-basis of Z_K, in Z^n, compute lift(x^n mod p) */
 GEN
 element_powid_mod_p(GEN nf, long I, GEN n, GEN p)
 {
-  ulong av = avma;
-  long s,N,i,j,m;
-  GEN y,p1;
+  gpmem_t av = avma;
+  eltmod_muldata D;
+  long s,N;
+  GEN y;
 
   if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
   nf=checknf(nf); N=degpol(nf[1]);
   s=signe(n);
   if (!s || I == 1) return gscalcol_i(gun,N);
-  p1 = n+2; m = *p1;
   y = zerocol(N); y[I] = un;
-  j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
-  for (i=lgefint(n)-2;;)
-  {
-    for (; j; m<<=1,j--)
-    {
-      y = element_sqri(nf, y);
-      if (m<0) y=element_mulid(nf, y, I);
-      y = FpV_red(y, p);
-    }
-    if (--i == 0) break;
-    m = *++p1, j = BITS_IN_LONG;
-  }
-  if (s<0)  y = FpV_red(element_inv(nf,y), p);
+  D.nf = nf;
+  D.p = p;
+  D.I = I;
+  y = leftright_pow(y,n, (void*)&D, &_sqrmod, &_mulidmod);
+  if (s < 0) y = FpV_red(element_inv(nf,y), p);
   return av==avma? gcopy(y): gerepileupto(av,y);
 }
 
+GEN
+element_mulidid(GEN nf, long i, long j)
+{
+  long N;
+  GEN tab = get_tab(nf, &N);
+  tab += (i-1)*N; return (GEN)tab[j];
+}
+
 /* Outputs x.w_i, where w_i is the i-th elt of the integral basis */
 GEN
 element_mulid(GEN nf, GEN x, long i)
 {
-  long av,j,k,N;
+  gpmem_t av;
+  long j,k,N;
   GEN s,v,c,p1, tab;
 
   if (i==1) return gcopy(x);
@@ -473,25 +493,46 @@ element_mulid(GEN nf, GEN x, long i)
   return v;
 }
 
+/* table of multiplication by wi in ZK = Z[w1,..., wN] */
+GEN
+eltmulid_get_table(GEN nf, long i)
+{
+  long k,N;
+  GEN m, tab;
+
+  tab = get_tab(nf, &N);
+  tab += (i-1)*N;
+  m = cgetg(N+1,t_COL);
+  for (k=1; k<=N; k++) m[k] = tab[k];
+  return m;
+}
+
+/* table of multiplication by x in ZK = Z[w1,..., wN] */
+GEN
+eltmul_get_table(GEN nf, GEN x)
+{
+  long i, N = degpol(nf[1]);
+  GEN mul = cgetg(N+1,t_MAT);
+  if (typ(x) != t_COL) x = algtobasis(nf,x);
+  mul[1] = (long)x; /* assume w_1 = 1 */
+  for (i=2; i<=N; i++) mul[i] = (long)element_mulid(nf,x,i);
+  return mul;
+}
+
 /* valuation of integer x, with resp. to prime ideal P above p.
  * p.P^(-1) = b Z_K, v >= val_p(norm(x)), and N = deg(nf)
  * [b may be given as the 'multiplication by b' matrix]
  */
 long
-int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx, long v)
+int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx)
 {
   long i,k,w, N = degpol(nf[1]);
   GEN r,a,y,mul;
   
-  if (typ(b) == t_MAT) mul = b;
-  else
-  {
-    mul = cgetg(N+1,t_MAT);
-    for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
-  }
+  mul = (typ(b) == t_MAT)? b: eltmul_get_table(nf, b);
   y = cgetg(N+1, t_COL); /* will hold the new x */
   x = dummycopy(x);
-  for(w=0; w<=v; w++)
+  for(w=0;; w++)
   {
     for (i=1; i<=N; i++)
     { /* compute (x.b)_i */
@@ -499,19 +540,21 @@ int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx, lo
       for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
       /* is it divisible by p ? */
       y[i] = ldvmdii(a,p,&r);
-      if (signe(r)) goto END;
+      if (signe(r))
+      {
+        if (newx) *newx = x;
+        return w;
+      }
     }
     r=x; x=y; y=r;
   }
-END:
-  if (newx) *newx = x;
-  return w;
 }
 
 long
 element_val(GEN nf, GEN x, GEN vp)
 {
-  long av = avma,N,w,vcx,e;
+  gpmem_t av = avma;
+  long N,w,vcx,e;
   GEN cx,p;
 
   if (gcmp0(x)) return VERYBIGINT;
@@ -523,7 +566,7 @@ element_val(GEN nf, GEN x, GEN vp)
       return ggval(x,p)*e;
     case t_POLMOD: x = (GEN)x[2]; /* fall through */
     case t_POL:
-      x = algtobasis_intern(nf,x); break;
+      x = algtobasis_i(nf,x); break;
     case t_COL:
       if (lg(x)==N+1) break;
     default: err(typeer,"element_val");
@@ -532,23 +575,10 @@ element_val(GEN nf, GEN x, GEN vp)
 
   cx = content(x);
   if (gcmp1(cx)) vcx=0; else { x = gdiv(x,cx); vcx = ggval(cx,p); }
-  w = int_elt_val(nf,x,p,(GEN)vp[5],NULL,VERYBIGINT);
+  w = int_elt_val(nf,x,p,(GEN)vp[5],NULL);
   avma=av; return w + vcx*e;
 }
 
-/* d = a multiple of norm(x), x integral */
-long
-element_val2(GEN nf, GEN x, GEN d, GEN vp)
-{
-  GEN p = (GEN)vp[1];
-  long av, v = ggval(d,p);
-
-  if (!v) return 0;
-  av=avma;
-  v = int_elt_val(nf,x,p,(GEN)vp[5],NULL,v);
-  avma=av; return v;
-}
-
 /* polegal without comparing variables */
 long
 polegal_spec(GEN x, GEN y)
@@ -564,7 +594,7 @@ polegal_spec(GEN x, GEN y)
 GEN
 basistoalg(GEN nf, GEN x)
 {
-  long tx=typ(x),lx=lg(x),i;
+  long tx=typ(x),lx=lg(x),i,j,l;
   GEN z;
 
   nf=checknf(nf);
@@ -584,9 +614,18 @@ basistoalg(GEN nf, GEN x)
       }
       /* fall through */
 
-    case t_VEC: case t_MAT: z=cgetg(lx,tx);
+    case t_VEC: z=cgetg(lx,tx);
       for (i=1; i<lx; i++) z[i]=(long)basistoalg(nf,(GEN)x[i]);
       return z;
+    case t_MAT: z=cgetg(lx,t_MAT);
+      if (lx == 1) return z;
+      l = lg(x[1]);
+      for (j=1; j<lx; j++)
+      {
+        z[j] = lgetg(l,t_COL);
+        for (i=1; i<l; i++) coeff(z,i,j) = (long)basistoalg(nf,gcoeff(x,i,j));
+      }
+      return z;
 
     case t_POLMOD:
       if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
@@ -598,32 +637,14 @@ basistoalg(GEN nf, GEN x)
   }
 }
 
-/* return the (N-dimensional) vector of coeffs of p */
+/* gmul(A, pol_to_vec(x)), A t_MAT (or t_VEC) of compatible dimensions */
 GEN
-pol_to_vec(GEN x, long N)
-{
-  long i, l;
-  GEN z = cgetg(N+1,t_COL);
-  if (typ(x) != t_POL)
-  {
-    z[1] = (long)x;
-    for (i=2; i<=N; i++) z[i]=zero;
-    return z;
-  }
-  l = lgef(x)-1; x++;
-  for (i=1; i<l ; i++) z[i]=x[i];
-  for (   ; i<=N; i++) z[i]=zero;
-  return z;
-}
-
-/* gmul(A, pol_to_vec(x)), A matrix of compatible dimensions */
-GEN
 mulmat_pol(GEN A, GEN x)
 {
   long i,l;
   GEN z;
-  if (typ(x) != t_POL) return gscalcol(x, lg(A[1])-1);
-  l=lgef(x)-1; if (l == 1) return zerocol(lg(A[1])-1);
+  if (typ(x) != t_POL) return gmul(x,(GEN)A[1]); /* scalar */
+  l=lgef(x)-1; if (l == 1) return typ(A)==t_VEC? gzero: zerocol(lg(A[1])-1);
   x++; z = gmul((GEN)x[1], (GEN)A[1]);
   for (i=2; i<l ; i++) 
     if (!gcmp0((GEN)x[i])) z = gadd(z, gmul((GEN)x[i], (GEN)A[i]));
@@ -631,10 +652,9 @@ mulmat_pol(GEN A, GEN x)
 }
 
 /* valid for scalars and polynomial, degree less than N.
- * No garbage collecting. No check (SEGV for vectors).
- */
+ * No garbage collecting. No check.  */
 GEN
-algtobasis_intern(GEN nf, GEN x)
+algtobasis_i(GEN nf, GEN x)
 {
   GEN P = (GEN)nf[1];
   long tx = typ(x), N = degpol(P);
@@ -653,7 +673,8 @@ algtobasis_intern(GEN nf, GEN x)
 GEN
 algtobasis(GEN nf, GEN x)
 {
-  long tx=typ(x),lx=lg(x),av=avma,i,N;
+  long tx=typ(x),lx=lg(x),i,N;
+  gpmem_t av=avma;
   GEN z;
 
   nf=checknf(nf);
@@ -668,7 +689,7 @@ algtobasis(GEN nf, GEN x)
 	err(talker,"not the same number field in algtobasis");
       x = (GEN)x[2]; /* fall through */
     case t_POL:
-      return gerepileupto(av,algtobasis_intern(nf,x));
+      return gerepileupto(av,algtobasis_i(nf,x));
 
     default: N=degpol(nf[1]); return gscalcol(x,N);
   }
@@ -680,7 +701,7 @@ algtobasis(GEN nf, GEN x)
 GEN
 nfdiveuc(GEN nf, GEN a, GEN b)
 {
-  long av=avma, tetpil;
+  gpmem_t av=avma, tetpil;
   a = element_div(nf,a,b); tetpil=avma;
   return gerepile(av,tetpil,ground(a));
 }
@@ -691,7 +712,7 @@ nfdiveuc(GEN nf, GEN a, GEN b)
 GEN
 nfmod(GEN nf, GEN a, GEN b)
 {
-  long av=avma,tetpil;
+  gpmem_t av=avma,tetpil;
   GEN p1=gneg_i(element_mul(nf,b,ground(element_div(nf,a,b))));
   tetpil=avma; return gerepile(av,tetpil,gadd(a,p1));
 }
@@ -702,7 +723,7 @@ nfmod(GEN nf, GEN a, GEN b)
 GEN
 nfdivres(GEN nf, GEN a, GEN b)
 {
-  long av=avma,tetpil;
+  gpmem_t av=avma,tetpil;
   GEN p1,z, y = ground(element_div(nf,a,b));
 
   p1=gneg_i(element_mul(nf,b,y)); tetpil=avma;
@@ -715,62 +736,94 @@ nfdivres(GEN nf, GEN a, GEN b)
 /**			      (Z_K/I)^*					**/
 /**									**/
 /*************************************************************************/
+/* return sign(sigma_k(x)), x t_COL (integral, primitive) */
+static long
+eval_sign(GEN M, GEN x, long k)
+{
+  long i, l = lg(x);
+  GEN z = mpmul(gcoeff(M,k,1), (GEN)x[1]);
+  for (i = 2; i < l; i++)
+    z = mpadd(z, mpmul(gcoeff(M,k,i), (GEN)x[i]));
+  if (lg(z) < DEFAULTPREC) err(precer,"zsigne");
+  return signe(z);
+}
 
+static void
+const_sign(GEN v, GEN arch, GEN s)
+{
+  long i,j, l = lg(v);
+  for (j=i=1; i<l; i++)
+    if (signe(arch[i])) v[j++] = (long)s;
+  setlg(v, j);
+}
+
 /* return (column) vector of R1 signatures of x (coeff modulo 2)
- * if arch = NULL, assume arch = [0,..0]
- */
+ * if arch = NULL, assume arch = [0,..0] */
 GEN
 zsigne(GEN nf,GEN x,GEN arch)
 {
-  GEN _0,_1, vecsign, rac = (GEN)nf[6];
-  long i,j,l,e,B;
-  ulong av0, av;
+  GEN _0, _1, M, vecsign;
+  long i,j,l,s;
+  gpmem_t av0, av;
 
   if (!arch) return cgetg(1,t_COL);
   av0 = avma;
   l = lg(arch); vecsign = cgetg(l,t_COL);
   _0 = gmodulss(0,2);
   _1 = gmodulss(1,2); av = avma;
+  nf = checknf(nf);
+  M = gmael(nf,5,1);
   switch(typ(x))
   {
     case t_MAT: /* factorisation */
     {
-      GEN t = (GEN)x[1], e = (GEN)x[2];
-      for (j=i=1; i<l; i++)
-        if (signe(arch[i])) vecsign[j++] = (long)_0;
-      setlg(vecsign,j);
-      if (j==1) { avma = av0; return cgetg(1, t_COL); }
-      for (i=1; i<lg(t); i++)
-      {
-        GEN p1 = (GEN)e[i];
-        if (mpodd(p1))
-          vecsign = gadd(vecsign, zsigne(nf,(GEN)t[i],arch));
-      }
+      GEN g = (GEN)x[1], e = (GEN)x[2];
+      const_sign(vecsign, arch, _0);
+      if (lg(vecsign)==1) { avma = av0; return cgetg(1, t_COL); }
+      for (i=1; i<lg(g); i++)
+        if (mpodd((GEN)e[i]))
+          vecsign = gadd(vecsign, zsigne(nf,(GEN)g[i],arch));
       return gerepileupto(av0, vecsign);
     }
-    case t_COL: x = gmul((GEN)nf[7],x); break;
-    case t_POLMOD: x = (GEN)x[2];
+    case t_POLMOD: x = (GEN)x[2];      /* fall through */
+    case t_POL: x = algtobasis(nf, x); /* fall through */
+    case t_COL: if (!isnfscalar(x)) break;
+                x = (GEN)x[1];         /* fall through */
+    case t_INT: case t_FRAC:
+    {
+      s = gsigne(x); if (!s) err(talker,"zero element in zsigne");
+      const_sign(vecsign, arch, (s > 0)? _0: _1);
+      return vecsign;
+    }
   }
-  if (gcmp0(x)) err(talker,"zero element in zsigne");
-
-  B = bit_accuracy(precision((GEN)rac[1]));
-  e = gexpo(x);
+  x = Q_primpart(x);
   for (j=i=1; i<l; i++)
     if (signe(arch[i]))
-    {
-      GEN y = poleval(x,(GEN)rac[i]);
-      if (e + gexpo((GEN)rac[i]) - gexpo(y) > B)
-        err(precer, "zsigne");
-      vecsign[j++] = (gsigne(y) > 0)? (long)_0: (long)_1;
-    }
+      vecsign[j++] = (eval_sign(M, x, i) > 0)? (long)_0: (long)_1;
   avma = av; setlg(vecsign,j); return vecsign;
 }
 
+/* return the t_COL vector of signs of x; the matrix of such if x is a vector
+ * of nf elements */
 GEN
+zsigns(GEN nf, GEN x)
+{
+  long r1, i, l;
+  GEN arch, S;
+
+  nf = checknf(nf); r1 = nf_get_r1(nf);
+  arch = cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i] = un;
+  if (typ(x) != t_VEC) return zsigne(nf, x, arch);
+  l = lg(x); S = cgetg(l, t_MAT);
+  for (i=1; i<l; i++) S[i] = (long)zsigne(nf, (GEN)x[i], arch);
+  return S;
+}
+
+GEN
 lllreducemodmatrix(GEN x,GEN y)
 {
-  ulong av = avma;
-  GEN z = gmul(y,lllint(y));
+  gpmem_t av = avma;
+  GEN z = lllint_ip(y,4);
   return gerepileupto(av, reducemodinvertible(x, z));
 }
 
@@ -783,7 +836,7 @@ reducemodinvertible(GEN x, GEN y)
 
 /* Reduce column x modulo y in HNF */
 GEN
-colreducemodmat(GEN x, GEN y, GEN *Q)
+colreducemodHNF(GEN x, GEN y, GEN *Q)
 {
   long i, l = lg(x);
   GEN q;
@@ -792,7 +845,7 @@ colreducemodmat(GEN x, GEN y, GEN *Q)
   if (l == 1) return cgetg(1,t_COL);
   for (i = l-1; i>0; i--)
   {
-    q = negi(gdivround((GEN)x[i], gcoeff(y,i,i)));
+    q = negi(diviiround((GEN)x[i], gcoeff(y,i,i)));
     if (Q) (*Q)[i] = (long)q;
     if (signe(q)) x = gadd(x, gmul(q, (GEN)y[i]));
   }
@@ -815,26 +868,31 @@ reducemodHNF(GEN x, GEN y, GEN *Q)
   if (Q)
   {
     GEN q = cgetg(lx, t_MAT); *Q = q;
-    for (i=1; i<lx; i++) R[i] = (long)colreducemodmat((GEN)x[i],y,(GEN*)(q+i));
+    for (i=1; i<lx; i++) R[i] = (long)colreducemodHNF((GEN)x[i],y,(GEN*)(q+i));
   }
   else
     for (i=1; i<lx; i++)
     {
-      long av = avma;
-      R[i] = lpileupto(av, colreducemodmat((GEN)x[i],y,NULL));
+      gpmem_t av = avma;
+      R[i] = lpileupto(av, colreducemodHNF((GEN)x[i],y,NULL));
     }
   return R;
 }
 
-/* For internal use. Reduce x modulo ideal. We want a non-zero result */
+/* For internal use. Reduce x modulo ideal (assumed non-zero, in HNF). We
+ * want a non-zero result */
 GEN
-nfreducemodideal(GEN nf,GEN x0,GEN ideal)
+nfreducemodideal_i(GEN x0,GEN ideal)
 {
-  GEN y = idealhermite(nf,ideal);
-  GEN x = colreducemodmat(x0, y, NULL);
-  if (gcmp0(x)) return (GEN)y[1];
+  GEN x = colreducemodHNF(x0, ideal, NULL);
+  if (gcmp0(x)) return (GEN)ideal[1];
   return x == x0? gcopy(x) : x;
 }
+GEN
+nfreducemodideal(GEN nf,GEN x0,GEN ideal)
+{
+  return nfreducemodideal_i(x0, idealhermite(nf,ideal));
+}
 
 /* multiply y by t = 1 mod^* f such that sign(x) = sign(y) at arch = idele[2].
  * If x == NULL, make y >> 0 at sarch */
@@ -902,12 +960,11 @@ element_powmodidele(GEN nf,GEN x,GEN k,GEN idele,GEN s
 GEN
 smithrel(GEN H, GEN *newU, GEN *newUi)
 { 
-  GEN z,U,Ui,D,p1;
+  GEN U,Ui,D,p1;
   long l,c;
 
-  z = smith2(H);
-  U = (GEN)z[1];
-  D = (GEN)z[3]; l = lg(D);
+  D = smithall(H, &U, NULL);
+  l = lg(D);
   for (c=1; c<l; c++)
   {
     p1 = gcoeff(D,c,c);
@@ -943,51 +1000,13 @@ zidealij(GEN x, GEN y)
   z[3] = lmul(U,xi); return z;
 }
 
-static GEN
-get_multab(GEN nf, GEN x)
-{
-  long lx = lg(x), i;
-  GEN multab = cgetg(lx, t_MAT);
-  for (i=1; i<lx; i++)
-    multab[i]=(long)element_mulid(nf,x,i);
-  return multab;
-}
-
-/* return mat * y mod prh */
-static GEN
-mul_matvec_mod_pr(GEN mat, GEN y, GEN prh)
-{
-  long av,i,j, lx = lg(mat);
-  GEN v, res = cgetg(lx,t_COL), p = gcoeff(prh,1,1);
-
-  av = avma; (void)new_chunk((lx-1) * lgefint(p));
-  v = zerocol(lx-1);
-  for (i=lx-1; i; i--)
-  {
-    GEN p1 = (GEN)v[i], t = (GEN)prh[i];
-    for (j=1; j<lx; j++)
-      p1 = addii(p1, mulii(gcoeff(mat,i,j), (GEN)y[j]));
-    p1 = modii(p1, p);
-    if (p1 != gzero && is_pm1(t[i]))
-    {
-      for (j=1; j<i; j++)
-        v[j] = lsubii((GEN)v[j], mulii(p1, (GEN)t[j]));
-      p1 = gzero;
-    }
-    if (p1 == gzero) /* intended */
-      res[i] = zero;
-    else
-      av = res[i] = (long)icopy_av(p1, (GEN)av);
-  }
-  avma = av; return res;
-}
-
 /* smallest integer n such that g0^n=x modulo p prime. Assume g0 has order q
  * (p-1 if q = NULL) */
 GEN
 Fp_shanks(GEN x,GEN g0,GEN p, GEN q)
 {
-  long av=avma,av1,lim,lbaby,i,k,c;
+  gpmem_t av=avma,av1,lim;
+  long lbaby,i,k,c;
   GEN p1,smalltable,giant,perm,v,g0inv;
 
   x = modii(x,p);
@@ -1005,7 +1024,7 @@ Fp_shanks(GEN x,GEN g0,GEN p, GEN q)
     av1 = avma;
     if (is_pm1(p1)) { avma=av; return stoi(i-1); }
     smalltable[i]=(long)p1; if (i==lbaby) break;
-    new_chunk(c); p1 = mulii(p1,g0inv);
+    (void)new_chunk(c); p1 = mulii(p1,g0inv); /* HACK */
     avma = av1; p1 = resii(p1, p);
   }
   giant = resii(mulii(x, mpinvmod(p1,p)), p);
@@ -1037,7 +1056,7 @@ Fp_shanks(GEN x,GEN g0,GEN p, GEN q)
 GEN
 Fp_PHlog(GEN a, GEN g, GEN p, GEN ord)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv;
   GEN fa, ex;
   long e,i,j,l;
@@ -1081,120 +1100,123 @@ Fp_PHlog(GEN a, GEN g, GEN p, GEN ord)
   return gerepileuptoint(av, lift(chinese(v,NULL)));
 }
 
+/* discrete log in Fq for a in Fp^*, g primitive root in Fq^* */
+static GEN
+ff_PHlog_Fp(GEN a, GEN g, GEN T, GEN p)
+{
+  gpmem_t av = avma;
+  GEN q,n_q,ord,ordp;
+
+  if (gcmp1(a) || egalii(p, gdeux)) { avma = av; return gzero; }
+
+  ordp = subis(p, 1);
+  ord = T? subis(gpowgs(p,degpol(T)), 1): p;
+  if (egalii(a, ordp)) /* -1 */
+    return gerepileuptoint(av, shifti(ord,-1));
+
+  if (!T) q = NULL;
+  else 
+  { /* we want < g > = Fp^* */
+    q = divii(ord,ordp);
+    g = FpXQ_pow(g,q,T,p);
+    if (typ(g) == t_POL) g = constant_term(g);
+  }
+  n_q = Fp_PHlog(a,g,p,NULL);
+  if (q) n_q = mulii(q, n_q);
+  return gerepileuptoint(av, n_q);
+}
+
 /* smallest n >= 0 such that g0^n=x modulo pr, assume g0 reduced
- * q = order of g0  (Npr - 1 if q = NULL)
- * TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
+ * q = order of g0 is prime (and != p) */
 static GEN
-nfshanks(GEN nf,GEN x,GEN g0,GEN pr,GEN prhall,GEN q)
+ffshanks(GEN x, GEN g0, GEN q, GEN T, GEN p)
 {
-  long av=avma,av1,lim,lbaby,i,k, f = itos((GEN)pr[4]);
-  GEN p1,smalltable,giant,perm,v,g0inv,prh = (GEN)prhall[1];
-  GEN multab, p = (GEN)pr[1];
+  gpmem_t av = avma, av1, lim;
+  long lbaby,i,k;
+  GEN p1,smalltable,giant,perm,v,g0inv;
 
-  x = lift_intern(nfreducemodpr(nf,x,prhall));
-  p1 = q? q: addsi(-1, gpowgs(p,f));
-  p1 = racine(p1);
-  if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in nfshanks");
-  lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
-  g0inv = lift_intern(element_invmodpr(nf,g0,prhall));
+  if (!degpol(x)) x = constant_term(x);
+  if (typ(x) == t_INT)
+  {
+    if (!gcmp1(modii(p,q))) return gzero;
+    /* g0 in Fp^*, order q | (p-1) */
+    if (typ(g0) == t_POL) g0 = constant_term(g0);
+    return Fp_PHlog(x,g0,p,q);
+  }
+
+  p1 = racine(q);
+  if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in ffshanks");
+  lbaby = itos(p1)+1; smalltable = cgetg(lbaby+1,t_VEC);
+  g0inv = FpXQ_inv(g0,T,p);
   p1 = x;
 
-  multab = get_multab(nf, g0inv);
-  for (i=lg(multab)-1; i; i--)
-    multab[i]=(long)FpV_red((GEN)multab[i], p);
-
   for (i=1;;i++)
   {
-    if (isnfscalar(p1) && gcmp1((GEN)p1[1])) { avma=av; return stoi(i-1); }
+    if (gcmp1(p1)) { avma = av; return stoi(i-1); }
 
     smalltable[i]=(long)p1; if (i==lbaby) break;
-    p1 = mul_matvec_mod_pr(multab,p1,prh);
+    p1 = FpXQ_mul(p1,g0inv, T,p);
   }
-  giant=lift_intern(element_divmodpr(nf,x,p1,prhall));
-  p1=cgetg(lbaby+1,t_VEC);
-  perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_vecint);
-  for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
-  smalltable=p1; p1=giant;
+  giant = FpXQ_div(x,p1,T,p);
+  perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_pol);
+  smalltable = vecextract_p(smalltable, perm);
+  p1 = giant;
 
-  multab = get_multab(nf, giant);
-  for (i=lg(multab)-1; i; i--)
-    multab[i]=(long)FpV_red((GEN)multab[i], p);
-
   av1 = avma; lim=stack_lim(av1,2);
   for (k=1;;k++)
   {
-    i=tablesearch(smalltable,p1,cmp_vecint);
+    i = tablesearch(smalltable,p1,cmp_pol);
     if (i)
     {
-      v=addis(mulss(lbaby-1,k),perm[i]);
-      return gerepileuptoint(av,addsi(-1,v));
+      v = addis(mulss(lbaby-1,k), perm[i]);
+      return gerepileuptoint(av, addsi(-1,v));
     }
-    p1 = mul_matvec_mod_pr(multab,p1,prh);
+    p1 = FpXQ_mul(p1, giant, T,p);
 
     if (low_stack(lim, stack_lim(av1,2)))
     {
-      if(DEBUGMEM>1) err(warnmem,"nfshanks");
+      if(DEBUGMEM>1) err(warnmem,"ffshanks");
       p1 = gerepileupto(av1, p1);
     }
   }
 }
 
-/* same in nf.zk / pr
- * TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
+/* same in Fp[X]/T */
 GEN
-nf_PHlog(GEN nf, GEN a, GEN g, GEN pr, GEN prhall)
+ff_PHlog(GEN a, GEN g, GEN T, GEN p)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv,ord,fa,ex;
   long e,i,j,l;
 
-  a = lift_intern(nfreducemodpr(nf,a,prhall));
-  if (isnfscalar(a))
-  { /* can be done in Fp^* */
-    GEN p = (GEN)pr[1], ordp = subis(p, 1);
-    a = (GEN)a[1];
-    if (gcmp1(a) || egalii(p, gdeux)) { avma = av; return gzero; }
-    ord = subis(idealnorm(nf,pr), 1);
-    if (egalii(a, ordp)) /* -1 */
-      return gerepileuptoint(av, shifti(ord,-1));
-
-    if (egalii(ord, ordp))
-      q = NULL;
-    else /* we want < g > = Fp^* */
-    {
-      q = divii(ord,ordp);
-      g = element_powmodpr(nf,g,q,prhall);
-    }
-    g = lift_intern((GEN)g[1]);
-    n_q = Fp_PHlog(a,g,p,NULL);
-    if (q) n_q = mulii(q, n_q);
-    return gerepileuptoint(av, n_q);
-  }
-  ord = subis(idealnorm(nf,pr),1);
-  fa = factor(ord); ex = (GEN)fa[2];
+  if (typ(a) == t_INT)
+    return gerepileuptoint(av, ff_PHlog_Fp(a,g,T,p));
+  /* f > 1 ==> T != NULL */
+  ord = subis(gpowgs(p, degpol(T)), 1);
+  fa = factor(ord);
+  ex = (GEN)fa[2];
   fa = (GEN)fa[1];
   l = lg(fa);
-  ginv = lift_intern(element_invmodpr(nf, g, prhall));
+  ginv = FpXQ_inv(g,T,p);
   v = cgetg(l, t_VEC);
   for (i=1; i<l; i++)
   {
     q = (GEN)fa[i];
     e = itos((GEN)ex[i]);
-    if (DEBUGLEVEL>5) 
-      fprintferr("nf_Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
+    if (DEBUGLEVEL>5) fprintferr("nf_Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
     qj = new_chunk(e+1); qj[0] = un;
     for (j=1; j<=e; j++) qj[j] = lmulii((GEN)qj[j-1], q);
     t0 = divii(ord, (GEN)qj[e]);
-    a0 = element_powmodpr(nf, a, t0, prhall); 
-    ginv0 = element_powmodpr(nf, ginv, t0, prhall); /* order q^e */
-    g_q = element_powmodpr(nf, g, divii(ord,q), prhall); /* order q */
+    a0 = FpXQ_pow(a, t0, T,p); 
+    ginv0 = FpXQ_pow(ginv, t0, T,p); /* order q^e */
+    g_q = FpXQ_pow(g, divii(ord,q), T,p); /* order q */
+
     n_q = gzero;
     for (j=0; j<e; j++)
     {
-      b = element_mulmodpr(nf,a0,
-                              element_powmodpr(nf, ginv0, n_q, prhall), prhall);
-      b = element_powmodpr(nf, b, (GEN)qj[e-1-j], prhall);
-      b = nfshanks(nf,b, g_q, pr,prhall, q);
+      b = FpXQ_mul(a0, FpXQ_pow(ginv0, n_q, T,p), T,p);
+      b = FpXQ_pow(b, (GEN)qj[e-1-j], T,p);
+      b = ffshanks(b, g_q, q, T,p);
       n_q = addii(n_q, mulii(b, (GEN)qj[j]));
     }
     v[i] = lmodulcp(n_q, (GEN)qj[e]);
@@ -1202,10 +1224,21 @@ nf_PHlog(GEN nf, GEN a, GEN g, GEN pr, GEN prhall)
   return gerepileuptoint(av, lift(chinese(v,NULL)));
 }
 
+/* same in nf.zk / pr */
 GEN
+nf_PHlog(GEN nf, GEN a, GEN g, GEN pr)
+{
+  gpmem_t av = avma;
+  GEN T,p, modpr = nf_to_ff_init(nf, &pr, &T, &p);
+  GEN A = nf_to_ff(nf,a,modpr);
+  GEN G = nf_to_ff(nf,g,modpr);
+  return gerepileuptoint(av, ff_PHlog(A,G,T,p));
+}
+
+GEN
 znlog(GEN x, GEN g0)
 {
-  long av=avma;
+  gpmem_t av=avma;
   GEN p = (GEN)g0[1];
   if (typ(g0) != t_INTMOD)
     err(talker,"not an element of (Z/pZ)* in znlog");
@@ -1224,7 +1257,7 @@ GEN
 dethnf(GEN mat)
 {
   long i,l = lg(mat);
-  ulong av;
+  gpmem_t av;
   GEN s;
 
   if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
@@ -1236,7 +1269,8 @@ dethnf(GEN mat)
 GEN
 dethnf_i(GEN mat)
 {
-  long av,i,l = lg(mat);
+  gpmem_t av;
+  long i,l = lg(mat);
   GEN s;
 
   if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
@@ -1249,7 +1283,8 @@ dethnf_i(GEN mat)
 GEN
 detcyc(GEN cyc)
 {
-  long av,i,l = lg(cyc);
+  gpmem_t av;
+  long i,l = lg(cyc);
   GEN s;
 
   if (l<3) return l<2? gun: icopy((GEN)cyc[1]);
@@ -1259,11 +1294,11 @@ detcyc(GEN cyc)
 }
 
 /* (U,V) = 1. Return y = x mod U, = 1 mod V (uv: u + v = 1, u in U, v in V) */
-static GEN
+GEN
 makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x)
 {
   GEN y = gadd((GEN)uv[1], element_mul(nf,x,(GEN)uv[2]));
-  return nfreducemodideal(nf,y,UV);
+  return nfreducemodideal_i(y,UV);
 }
 
 static GEN
@@ -1276,6 +1311,27 @@ makeprimetoidealvec(GEN nf,GEN UV,GEN uv,GEN gen)
   return y;
 }
 
+GEN
+FpXQ_gener(GEN T, GEN p)
+{
+  long i,j, k, vT = varn(T), f = degpol(T);
+  GEN g, list, pf_1 = subis(gpowgs(p, f), 1);
+  gpmem_t av0 = avma, av;
+
+  list = (GEN)factor(pf_1)[1];
+  k = lg(list)-1;
+
+  for (i=1; i<=k; i++) list[i] = (long)diviiexact(pf_1, (GEN)list[i]);
+  for (av = avma;; avma = av)
+  {
+    g = FpX_rand(f, vT, p);
+    if (degpol(g) < 1) continue;
+    for (j=1; j<=k; j++)
+      if (gcmp1(FpXQ_pow(g, (GEN)list[j], T, p))) break;
+    if (j > k) return gerepilecopy(av0, g);
+  }
+}
+
 /* Given an ideal pr^ep, and an integral ideal x (in HNF form) compute a list
  * of vectors, each with 5 components as follows :
  * [[clh],[gen1],[gen2],[signat2],U.X^-1]. Each component corresponds to
@@ -1287,51 +1343,39 @@ makeprimetoidealvec(GEN nf,GEN UV,GEN uv,GEN gen)
 static GEN
 zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
 {
-  long av=avma,av1,N,f,nbp,j,n,m,tetpil,i,e,a,b;
+  gpmem_t av = avma, av1, tetpil;
+  long N, f, i, e, a, b;
   GEN prh,p,pf_1,list,v,p1,p3,p4,prk,uv,g0,newgen,pra,prb;
   GEN *gptr[2];
 
   if(DEBUGLEVEL>3) { fprintferr("treating pr = %Z ^ %Z\n",pr,ep); flusherr(); }
-  prh=prime_to_ideal(nf,pr); N=lg(prh)-1;
-  f=itos((GEN)pr[4]); p=(GEN)pr[1];
+  prh = prime_to_ideal(nf,pr); N = degpol(nf[1]);
+  f = itos((GEN)pr[4]);
+  p = (GEN)pr[1];
 
   pf_1 = addis(gpowgs(p,f), -1);
-  v = zerocol(N);
-  if (f==1) v[1]=gener(p)[2];
+  if (f==1)
+  {
+    v = zerocol(N);
+    v[1] = gener(p)[2];
+  }
   else
   {
-    GEN prhall = cgetg(3,t_VEC);
-    long psim;
-    if (is_bigint(p)) err(talker,"prime too big in zprimestar");
-    psim = itos(p);
-    list = (GEN)factor(pf_1)[1]; nbp=lg(list)-1;
-    prhall[1]=(long)prh; prhall[2]=zero;
-    for (n=psim; n>=0; n++)
-    {
-      m=n;
-      for (i=1; i<=N; i++)
-	if (!gcmp1(gcoeff(prh,i,i))) { v[i]=lstoi(m%psim); m/=psim; }
-      for (j=1; j<=nbp; j++)
-      {
-        p1 = divii(pf_1,(GEN)list[j]);
-	p1 = lift_intern(element_powmodpr(nf,v,p1,prhall));
-        if (isnfscalar(p1) && gcmp1((GEN)p1[1])) break;
-      }
-      if (j>nbp) break;
-    }
-    if (n < 0) err(talker,"prime too big in zprimestar");
+    GEN T, modpr = zk_to_ff_init(nf, &pr, &T, &p);
+    v = ff_to_nf(FpXQ_gener(T,p), modpr);
+    v = algtobasis(nf, v);
   }
   /* v generates  (Z_K / pr)^* */
-  if(DEBUGLEVEL>3) {fprintferr("v computed\n");flusherr();}
+  if(DEBUGLEVEL>3) fprintferr("v computed\n");
   e = itos(ep); prk=(e==1)? pr: idealpow(nf,pr,ep);
-  if(DEBUGLEVEL>3) {fprintferr("prk computed\n");flusherr();}
+  if(DEBUGLEVEL>3) fprintferr("prk computed\n");
   g0 = v;
   uv = NULL; /* gcc -Wall */
   if (x)
   {
-    uv = idealaddtoone(nf,prk, idealdivexact(nf,x,prk));
+    uv = idealaddtoone(nf,prk, idealdivpowprime(nf,x,pr,ep));
     g0 = makeprimetoideal(nf,x,uv,v);
-    if(DEBUGLEVEL>3) {fprintferr("g0 computed\n");flusherr();}
+    if(DEBUGLEVEL>3) fprintferr("g0 computed\n");
   }
 
   p1 = cgetg(6,t_VEC); list = _vec(p1);
@@ -1346,12 +1390,11 @@ zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
   pra = prh; prb = (e==2)? prk: idealpow(nf,pr,gdeux);
   for(;;)
   {
-    if(DEBUGLEVEL>3)
-      {fprintferr("  treating a = %ld, b = %ld\n",a,b); flusherr();}
+    if(DEBUGLEVEL>3) fprintferr("  treating a = %ld, b = %ld\n",a,b);
     p1 = zidealij(pra,prb);
     newgen = dummycopy((GEN)p1[2]);
     p3 = cgetg(lg(newgen),t_VEC);
-    if(DEBUGLEVEL>3) {fprintferr("zidealij done\n"); flusherr();}
+    if(DEBUGLEVEL>3) fprintferr("zidealij done\n");
     for (i=1; i<lg(newgen); i++)
     {
       if (x) newgen[i]=(long)makeprimetoideal(nf,x,uv,(GEN)newgen[i]);
@@ -1371,17 +1414,24 @@ zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
     pra = gcopy(prb);
     gptr[0]=&pra; gptr[1]=&list;
     gerepilemanysp(av1,tetpil,gptr,2);
-    prb = (b==(a<<1))? idealpow(nf,pra,gdeux): prk;
+    prb = (b==(a<<1))? idealpows(nf,pr,b): prk;
   }
 }
 
-/* x ideal, compute elements in 1+x whose sign matrix is invertible */
+/* x integral ideal, compute elements in 1+x whose sign matrix is invertible */
 GEN
 zarchstar(GEN nf,GEN x,GEN arch,long nba)
 {
-  long av,N,i,j,k,r,rr,limr,zk,lgmat;
+  long limr, N, i, j, k, r, rr, zk, lgmat;
+  gpmem_t av;
   GEN p1,y,bas,genarch,mat,lambda,nfun,vun;
 
+  if (nba < 0)
+  {
+    nba = 0;
+    for (i=1; i<lg(arch); i++)
+      if (signe(arch[i])) nba++;
+  }
   y = cgetg(4,t_VEC);
   if (!nba)
   {
@@ -1398,10 +1448,10 @@ zarchstar(GEN nf,GEN x,GEN arch,long nba)
     y[2] = (long)_vec(_col(p1));
     y[3] = (long)gscalmat(gun,1); return y;
   }
-  zk = ideal_is_zk(x,N);
+  zk = gcmp1(gcoeff(x,1,1));
   genarch = cgetg(nba+1,t_VEC);
   mat = cgetg(nba+1,t_MAT); setlg(mat,2);
-  for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_MAT);
+  for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_VECSMALL);
   lambda = cgetg(N+1, t_VECSMALL);
 
   bas = dummycopy(gmael(nf,5,1)); r = lg(arch);
@@ -1415,7 +1465,7 @@ zarchstar(GEN nf,GEN x,GEN arch,long nba)
   r = nba+1; for (i=1; i<=N; i++) setlg(bas[i], r);
   if (!zk)
   {
-    x = gmul(x,lllintpartial(x));
+    x = lllint_ip(x,4);
     bas = gmul(bas, x);
   }
 
@@ -1426,12 +1476,13 @@ zarchstar(GEN nf,GEN x,GEN arch,long nba)
   {
     if (DEBUGLEVEL>5) fprintferr("zarchstar: r = %ld\n",r);
     p1 = gpowgs(stoi(rr),N);
-    limr = is_bigint(p1)? BIGINT: p1[2];
+    limr = (cmpis(p1,BIGINT) >= 0)? BIGINT: p1[2];
     limr = (limr-1)>>1;
     k = zk? rr: 1; /* to avoid 0 */
     for ( ; k<=limr; k++)
     {
-      long av1=avma, kk = k;
+      gpmem_t av1=avma;
+      long kk = k;
       GEN alpha = vun;
       for (i=1; i<=N; i++)
       {
@@ -1441,9 +1492,9 @@ zarchstar(GEN nf,GEN x,GEN arch,long nba)
       }
       p1 = (GEN)mat[lgmat];
       for (i=1; i<=nba; i++)
-        p1[i] = (gsigne((GEN)alpha[i]) > 0)? zero: un;
+        p1[i] = (gsigne((GEN)alpha[i]) < 0)? 1: 0;
 
-      if (ker_trivial_mod_p(mat, gdeux)) avma = av1;
+      if (u_FpM_deplin(mat, 2)) avma = av1;
       else
       { /* new vector indep. of previous ones */
         avma = av1; alpha = nfun;
@@ -1454,7 +1505,8 @@ zarchstar(GEN nf,GEN x,GEN arch,long nba)
 	if (lgmat > nba)
 	{
           GEN *gptr[2];
-          mat = ginv(mat); gptr[0]=&genarch; gptr[1]=&mat;
+          mat = small_to_mat( u_FpM_inv(mat, 2) );
+          gptr[0]=&genarch; gptr[1]=&mat;
           gerepilemany(av,gptr,2);
 	  y[2]=(long)genarch;
           y[3]=(long)mat; return y;
@@ -1476,7 +1528,7 @@ zinternallog_pk(GEN nf, GEN a0, GEN y, GEN pr, GEN prk
   {
     L = (GEN)list[j]; cyc=(GEN)L[1]; gen=(GEN)L[2];
     if (j==1)
-      p1 = nf_PHlog(nf,a,(GEN)gen[1],pr, nfmodprinit(nf,pr));
+      p1 = nf_PHlog(nf,a,(GEN)gen[1],pr);
     else
     {
       p1 = (GEN)a[1]; a[1] = laddsi(-1,(GEN)a[1]);
@@ -1510,7 +1562,8 @@ static GEN
 zinternallog(GEN nf,GEN a,GEN list_set,long nbgen,GEN arch,GEN fa,long index)
 {
   GEN prlist,ep,y0,y,ainit,list,pr,prk,cyc,psigne,p1,p2;
-  long av,nbp,i,j,k;
+  gpmem_t av;
+  long nbp,i,j,k;
 
   y0 = y = cgetg(nbgen+1,t_COL); av=avma;
   prlist=(GEN)fa[1]; ep=(GEN)fa[2]; nbp=lg(ep)-1;
@@ -1560,7 +1613,8 @@ compute_gen(GEN nf, GEN u1, GEN gen, GEN bid)
 GEN
 zidealstarinitall(GEN nf, GEN ideal,long add_gen)
 {
-  long av = avma,i,j,k,nba,nbp,R1,nbgen,cp;
+  gpmem_t av = avma;
+  long i,j,k,nba,nbp,R1,nbgen,cp;
   GEN p1,y,h,cyc,U,u1,fa,fa2,ep,x,arch,list,sarch,gen;
 
   nf = checknf(nf); R1=itos(gmael(nf,2,1));
@@ -1649,7 +1703,7 @@ zidealstarinit(GEN nf, GEN ideal)
 GEN
 zidealstar(GEN nf, GEN ideal)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN y = zidealstarinitall(nf,ideal,1);
   return gerepilecopy(av,(GEN)y[2]);
 }
@@ -1694,7 +1748,8 @@ check_nfelt(GEN x, GEN *den)
 GEN
 zideallog(GEN nf, GEN x, GEN bid)
 {
-  long av,l,i,N,c;
+  gpmem_t av;
+  long l,i,N,c;
   GEN fa,fa2,ideal,arch,den,p1,cyc,y;
 
   nf=checknf(nf); checkbid(bid);
@@ -1712,6 +1767,7 @@ zideallog(GEN nf, GEN x, GEN bid)
     case t_POLMOD: case t_POL:
       x = algtobasis(nf,x); break;
     case t_MAT:
+      if (lg(x) == 1) return zerocol(c-1);
       return famat_ideallog(nf,(GEN)x[1],(GEN)x[2],bid);
     case t_COL: break;
     default: err(talker,"not an element in zideallog");
@@ -1722,8 +1778,8 @@ zideallog(GEN nf, GEN x, GEN bid)
   {
     GEN g = cgetg(3, t_COL);
     GEN e = cgetg(3, t_COL);
-    g[1] = lmul(x,den); e[1] = un;
-    g[2] = (long)den;   e[2] = lstoi(-1);
+    g[1] = (long)Q_muli_to_int(x,den); e[1] = un;
+    g[2] = (long)den;                  e[2] = lstoi(-1);
     p1 = famat_ideallog(nf, g, e, bid);
   }
   else
@@ -1748,7 +1804,8 @@ zideallog(GEN nf, GEN x, GEN bid)
 GEN
 zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2, long add_gen)
 {
-  long av=avma,i,nbp,nbgen,lx1,lx2,l1,l2,lx;
+  gpmem_t av=avma;
+  long i,nbp,nbgen,lx1,lx2,l1,l2,lx;
   GEN module1,module2,struct1,struct2,fact1,fact2,liste1,liste2;
   GEN module,liste,fact,U,U1,U2,cyc,cyc1,cyc2,uv;
   GEN p1,p2,y,u1,x,fa1,fa2, gen = add_gen? gun: NULL;
@@ -1824,7 +1881,8 @@ zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2, long ad
 GEN
 zidealstarinitjoinarch(GEN nf, GEN bid1, GEN arch, long nba, long add_gen)
 {
-  long av=avma,i,nbp,lx1;
+  gpmem_t av=avma;
+  long i,nbp,lx1;
   GEN module1,struct1,fact1,liste1,U1,U;
   GEN module,liste,cyc,p1,y,u1,x,sarch, gen = add_gen? gun: NULL;
 
@@ -1921,7 +1979,8 @@ static GEN
 ideallistzstarall(GEN bnf,long bound,long flag)
 {
   byteptr ptdif=diffptr;
-  long lim,av0=avma,av,i,j,k,l,q2,lp1,q;
+  gpmem_t lim,av0=avma,av;
+  long i,j,k,l,q2,lp1,q;
   long do_gen = flag & 1, do_units = flag & 2, big_id = !(flag & 4);
   GEN y,nf,p,z,z2,p1,p2,p3,fa,pr,ideal,lu,lu2,funits,racunit,embunit;
 
@@ -1947,7 +2006,7 @@ ideallistzstarall(GEN bnf,long bound,long flag)
   if (bound > (long)maxprime()) err(primer1);
   for (p[2]=0; p[2]<=bound; )
   {
-    p[2] += *ptdif++;
+    NEXT_PRIME_VIADIFF(p[2], ptdif);
     if (DEBUGLEVEL>1) { fprintferr("%ld ",p[2]); flusherr(); }
     fa = primedec(nf,p);
     for (j=1; j<lg(fa); j++)
@@ -2050,6 +2109,7 @@ ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
   long nba,i,j,lx,ly;
   GEN p1,z,p2;
 
+  if (typ(arch) != t_VEC) err(typeer,"ideallistarch");
   nba=0; for (i=1; i<lg(arch); i++) if (signe(arch[i])) nba++;
   lx=lg(list); z=cgetg(lx,t_VEC);
   for (i=1; i<lx; i++)
@@ -2065,7 +2125,7 @@ ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
 static GEN
 ideallistarchall(GEN bnf,GEN list,GEN arch,long flag)
 {
-  ulong av;
+  gpmem_t av;
   long i,j,lp1, do_units = flag & 2;
   GEN nf = checknf(bnf), p1,p2,p3,racunit,funits,lu2,lu,embunit,z,y;