OpenXM_contrib/pari-2.2/src/basemath/base4.c - diff

Return to base4.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/base4.c between version 1.1.1.1 and 1.2

-version 1.1.1.1, 2001/10/02 11:17:02
+version 1.2, 2002/09/11 07:26:49
 Line 22  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 22  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 22  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  #include "pari.h"
  #include "parinf.h"
- #define principalideal_aux(nf,x) (principalideal0((nf),(x),0))
+ extern GEN makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x);
+ extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
+ extern GEN hnf_invimage(GEN A, GEN b);
+ extern int hnfdivide(GEN A, GEN B);
+ extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
+ extern GEN zinternallog_pk(GEN nf,GEN a0,GEN y,GEN pr,GEN prk,GEN list,GEN *psigne);
+ extern GEN special_anti_uniformizer(GEN nf, GEN pr);
+ extern GEN set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch);
+ extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
- extern GEN element_muli(GEN nf, GEN x, GEN y);
+ static GEN mat_ideal_two_elt2(GEN nf, GEN x, GEN a);
- extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);
- static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN idb, GEN *u, GEN *v, GEN *w, GEN *di);
  /*******************************************************************/
  /*                                                                 */
  /*                     IDEAL OPERATIONS                            */
-Line 47  static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN
+Line 52  static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN
 Line 47  static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN
 Line 52  static GEN nfbezout(GEN nf, GEN a, GEN b, GEN ida, GEN
   * (depends on the chosen generator a). All subroutines work with either
   * ideles or ideals (an omitted V is assumed to be 0).
   *
-  * All the output ideals will be in HNF form.
+  * All ideals are output in HNF form. */
-  */
  /* types and conversions */
-Line 85  idealtyp(GEN *ideal, GEN *arch)
+Line 89  idealtyp(GEN *ideal, GEN *arch)
 Line 85  idealtyp(GEN *ideal, GEN *arch)
 Line 89  idealtyp(GEN *ideal, GEN *arch)
    *ideal = x; return t;
  }
- /* Assume ideal in HNF form */
- long
- ideal_is_zk(GEN ideal,long N)
- {
-   long i,j, lx = lg(ideal);
-   if (typ(ideal) != t_MAT || lx==1) return 0;
-   N++; if (lx != N || lg(ideal[1]) != N) return 0;
-   for (i=1; i<N; i++)
-   {
-     if (!gcmp1(gcoeff(ideal,i,i))) return 0;
-     for (j=i+1; j<N; j++)
-       if (!gcmp0(gcoeff(ideal,i,j))) return 0;
-   }
-   return 1;
- }
  static GEN
  prime_to_ideal_aux(GEN nf, GEN vp)
  {
-   GEN m,el;
+   GEN m = eltmul_get_table(nf, (GEN)vp[2]);
-   long i, N = degpol(nf[1]);
+   return hnfmodid(m, (GEN)vp[1]);
-   m = cgetg(N+1,t_MAT); el = (GEN)vp[2];
-   for (i=1; i<=N; i++) m[i] = (long) element_mulid(nf,el,i);
-   return hnfmodid(m,(GEN)vp[1]);
  }
+ /* vp = [a,x,...], a in Z. Return (a,x)  [HACK: vp need not be prime] */
  GEN
  prime_to_ideal(GEN nf, GEN vp)
  {
-   long av=avma;
+   gpmem_t av=avma;
    if (typ(vp) == t_INT) return gscalmat(vp, degpol(nf[1]));
    return gerepileupto(av, prime_to_ideal_aux(nf,vp));
  }
-Line 126  static GEN
+Line 110  static GEN
 Line 126  static GEN
 Line 110  static GEN
  idealmat_to_hnf(GEN nf, GEN x)
  {
    long rx,i,j,N;
-   GEN m,dx;
+   GEN m, cx;
    N=degpol(nf[1]); rx=lg(x)-1;
    if (!rx) return gscalmat(gzero,N);
-   dx=denom(x); if (gcmp1(dx)) dx = NULL; else x=gmul(dx,x);
+   x = Q_primitive_part(x, &cx);
    if (rx >= N) m = x;
    else
    {
-Line 140  idealmat_to_hnf(GEN nf, GEN x)
+Line 124  idealmat_to_hnf(GEN nf, GEN x)
 Line 140  idealmat_to_hnf(GEN nf, GEN x)
 Line 124  idealmat_to_hnf(GEN nf, GEN x)
        for (j=1; j<=N; j++)
          m[(i-1)*N + j] = (long) element_mulid(nf,(GEN)x[i],j);
    }
-   x = hnfmod(m,detint(m));
+   x = hnfmod(m, detint(m));
-   return dx? gdiv(x,dx): x;
+   return cx? gmul(x,cx): x;
  }
  int
-Line 156  ishnfall(GEN x)
+Line 140  ishnfall(GEN x)
 Line 156  ishnfall(GEN x)
 Line 140  ishnfall(GEN x)
    }
    return (gsigne(gcoeff(x,1,1)) > 0);
  }
+ /* same x is known to be integral */
+ int
+ Z_ishnfall(GEN x)
+ {
+   long i,j, lx = lg(x);
+   for (i=2; i<lx; i++)
+   {
+     if (signe(gcoeff(x,i,i)) <= 0) return 0;
+     for (j=1; j<i; j++)
+       if (signe(gcoeff(x,i,j))) return 0;
+   }
+   return (signe(gcoeff(x,1,1)) > 0);
+ }
  GEN
  idealhermite_aux(GEN nf, GEN x)
  {
    long N,tx,lx;
-   GEN z;
+   GEN z, cx;
    tx = idealtyp(&x,&z);
    if (tx == id_PRIME) return prime_to_ideal_aux(nf,x);
    if (tx == id_PRINCIPAL)
    {
-     x = principalideal(nf,x);
+     x = eltmul_get_table(nf, x);
      return idealmat_to_hnf(nf,x);
    }
    N=degpol(nf[1]); lx = lg(x);
-Line 175  idealhermite_aux(GEN nf, GEN x)
+Line 172  idealhermite_aux(GEN nf, GEN x)
 Line 175  idealhermite_aux(GEN nf, GEN x)
 Line 172  idealhermite_aux(GEN nf, GEN x)
    if (lx == N+1 && ishnfall(x)) return x;
    if (lx <= N) return idealmat_to_hnf(nf,x);
-   z=denom(x); if (gcmp1(z)) z=NULL; else x = gmul(z,x);
+   x = Q_primitive_part(x, &cx);
    x = hnfmod(x,detint(x));
-   return z? gdiv(x,z): x;
+   return cx? gmul(x, cx): x;
  }
  GEN
  idealhermite(GEN nf, GEN x)
  {
-   long av=avma;
+   gpmem_t av=avma;
    GEN p1;
    nf = checknf(nf); p1 = idealhermite_aux(nf,x);
    if (p1==x || p1==(GEN)x[1]) return gcopy(p1);
    return gerepileupto(av,p1);
  }
- static GEN
+ GEN
- principalideal0(GEN nf, GEN x, long copy)
+ principalideal(GEN nf, GEN x)
  {
    GEN z = cgetg(2,t_MAT);
+   nf = checknf(nf);
    switch(typ(x))
    {
      case t_INT: case t_FRAC: case t_FRACN:
-       if (copy) x = gcopy(x);
+       x = gscalcol(x, degpol(nf[1])); break;
-       x = gscalcol_i(x, degpol(nf[1])); break;
      case t_POLMOD:
        x = checknfelt_mod(nf,x,"principalideal");
        /* fall through */
      case t_POL:
-       x = copy? algtobasis(nf,x): algtobasis_intern(nf,x);
+       x = algtobasis(nf,x);
        break;
      case t_MAT:
        if (lg(x)!=2) err(typeer,"principalideal");
        x = (GEN)x[1];
      case t_COL:
-       if (lg(x)==lgef(nf[1])-2)
+       if (lg(x)-1==degpol(nf[1])) { x = gcopy(x); break; }
-       {
-         if (copy) x = gcopy(x);
-         break;
-       }
      default: err(typeer,"principalideal");
    }
    z[1]=(long)x; return z;
  }
- GEN
+ static GEN
- principalideal(GEN nf, GEN x)
+ mylog(GEN x, long prec)
  {
-   nf = checknf(nf); return principalideal0(nf,x,1);
+   if (gcmp0(x)) err(precer,"get_arch");
+   return glog(x,prec);
  }
+ GEN get_arch(GEN nf,GEN x,long prec);
  static GEN
- mylog(GEN x, long prec)
+ famat_get_arch(GEN nf, GEN x, long prec)
  {
-   if (gcmp0(x))
+   GEN A, a, g = (GEN)x[1], e = (GEN)x[2];
-     err(precer,"get_arch");
+   long i, l = lg(e);
-   return glog(x,prec);
+   if (l <= 1)
+     return get_arch(nf, gun, prec);
+   A = NULL; /* -Wall */
+   for (i=1; i<l; i++)
+   {
+     a = get_arch(nf, (GEN)g[i], prec);
+     a = gmul((GEN)e[i], a);
+     A = (i == 1)? a: gadd(A, a);
+   }
+   return A;
  }
- /* for internal use */
+ static GEN
+ scalar_get_arch(long R1, long RU, GEN x, long prec)
+ {
+   GEN v = cgetg(RU+1,t_VEC);
+   GEN p1 = glog(x, prec);
+   long i;
+   for (i=1; i<=R1; i++) v[i] = (long)p1;
+   if (i <= RU)
+   {
+     p1 = gmul2n(p1,1);
+     for (   ; i<=RU; i++) v[i] = (long)p1;
+   }
+   return v;
+ }
+ /* For internal use. Get archimedean components: [e_i log( sigma_i(x) )],
+  * with e_i = 1 (resp 2.) for i <= R1 (resp. > R1) */
  GEN
  get_arch(GEN nf,GEN x,long prec)
  {
-   long i,R1,RU;
+   long i, RU, R1 = nf_get_r1(nf);
-   GEN v,p1,p2;
+   GEN v;
-   R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
+   RU = lg(nf[6]) - 1;
-   if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
+   switch(typ(x))
+   {
+     case t_MAT: return famat_get_arch(nf,x,prec);
+     case t_POLMOD:
+     case t_POL: x = algtobasis_i(nf,x);   /* fall through */
+     case t_COL: if (!isnfscalar(x)) break;
+       x = (GEN)x[1]; /* fall through */
+     default: /* rational number */
+       return scalar_get_arch(R1, RU, x, prec);
+   }
    v = cgetg(RU+1,t_VEC);
    if (isnfscalar(x)) /* rational number */
    {
-     p1 = glog((GEN)x[1],prec);
+     GEN p1 = glog((GEN)x[1], prec);
-     p2 = (RU > R1)? gmul2n(p1,1): NULL;
-     for (i=1; i<=R1; i++) v[i]=(long)p1;
+     for (i=1; i<=R1; i++) v[i] = (long)p1;
-     for (   ; i<=RU; i++) v[i]=(long)p2;
+     if (i <= RU)
+     {
+       p1 = gmul2n(p1,1);
+       for (   ; i<=RU; i++) v[i] = (long)p1;
+     }
    }
-   else
+   x = gmul(gmael(nf,5,1),x);
+   for (i=1; i<=R1; i++) v[i] = (long)mylog((GEN)x[i],prec);
+   for (   ; i<=RU; i++) v[i] = lmul2n(mylog((GEN)x[i],prec),1);
+   return v;
+ }
+ GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
+ static GEN
+ famat_get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
+ {
+   GEN A, T, a, t, g = (GEN)x[1], e = (GEN)x[2];
+   long i, l = lg(e);
+   if (l <= 1)
+     return get_arch_real(nf, gun, emb, prec);
+   A = T = NULL; /* -Wall */
+   for (i=1; i<l; i++)
    {
-     x = gmul(gmael(nf,5,1),x);
+     a = get_arch_real(nf, (GEN)g[i], &t, prec);
-     for (i=1; i<=R1; i++) v[i] = (long)mylog((GEN)x[i],prec);
+     if (!a) return NULL;
-     for (   ; i<=RU; i++) v[i] = lmul2n(mylog((GEN)x[i],prec),1);
+     a = gmul((GEN)e[i], a);
+     t = vecpow(t, (GEN)e[i]);
+     if (i == 1) { A = a;          T = t; }
+     else        { A = gadd(A, a); T = vecmul(T, t); }
    }
-   return v;
+   *emb = T; return A;
  }
+ static GEN
+ scalar_get_arch_real(long R1, long RU, GEN u, GEN *emb, long prec)
+ {
+   GEN v, x, p1;
+   long i, s;
+   s = gsigne(u);
+   if (!s) err(talker,"0 in get_arch_real");
+   x = cgetg(RU+1, t_COL);
+   for (i=1; i<=RU; i++) x[i] = (long)u;
+   v = cgetg(RU+1, t_COL);
+   p1 = (s > 0)? glog(u,prec): gzero;
+   for (i=1; i<=R1; i++) v[i] = (long)p1;
+   if (i <= RU)
+   {
+     p1 = gmul2n(p1,1);
+     for (   ; i<=RU; i++) v[i] = (long)p1;
+   }
+   *emb = x; return v;
+ }
+ static int
+ low_prec(GEN x)
+ {
+   return gcmp0(x) || (typ(x) == t_REAL && lg(x) == 3);
+ }
  /* as above but return NULL if precision problem, and set *emb to the
   * embeddings of x */
  GEN
- get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
+ get_arch_real(GEN nf, GEN x, GEN *emb, long prec)
  {
-   long i,R1,RU;
+   long i, RU, R1 = nf_get_r1(nf);
-   GEN v,p1,p2;
+   GEN v, t;
-   R1=itos(gmael(nf,2,1)); RU = R1+itos(gmael(nf,2,2));
+   RU = lg(nf[6])-1;
-   if (typ(x)!=t_COL) x = algtobasis_intern(nf,x);
+   switch(typ(x))
+   {
+     case t_MAT: return famat_get_arch_real(nf,x,emb,prec);
+     case t_POLMOD:
+     case t_POL: x = algtobasis_i(nf,x);   /* fall through */
+     case t_COL: if (!isnfscalar(x)) break;
+       x = (GEN)x[1]; /* fall through */
+     default: /* rational number */
+       return scalar_get_arch_real(R1, RU, x, emb, prec);
+   }
    v = cgetg(RU+1,t_COL);
-   if (isnfscalar(x)) /* rational number */
+   x = gmul(gmael(nf,5,1), x);
+   for (i=1; i<=R1; i++)
    {
-     GEN u = (GEN)x[1];
+     t = gabs((GEN)x[i],prec); if (low_prec(t)) return NULL;
-     i = signe(u);
+     v[i] = llog(t,prec);
-     if (!i) err(talker,"0 in get_arch_real");
-     p1= (i > 0)? glog(u,prec): gzero;
-     p2 = (RU > R1)? gmul2n(p1,1): NULL;
-     for (i=1; i<=R1; i++) v[i]=(long)p1;
-     for (   ; i<=RU; i++) v[i]=(long)p2;
    }
-   else
+   for (   ; i<=RU; i++)
    {
-     GEN t;
+     t = gnorm((GEN)x[i]); if (low_prec(t)) return NULL;
-     x = gmul(gmael(nf,5,1),x);
+     v[i] = llog(t,prec);
-     for (i=1; i<=R1; i++)
-     {
-       t = gabs((GEN)x[i],prec); if (gcmp0(t)) return NULL;
-       v[i] = llog(t,prec);
-     }
-     for (   ; i<=RU; i++)
-     {
-       t = gnorm((GEN)x[i]); if (gcmp0(t)) return NULL;
-       v[i] = llog(t,prec);
-     }
    }
    *emb = x; return v;
  }
-Line 303  get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
+Line 385  get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
 Line 303  get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
 Line 385  get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
  GEN
  principalidele(GEN nf, GEN x, long prec)
  {
-   GEN p1,y = cgetg(3,t_VEC);
+   GEN p1, y = cgetg(3,t_VEC);
-   long av;
+   gpmem_t av;
-   nf = checknf(nf);
+   p1 = principalideal(nf,x);
-   p1 = principalideal0(nf,x,1);
    y[1] = (long)p1;
-   av =avma; p1 = get_arch(nf,(GEN)p1[1],prec);
+   av = avma; p1 = get_arch(nf,(GEN)p1[1],prec);
    y[2] = lpileupto(av,p1); return y;
  }
-Line 355  idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
+Line 436  idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
 Line 355  idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
 Line 436  idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
    return idealaddtoone(nf,arg1,arg2);
  }
- static GEN
- two_to_hnf(GEN nf, GEN a, GEN b)
- {
-   a = principalideal_aux(nf,a);
-   b = principalideal_aux(nf,b);
-   a = concatsp(a,b);
-   if (lgef(nf[1])==5) /* quadratic field: a has to be turned into idealmat */
-     a = idealmul(nf,idmat(2),a);
-   return idealmat_to_hnf(nf, a);
- }
  GEN
  idealhnf0(GEN nf, GEN a, GEN b)
  {
-   long av;
+   gpmem_t av;
+   GEN x;
    if (!b) return idealhermite(nf,a);
    /* HNF of aZ_K+bZ_K */
-   av = avma; nf=checknf(nf);
+   av = avma; nf = checknf(nf);
-   return gerepileupto(av, two_to_hnf(nf,a,b));
+   x = concatsp(eltmul_get_table(nf,a), eltmul_get_table(nf,b));
+   return gerepileupto(av, idealmat_to_hnf(nf, x));
  }
  GEN
-Line 386  idealhermite2(GEN nf, GEN a, GEN b)
+Line 458  idealhermite2(GEN nf, GEN a, GEN b)
 Line 386  idealhermite2(GEN nf, GEN a, GEN b)
 Line 458  idealhermite2(GEN nf, GEN a, GEN b)
  static int
  ok_elt(GEN x, GEN xZ, GEN y)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
    int r = gegal(x, hnfmodid(y, xZ));
    avma = av; return r;
  }
- static void
- setprec(GEN x, long prec)
- {
-   long i,j, n=lg(x);
-   for (i=1;i<n;i++)
-   {
-     GEN p2,p1 = (GEN)x[i];
-     for (j=1;j<n;j++)
-     {
-       p2 = (GEN)p1[j];
-       if (typ(p2) == t_REAL) setlg(p2, prec);
-     }
-   }
- }
- /* find a basis of x whose elements have small norm
-  * M a bound for the size of coeffs of x */
- GEN
- ideal_better_basis(GEN nf, GEN x, GEN M)
- {
-   GEN a,b;
-   long nfprec = nfgetprec(nf);
-   long prec = DEFAULTPREC + (expi(M) >> TWOPOTBITS_IN_LONG);
-   if (typ(nf[5]) != t_VEC) return x;
-   if ((prec<<1) < nfprec) prec = (prec+nfprec) >> 1;
-   a = qf_base_change(gmael(nf,5,3),x,1);
-   setprec(a,prec);
-   b = lllgramintern(a,4,1,prec);
-   if (!b)
-   {
-     if (DEBUGLEVEL)
-       err(warner, "precision too low in ideal_better_basis (1)");
-     if (nfprec > prec)
-     {
-       setprec(a,nfprec);
-       b = lllgramintern(a,4,1,nfprec);
-     }
-   }
-   if (!b)
-   {
-     if (DEBUGLEVEL)
-       err(warner, "precision too low in ideal_better_basis (2)");
-     b = lllint(x); /* better than nothing */
-   }
-   return gmul(x, b);
- }
  static GEN
  addmul_col(GEN a, long s, GEN b)
  {
-Line 468  addmul_mat(GEN a, long s, GEN b)
+Line 492  addmul_mat(GEN a, long s, GEN b)
 Line 468  addmul_mat(GEN a, long s, GEN b)
 Line 492  addmul_mat(GEN a, long s, GEN b)
  static GEN
  mat_ideal_two_elt(GEN nf, GEN x)
  {
-   GEN y,a,beta,cx,xZ,mul, pol = (GEN)nf[1];
+   GEN y,a,beta,cx,xZ,mul;
-   long i,j,lm, N = degpol(pol);
+   long i,lm, N = degpol(nf[1]);
-   ulong av,tetpil;
+   gpmem_t av0,av,tetpil;
-   y=cgetg(3,t_VEC); av=avma;
+   y = cgetg(3,t_VEC); av = avma;
-   if (lg(x[1])!=N+1) err(typeer,"ideal_two_elt");
+   if (lg(x[1]) != N+1) err(typeer,"ideal_two_elt");
    if (N == 2)
    {
      y[1] = lcopy(gcoeff(x,1,1));
      y[2] = lcopy((GEN)x[2]); return y;
    }
-   cx = content(x); if (!gcmp1(cx)) x = gdiv(x,cx);
+   x = Q_primitive_part(x, &cx); if (!cx) cx = gun;
    if (lg(x) != N+1) x = idealhermite_aux(nf,x);
    xZ = gcoeff(x,1,1);
    if (gcmp1(xZ))
    {
      y[1] = lpilecopy(av,cx);
-     y[2] = (long)gscalcol(cx,N); return y;
+     y[2] = (long)gscalcol(cx, N); return y;
    }
+   av0 = avma;
    a = NULL; /* gcc -Wall */
    beta= cgetg(N+1, t_VEC);
    mul = cgetg(N+1, t_VEC); lm = 1; /* = lg(mul) */
    /* look for a in x such that a O/xZ = x O/xZ */
    for (i=2; i<=N; i++)
    {
-     GEN t, y = cgetg(N+1,t_MAT);
+     gpmem_t av1 = avma;
-     a = (GEN)x[i];
+     GEN t, y = eltmul_get_table(nf, (GEN)x[i]);
-     for (j=1; j<=N; j++) y[j] = (long)element_mulid(nf,a,j);
+     t = FpM_red(y, xZ);
-     /* columns of mul[i] = canonical generators for x[i] O/xZ as Z-module */
+     if (gcmp0(t)) { avma = av1; continue; }
-     t = gmod(y, xZ);
+     if (ok_elt(x,xZ, t)) { a = (GEN)x[i]; break; }
-     if (gcmp0(t)) continue;
-     if (ok_elt(x,xZ, t)) break;
      beta[lm]= x[i];
+     /* mul[i] = { canonical generators for x[i] O/xZ as Z-module } */
      mul[lm] = (long)t; lm++;
    }
-   if (i>N)
+   if (i > N)
    {
      GEN z = cgetg(lm, t_VECSMALL);
-     ulong av1, c = 0;
+     gpmem_t av1;
+     ulong c = 0;
      setlg(mul, lm);
      setlg(beta,lm);
-     if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case: ");
+     if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case:\n");
      for(av1=avma;;avma=av1)
      {
-       if (DEBUGLEVEL>3) fprintferr("%ld ", ++c);
+       if (++c == 100) {
+         if (DEBUGLEVEL>3) fprintferr("using approximation theorem\n");
+         a = mat_ideal_two_elt2(nf, x, xZ); goto END;
+       }
        for (a=NULL,i=1; i<lm; i++)
        {
          long t = (mymyrand() >> (BITS_IN_RANDOM-5)) - 7; /* in [-7,8] */
-Line 526  mat_ideal_two_elt(GEN nf, GEN x)
+Line 555  mat_ideal_two_elt(GEN nf, GEN x)
 Line 526  mat_ideal_two_elt(GEN nf, GEN x)
 Line 555  mat_ideal_two_elt(GEN nf, GEN x)
      }
      for (a=NULL,i=1; i<lm; i++)
        a = addmul_col(a, z[i], (GEN)beta[i]);
-     if (DEBUGLEVEL>3) fprintferr("\n");
    }
-   a = centermod(a, xZ); tetpil=avma;
+ END:
+   a = centermod(a, xZ);
+   tetpil = avma;
    y[1] = lmul(xZ,cx);
    y[2] = lmul(a, cx);
    gerepilemanyvec(av,tetpil,y+1,2); return y;
  }
- /* Etant donne un ideal ix, ressort un vecteur [a,alpha] a deux composantes
+ /* Given an ideal x, returns [a,alpha] such that a is in Q,
-  * tel que a soit rationnel et ix=aZ_K+alpha Z_K, alpha etant un vecteur
+  * x = a Z_K + alpha Z_K, alpha in K^*
-  * colonne sur la base d'entiers. On peut avoir a=0 ou alpha=0, mais on ne
+  * a = 0 or alpha = 0 are possible, but do not try to determine whether
-  * cherche pas a determiner si ix est principal.
+  * x is principal. */
-  */
  GEN
  ideal_two_elt(GEN nf, GEN x)
  {
-Line 604  factor_norm(GEN x)
+Line 633  factor_norm(GEN x)
 Line 604  factor_norm(GEN x)
 Line 633  factor_norm(GEN x)
  GEN
  idealfactor(GEN nf, GEN x)
  {
-   long av,tx, tetpil,i,j,k,lf,lc,N,l,v,vc,e;
+   gpmem_t av,tetpil;
+   long tx,i,j,k,lf,lc,N,l,v,vc,e;
    GEN f,f1,f2,c1,c2,y1,y2,y,p1,p2,cx,P;
    tx = idealtyp(&x,&y);
    if (tx == id_PRIME)
    {
      y=cgetg(3,t_MAT);
-     y[1]=lgetg(2,t_COL); mael(y,1,1)=lcopy(x);
+     y[1]=(long)_col(gcopy(x));
-     y[2]=lgetg(2,t_COL); mael(y,2,1)=un; return y;
+     y[2]=(long)_col(gun); return y;
    }
-   nf=checknf(nf); av=avma;
+   av = avma;
-   if (tx == id_PRINCIPAL) x = principalideal_aux(nf,x);
+   nf = checknf(nf);
+   N = degpol(nf[1]);
-   N=degpol(nf[1]); if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
+   if (tx == id_PRINCIPAL) x = idealhermite_aux(nf, x);
+   else
+     if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
    if (lg(x)==1) err(talker,"zero ideal in idealfactor");
-   cx = content(x);
+   x = Q_primitive_part(x, &cx);
-   if (gcmp1(cx))
+   if (!cx)
    {
      c1 = c2 = NULL; /* gcc -Wall */
-     lc=1;
+     lc = 1;
    }
    else
    {
-     f = factor(cx); x = gdiv(x,cx);
+     f = factor(cx);
      c1 = (GEN)f[1];
      c2 = (GEN)f[2]; lc = lg(c1);
    }
-Line 641  idealfactor(GEN nf, GEN x)
+Line 673  idealfactor(GEN nf, GEN x)
 Line 641  idealfactor(GEN nf, GEN x)
 Line 673  idealfactor(GEN nf, GEN x)
    {
      p1 = primedec(nf,(GEN)f1[i]);
      l = f2[i]; /* = v_p(Nx) */
-     vc = ggval(cx,(GEN)f1[i]);
+     vc = cx? ggval(cx,(GEN)f1[i]): 0;
      for (j=1; j<lg(p1); j++)
      {
        P = (GEN)p1[j]; e = itos((GEN)P[3]);
-Line 685  idealfactor(GEN nf, GEN x)
+Line 717  idealfactor(GEN nf, GEN x)
 Line 685  idealfactor(GEN nf, GEN x)
 Line 717  idealfactor(GEN nf, GEN x)
  long
  idealval(GEN nf, GEN ix, GEN P)
  {
-   long N,v,vd,w,av=avma,av1,lim,e,f,i,j,k, tx = typ(ix);
+   gpmem_t av=avma,av1,lim;
+   long N,v,vd,w,e,f,i,j,k, tx = typ(ix);
    GEN mul,mat,a,x,y,r,bp,p,pk,cx;
    nf=checknf(nf); checkprimeid(P);
    if (is_extscalar_t(tx) || tx==t_COL) return element_val(nf,ix,P);
    p=(GEN)P[1]; N=degpol(nf[1]);
    tx = idealtyp(&ix,&a);
-   cx = content(ix); if (!gcmp1(cx)) ix = gdiv(ix,cx);
+   ix = Q_primitive_part(ix, &cx);
    if (tx != id_MAT)
      ix = idealhermite_aux(nf,ix);
    else
-Line 707  idealval(GEN nf, GEN ix, GEN P)
+Line 740  idealval(GEN nf, GEN ix, GEN P)
 Line 707  idealval(GEN nf, GEN ix, GEN P)
 Line 740  idealval(GEN nf, GEN ix, GEN P)
    /* 0 <= ceil[v_P(ix) / e] <= v_p(ix \cap Z) --> v_P <= e * v_p */
    j = k * e;
    v = min(i,j); /* v_P(ix) <= v */
-   vd = ggval(cx,p) * e;
+   vd = cx? ggval(cx,p) * e: 0;
    if (!v) { avma = av; return vd; }
    mul = cgetg(N+1,t_MAT); bp=(GEN)P[5];
-Line 762  idealval(GEN nf, GEN ix, GEN P)
+Line 795  idealval(GEN nf, GEN ix, GEN P)
 Line 762  idealval(GEN nf, GEN ix, GEN P)
 Line 795  idealval(GEN nf, GEN ix, GEN P)
  GEN
  idealadd(GEN nf, GEN x, GEN y)
  {
-   long av=avma,N,tx,ty;
+   gpmem_t av=avma;
+   long N,tx,ty;
    GEN z,p1,dx,dy,dz;
    int modid;
-Line 774  idealadd(GEN nf, GEN x, GEN y)
+Line 808  idealadd(GEN nf, GEN x, GEN y)
 Line 774  idealadd(GEN nf, GEN x, GEN y)
 Line 808  idealadd(GEN nf, GEN x, GEN y)
    if (ty != id_MAT || lg(y)!=N+1) y = idealhermite_aux(nf,y);
    if (lg(x) == 1) return gerepileupto(av,y);
    if (lg(y) == 1) return gerepileupto(av,x); /* check for 0 ideal */
-   dx=denom(x);
+   dx = denom(x);
-   dy=denom(y); dz=mulii(dx,dy);
+   dy = denom(y); dz = mulii(dx,dy);
-   if (gcmp1(dz)) dz = NULL; else { x=gmul(x,dz); y=gmul(y,dz); }
+   if (gcmp1(dz)) dz = NULL; else {
+     x = Q_muli_to_int(x, dz);
+     y = Q_muli_to_int(y, dz);
+   }
    if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
    {
-     p1 = mppgcd(gcoeff(x,1,1),gcoeff(y,1,1));
+     p1 = mppgcd(gcoeff(x,1,1), gcoeff(y,1,1));
      modid = 1;
    }
    else
    {
-     p1 = mppgcd(detint(x),detint(y));
+     p1 = mppgcd(detint(x), detint(y));
      modid = 0;
    }
    if (gcmp1(p1))
    {
      long i,j;
      if (!dz) { avma=av; return idmat(N); }
-     avma = (long)dz; dz = gerepileupto((long)z, ginv(dz));
+     avma = (gpmem_t)dz; dz = gerepileupto((gpmem_t)z, ginv(dz));
      for (i=1; i<=N; i++)
      {
        z[i]=lgetg(N+1,t_COL);
-Line 806  idealadd(GEN nf, GEN x, GEN y)
+Line 843  idealadd(GEN nf, GEN x, GEN y)
 Line 806  idealadd(GEN nf, GEN x, GEN y)
 Line 843  idealadd(GEN nf, GEN x, GEN y)
    return gerepileupto(av,z);
  }
+ /* Assume x,y integral non zero (presumably coprime, which is checked).
+  * Return a in x such that 1-a in y. If (x + y) \cap Z is 1, then addone_aux
+  * is more efficient. */
+ GEN
+ addone_aux2(GEN x, GEN y)
+ {
+   GEN U, H;
+   long i, l;
+   H = hnfall_i(concatsp(x,y), &U, 1);
+   l = lg(H);
+   for (i=1; i<l; i++)
+     if (!gcmp1(gcoeff(H,i,i)))
+       err(talker,"ideals don't sum to Z_K in idealaddtoone");
+   U = (GEN)U[l]; setlg(U, lg(x));
+   return gmul(x, U);
+ }
+ /* assume x,y integral, non zero in HNF */
  static GEN
- get_p1(GEN nf, GEN x, GEN y,long fl)
+ addone_aux(GEN x, GEN y)
  {
-   GEN u,v,v1,v2,v3,v4;
+   GEN a, b, u, v;
-   long i,j,N;
-   switch(fl)
+   a = gcoeff(x,1,1);
-   {
+   b = gcoeff(y,1,1);
-     case 1:
+   if (typ(a) != t_INT || typ(b) != t_INT)
-       v1 = gcoeff(x,1,1);
+     err(talker,"ideals don't sum to Z_K in idealaddtoone");
-       v2 = gcoeff(y,1,1);
+   if (gcmp1(bezout(a,b,&u,&v))) return gmul(u,(GEN)x[1]);
-       if (typ(v1)!=t_INT || typ(v2)!=t_INT)
+   return addone_aux2(x, y);
-         err(talker,"ideals don't sum to Z_K in idealaddtoone");
-       if (gcmp1(bezout(v1,v2,&u,&v)))
-         return gmul(u,(GEN)x[1]);
-     default:
-       v=hnfperm(concatsp(x,y));
-       v1=(GEN)v[1]; v2=(GEN)v[2]; v3=(GEN)v[3];
-       j=0; N = degpol(nf[1]);
-       for (i=1; i<=N; i++)
-       {
-         if (!gcmp1(gcoeff(v1,i,i)))
-           err(talker,"ideals don't sum to Z_K in idealaddtoone");
-         if (gcmp1((GEN)v3[i])) j=i;
-       }
-       v4=(GEN)v2[N+j]; setlg(v4,N+1);
-       return gmul(x,v4);
-   }
  }
  GEN
-Line 849  idealaddtoone_i(GEN nf, GEN x, GEN y)
+Line 888  idealaddtoone_i(GEN nf, GEN x, GEN y)
 Line 849  idealaddtoone_i(GEN nf, GEN x, GEN y)
 Line 888  idealaddtoone_i(GEN nf, GEN x, GEN y)
      fprintferr(" y = %Z\n",y);
    }
    t = idealtyp(&x,&p1);
-   if (t != id_MAT || lg(x) > 1 || lg(x) != lg(x[1]))
+   if (t != id_MAT || lg(x) == 1 || lg(x) != lg(x[1]))
      xh = idealhermite_aux(nf,x);
    else
-     { xh=x; fl = isnfscalar((GEN)x[1]); }
+     { xh = x; fl = isnfscalar((GEN)x[1]); }
    t = idealtyp(&y,&p1);
    if (t != id_MAT || lg(y) == 1 || lg(y) != lg(y[1]))
      yh = idealhermite_aux(nf,y);
    else
-     { yh=y; if (fl) fl = isnfscalar((GEN)y[1]); }
+     { yh = y; if (fl) fl = isnfscalar((GEN)y[1]); }
    if (lg(xh) == 1)
    {
      if (lg(yh) == 1 || !gcmp1(gabs(gcoeff(yh,1,1),0)))
-Line 872  idealaddtoone_i(GEN nf, GEN x, GEN y)
+Line 911  idealaddtoone_i(GEN nf, GEN x, GEN y)
 Line 872  idealaddtoone_i(GEN nf, GEN x, GEN y)
 Line 911  idealaddtoone_i(GEN nf, GEN x, GEN y)
      return algtobasis(nf, gneg(p1));
    }
-   p1 = get_p1(nf,xh,yh,fl);
+   if (fl) p1 = addone_aux (xh,yh); /* xh[1], yh[1] scalar */
+   else    p1 = addone_aux2(xh,yh);
    p1 = element_reduce(nf,p1, idealmullll(nf,x,y));
    if (DEBUGLEVEL>4 && !gcmp0(p1))
      fprintferr(" leaving idealaddtoone: %Z\n",p1);
-Line 909  ideleaddone_aux(GEN nf,GEN x,GEN ideal)
+Line 949  ideleaddone_aux(GEN nf,GEN x,GEN ideal)
 Line 909  ideleaddone_aux(GEN nf,GEN x,GEN ideal)
 Line 949  ideleaddone_aux(GEN nf,GEN x,GEN ideal)
    return nba? p3: gcopy(p3);
  }
- static GEN
+ GEN
  unnf_minus_x(GEN x)
  {
    long i, N = lg(x);
-Line 923  unnf_minus_x(GEN x)
+Line 963  unnf_minus_x(GEN x)
 Line 923  unnf_minus_x(GEN x)
 Line 963  unnf_minus_x(GEN x)
  static GEN
  addone(GEN f(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
  {
-   GEN z = cgetg(3,t_VEC);
+   GEN z = cgetg(3,t_VEC), a;
-   long av=avma;
+   gpmem_t av = avma;
+   nf = checknf(nf);
+   a = gerepileupto(av, f(nf,x,y));
+   z[1]=(long)a;
+   z[2]=(long)unnf_minus_x(a); return z;
+ }
-   nf=checknf(nf); x = gerepileupto(av, f(nf,x,y));
+ /* assume x,y HNF, non-zero */
-   z[1]=(long)x; z[2]=(long)unnf_minus_x(x); return z;
+ static GEN
+ addone_nored(GEN x, GEN y)
+ {
+   GEN z = cgetg(3,t_VEC), a;
+   gpmem_t av = avma;
+   a = gerepileupto(av, addone_aux(x,y));
+   z[1] = (long)a;
+   z[2] = (long)unnf_minus_x(a); return z;
  }
  GEN
-Line 942  ideleaddone(GEN nf,GEN x,GEN idele)
+Line 994  ideleaddone(GEN nf,GEN x,GEN idele)
 Line 942  ideleaddone(GEN nf,GEN x,GEN idele)
 Line 994  ideleaddone(GEN nf,GEN x,GEN idele)
    return addone(ideleaddone_aux,nf,x,idele);
  }
- /* return integral x = 0 mod p/pr^e, (x,pr) = 1.
-  * Don't reduce mod p here: caller may need result mod pr^k */
- GEN
- special_anti_uniformizer(GEN nf, GEN pr)
- {
-   GEN p = (GEN)pr[1], e = (GEN)pr[3];
-   return gdivexact(element_pow(nf,(GEN)pr[5],e), gpuigs(p,itos(e)-1));
- }
- GEN
- nfmodprinit(GEN nf, GEN pr)
- {
-   long av;
-   GEN p,p1,prhall;
-   nf = checknf(nf); checkprimeid(pr);
-   prhall = cgetg(3,t_VEC);
-   prhall[1] = (long) prime_to_ideal(nf,pr);
-   av = avma; p = (GEN)pr[1];
-   p1 = cgetg(2,t_MAT);
-   p1[1] = (long)gmod(special_anti_uniformizer(nf, pr), p);
-   p1 = hnfmodid(idealhermite_aux(nf,p1), p);
-   p1 = idealaddtoone_i(nf,pr,p1);
-   /* p1 = 1 mod pr, p1 = 0 mod q^{e_q} for all other primes q | p */
-   prhall[2] = lpileupto(av, unnf_minus_x(p1)); return prhall;
- }
  /* given an element x in Z_K and an integral ideal y with x, y coprime,
     outputs an element inverse of x modulo y */
  GEN
  element_invmodideal(GEN nf, GEN x, GEN y)
  {
-   long av=avma,N,i, fl = 1;
+   gpmem_t av=avma;
+   long N,i, fl = 1;
    GEN v,p1,xh,yh;
    nf=checknf(nf); N=degpol(nf[1]);
-   if (ideal_is_zk(y,N)) return zerocol(N);
+   if (gcmp1(gcoeff(y,1,1))) return zerocol(N);
-   if (DEBUGLEVEL>4)
-   {
-     fprintferr(" entree dans element_invmodideal() :\n");
-     fprintferr(" x = "); outerr(x);
-     fprintferr(" y = "); outerr(y);
-   }
    i = lg(y);
    if (typ(y)!=t_MAT || i==1 || i != lg(y[1])) yh=idealhermite_aux(nf,y);
    else
-Line 998  element_invmodideal(GEN nf, GEN x, GEN y)
+Line 1016  element_invmodideal(GEN nf, GEN x, GEN y)
 Line 998  element_invmodideal(GEN nf, GEN x, GEN y)
 Line 1016  element_invmodideal(GEN nf, GEN x, GEN y)
      default: err(typeer,"element_invmodideal");
        return NULL; /* not reached */
    }
-   p1 = get_p1(nf,xh,yh,fl);
+   if (fl) p1 = addone_aux (xh,yh);
+   else    p1 = addone_aux2(xh,yh);
    p1 = element_div(nf,p1,x);
    v = gerepileupto(av, nfreducemodideal(nf,p1,y));
-   if (DEBUGLEVEL>2)
-     { fprintferr(" sortie de element_invmodideal : v = "); outerr(v); }
    return v;
  }
  GEN
  idealaddmultoone(GEN nf, GEN list)
  {
-   long av=avma,tetpil,N,i,i1,j,k;
+   gpmem_t av=avma,tetpil;
+   long N,i,i1,j,k;
    GEN z,v,v1,v2,v3,p1;
    nf=checknf(nf); N=degpol(nf[1]);
-Line 1053  idealaddmultoone(GEN nf, GEN list)
+Line 1071  idealaddmultoone(GEN nf, GEN list)
 Line 1053  idealaddmultoone(GEN nf, GEN list)
 Line 1071  idealaddmultoone(GEN nf, GEN list)
  /* multiplication */
  /* x integral ideal (without archimedean component) in HNF form
-  * [a,alpha,n] corresponds to the ideal aZ_K+alpha Z_K (a is a
+  * y = [a,alpha] corresponds to the ideal aZ_K+alpha Z_K (a is a
   * rational integer). Multiply them
   */
  static GEN
- idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
+ idealmulspec(GEN nf, GEN x, GEN y)
  {
    long i, N=lg(x)-1;
-   GEN m, mod;
+   GEN m, mod, a = (GEN)y[1], alpha = (GEN)y[2];
    if (isnfscalar(alpha))
-     return gmul(mppgcd(a,(GEN)alpha[1]),x);
+     return gmul(mppgcd(a, (GEN)alpha[1]),x);
    mod = mulii(a, gcoeff(x,1,1));
    m = cgetg((N<<1)+1,t_MAT);
    for (i=1; i<=N; i++) m[i]=(long)element_muli(nf,alpha,(GEN)x[i]);
-Line 1074  idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
+Line 1092  idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
 Line 1074  idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
 Line 1092  idealmulspec(GEN nf, GEN x, GEN a, GEN alpha)
  /* x ideal (matrix form,maximal rank), vp prime ideal (primedec). Output the
   * product. Can be used for arbitrary vp of the form [p,a,e,f,b], IF vp
   * =pZ_K+aZ_K, p is an integer, and norm(vp) = p^f; e and b are not used.
-  * For internal use.
+  * For internal use. */
-  */
  GEN
  idealmulprime(GEN nf, GEN x, GEN vp)
  {
-   GEN denx = denom(x);
+   GEN cx;
+   x = Q_primitive_part(x, &cx);
+   x = idealmulspec(nf,x,vp);
+   return cx? gmul(x,cx): x;
+ }
-   if (gcmp1(denx)) denx = NULL; else x = gmul(denx,x);
+ GEN arch_mul(GEN x, GEN y);
-   x = idealmulspec(nf,x, (GEN)vp[1], (GEN)vp[2]);
+ GEN vecpow(GEN x, GEN n);
-   return denx? gdiv(x,denx): x;
+ GEN vecmul(GEN x, GEN y);
+ static GEN
+ mul(GEN nf, GEN x, GEN y)
+ {
+   GEN yZ = gcoeff(y,1,1);
+   if (is_pm1(yZ)) return gcopy(x);
+   y = mat_ideal_two_elt(nf,y);
+   return idealmulspec(nf, x, y);
  }
  /* Assume ix and iy are integral in HNF form (or ideles of the same form).
-Line 1099  idealmulh(GEN nf, GEN ix, GEN iy)
+Line 1128  idealmulh(GEN nf, GEN ix, GEN iy)
 Line 1099  idealmulh(GEN nf, GEN ix, GEN iy)
 Line 1128  idealmulh(GEN nf, GEN ix, GEN iy)
    if (typ(iy)==t_VEC && typ(iy[1]) != t_INT) {f+=2; y=(GEN)iy[1];} else y=iy;
    res = f? cgetg(3,t_VEC): NULL;
-   if (typ(y) != t_VEC) y = ideal_two_elt(nf,y);
+   if (typ(y) == t_VEC)
-   y = idealmulspec(nf,x,(GEN)y[1],(GEN)y[2]);
+     y = idealmulspec(nf,x,y);
+   else
+   {
+     GEN xZ = gcoeff(x,1,1);
+     GEN yZ = gcoeff(x,1,1);
+     y = (cmpii(xZ, yZ) < 0)? mul(nf,y,x): mul(nf,x,y);
+   }
    if (!f) return y;
-   res[1]=(long)y;
+   res[1] = (long)y;
-   if (f==3) y = gadd((GEN)ix[2],(GEN)iy[2]);
+   if (f==3) y = arch_mul((GEN)ix[2], (GEN)iy[2]);
    else
    {
      y = (f==2)? (GEN)iy[2]: (GEN)ix[2];
      y = gcopy(y);
    }
-   res[2]=(long)y; return res;
+   res[2] = (long)y; return res;
  }
+ GEN
+ mul_content(GEN cx, GEN cy)
+ {
+   if (!cx) return cy;
+   if (!cy) return cx;
+   return gmul(cx,cy);
+ }
  /* x and y are ideals in matrix form */
  static GEN
  idealmat_mul(GEN nf, GEN x, GEN y)
  {
-   long i,j, rx=lg(x)-1, ry=lg(y)-1;
+   long i,j, rx = lg(x)-1, ry = lg(y)-1;
-   GEN dx,dy,m;
+   GEN cx, cy, m;
-   dx=denom(x); if (!gcmp1(dx)) x=gmul(dx,x);
+   x = Q_primitive_part(x, &cx);
-   dy=denom(y); if (!gcmp1(dy)) y=gmul(dy,y);
+   y = Q_primitive_part(y, &cy); cx = mul_content(cx,cy);
-   dx = mulii(dx,dy);
+   if (rx <= 2 || ry <= 2)
-   if (rx<=2 || ry<=2)
    {
-     m=cgetg(rx*ry+1,t_MAT);
+     m = cgetg(rx*ry+1,t_MAT);
      for (i=1; i<=rx; i++)
        for (j=1; j<=ry; j++)
-         m[(i-1)*ry+j]=(long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
+         m[(i-1)*ry+j] = (long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
      if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
      {
        GEN p1 = mulii(gcoeff(x,1,1),gcoeff(y,1,1));
        y = hnfmodid(m, p1);
      }
      else
-       y=hnfmod(m, detint(m));
+       y = hnfmod(m, detint(m));
    }
    else
    {
-     x=idealmat_to_hnf(nf,x);
+     if (!Z_ishnfall(x)) x = idealmat_to_hnf(nf,x);
-     y=idealmat_to_hnf(nf,y); y=idealmulh(nf,x,y);
+     if (!Z_ishnfall(y)) y = idealmat_to_hnf(nf,y);
+     y = idealmulh(nf,x,y);
    }
-   return gcmp1(dx)? y: gdiv(y,dx);
+   return cx? gmul(y,cx): y;
  }
- #if 0
+ /* operations on elements in factored form */
- /* y is principal */
- static GEN
+ GEN
- add_arch(GEN nf, GEN ax, GEN y)
+ to_famat(GEN g, GEN e)
  {
-   long tetpil, av=avma, prec=precision(ax);
+   GEN h = cgetg(3,t_MAT);
+   if (typ(g) != t_COL) { g = dummycopy(g); settyp(g, t_COL); }
+   if (typ(e) != t_COL) { e = dummycopy(e); settyp(e, t_COL); }
+   h[1] = (long)g;
+   h[2] = (long)e; return h;
+ }
-   y = get_arch(nf,y,prec); tetpil=avma;
+ GEN
-   return gerepile(av,tetpil,gadd(ax,y));
+ to_famat_all(GEN x, GEN y) { return to_famat(_col(x), _col(y)); }
+ static GEN
+ append(GEN v, GEN x)
+ {
+   long i, l = lg(v);
+   GEN w = cgetg(l+1, typ(v));
+   for (i=1; i<l; i++) w[i] = lcopy((GEN)v[i]);
+   w[i] = lcopy(x); return w;
  }
- #endif
  /* add x^1 to factorisation f */
  static GEN
-Line 1169  famat_add(GEN f, GEN x)
+Line 1224  famat_add(GEN f, GEN x)
 Line 1169  famat_add(GEN f, GEN x)
 Line 1224  famat_add(GEN f, GEN x)
    }
    else
    {
-     h[1] = (long)concat((GEN)f[1], x);
+     h[1] = (long)append((GEN)f[1], x); /* x may be a t_COL */
      h[2] = (long)concat((GEN)f[2], gun);
    }
    return h;
  }
  /* cf merge_factor_i */
- static GEN
+ GEN
  famat_mul(GEN f, GEN g)
  {
    GEN h;
-Line 1199  famat_sqr(GEN f)
+Line 1254  famat_sqr(GEN f)
 Line 1199  famat_sqr(GEN f)
 Line 1254  famat_sqr(GEN f)
    h[2] = lmul2n((GEN)f[2],1);
    return h;
  }
- static GEN
+ GEN
  famat_inv(GEN f)
  {
    GEN h;
-Line 1209  famat_inv(GEN f)
+Line 1264  famat_inv(GEN f)
 Line 1209  famat_inv(GEN f)
 Line 1264  famat_inv(GEN f)
    h[2] = lneg((GEN)f[2]);
    return h;
  }
- static GEN
+ GEN
  famat_pow(GEN f, GEN n)
  {
    GEN h;
    if (lg(f) == 1) return cgetg(1,t_MAT);
+   if (typ(f) != t_MAT) return to_famat_all(f,n);
    h = cgetg(3,t_MAT);
    h[1] = lcopy((GEN)f[1]);
    h[2] = lmul((GEN)f[2],n);
-Line 1235  famat_to_nf(GEN nf, GEN f)
+Line 1291  famat_to_nf(GEN nf, GEN f)
 Line 1235  famat_to_nf(GEN nf, GEN f)
 Line 1291  famat_to_nf(GEN nf, GEN f)
    return t;
  }
- GEN
+ /* "compare" two nf elt. Goal is to quickly sort for uniqueness of
- to_famat(GEN g, GEN e)
+  * representation, not uniqueness of represented element ! */
+ static int
+ elt_cmp(GEN x, GEN y)
  {
-   GEN h = cgetg(3,t_MAT);
+   long tx = typ(x), ty = typ(y);
-   h[1] = (long)g;
+   if (ty == tx)
-   h[2] = (long)e; return h;
+     return (tx == t_POL || tx == t_POLMOD)? cmp_pol(x,y): lexcmp(x,y);
+   return tx - ty;
  }
+ static int
+ elt_egal(GEN x, GEN y)
+ {
+   if (typ(x) == typ(y)) return gegal(x,y);
+   return 0;
+ }
  GEN
- to_famat_all(GEN x, GEN y) { return to_famat(_col(x), _col(y)); }
+ famat_reduce(GEN fa)
+ {
+   GEN E, F, G, L, g, e;
+   long i, k, l;
- /* assume (num(g[i]), id) = 1 and for all i. return prod g[i]^e[i] mod id */
+   if (lg(fa) == 1) return fa;
+   g = (GEN)fa[1]; l = lg(g);
+   e = (GEN)fa[2];
+   L = gen_sort(g, cmp_IND|cmp_C, &elt_cmp);
+   G = cgetg(l, t_COL);
+   E = cgetg(l, t_COL);
+   /* merge */
+   for (k=i=1; i<l; i++,k++)
+   {
+     G[k] = g[L[i]];
+     E[k] = e[L[i]];
+     if (k > 1 && elt_egal((GEN)G[k], (GEN)G[k-1]))
+     {
+       E[k-1] = laddii((GEN)E[k], (GEN)E[k-1]);
+       k--;
+     }
+   }
+   /* kill 0 exponents */
+   l = k;
+   for (k=i=1; i<l; i++)
+     if (!gcmp0((GEN)E[i]))
+     {
+       G[k] = G[i];
+       E[k] = E[i]; k++;
+     }
+   F = cgetg(3, t_MAT);
+   setlg(G, k); F[1] = (long)G;
+   setlg(E, k); F[2] = (long)E; return F;
+ }
+ /* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id */
  GEN
  famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id)
  {
-Line 1268  famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN
+Line 1366  famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN
 Line 1268  famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN
 Line 1366  famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN
    return t? t: gscalcol(gun, lg(id)-1);
  }
- /* assume prh has degree 1 and coprime to numerator(x) */
+ /* assume pr has degree 1 and coprime to numerator(x) */
- GEN
+ static GEN
- nf_to_Fp_simple(GEN x, GEN prh)
+ nf_to_Fp_simple(GEN x, GEN modpr, GEN p)
  {
-   GEN ch = denom(x), p = gcoeff(prh,1,1);
+   GEN c, r = zk_to_ff(Q_primitive_part(x, &c), modpr);
-   if (!is_pm1(ch))
+   if (c) r = gmod(gmul(r, c), p);
-   {
+   return r;
-     x = gmul(gmul(x,ch), mpinvmod(ch,p));
-   }
-   ch = colreducemodmat(gmod(x, p), prh, NULL);
-   return (GEN)ch[1]; /* in Fp^* */
  }
- GEN
+ static GEN
- famat_to_Fp_simple(GEN g, GEN e, GEN prh)
+ famat_to_Fp_simple(GEN nf, GEN x, GEN modpr, GEN p)
  {
-   GEN t = gun, h,n, p = gcoeff(prh,1,1), q = subis(p,1);
+   GEN h, n, t = gun, g = (GEN)x[1], e = (GEN)x[2], q = subis(p,1);
-   long i, lx = lg(g);
+   long i, l = lg(g);
-   for (i=1; i<lx; i++)
+   for (i=1; i<l; i++)
    {
      n = (GEN)e[i]; n = modii(n,q);
      if (!signe(n)) continue;
-     h = nf_to_Fp_simple((GEN)g[i], prh);
+     h = (GEN)g[i];
+     switch(typ(h))
+     {
+       case t_POL: case t_POLMOD: h = algtobasis(nf, h);  /* fall through */
+       case t_COL: h = nf_to_Fp_simple(h, modpr, p); break;
+       default: h = gmod(h, p);
+     }
      t = mulii(t, powmodulo(h, n, p)); /* not worth reducing */
    }
    return modii(t, p);
  }
- /* cf famat_to_nf_modideal_coprime, but id is a prime of degree 1 (=prh) */
+ /* cf famat_to_nf_modideal_coprime, but id is a prime of degree 1 (=pr) */
  GEN
- to_Fp_simple(GEN x, GEN prh)
+ to_Fp_simple(GEN nf, GEN x, GEN pr)
  {
+   GEN T, p, modpr = zk_to_ff_init(nf, &pr, &T, &p);
    switch(typ(x))
    {
-     case t_COL: return nf_to_Fp_simple(x,prh);
+     case t_COL: return nf_to_Fp_simple(x,modpr,p);
-     case t_MAT: return famat_to_Fp_simple((GEN)x[1],(GEN)x[2],prh);
+     case t_MAT: return famat_to_Fp_simple(nf,x,modpr,p);
      default: err(impl,"generic conversion to finite field");
    }
    return NULL;
  }
- extern GEN zinternallog_pk(GEN nf,GEN a0,GEN y,GEN pr,GEN prk,GEN list,GEN *psigne);
+ /* Compute t = prod g[i]^e[i] mod pr^k, assuming (t, pr) = 1.
- extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);
- extern GEN special_anti_uniformizer(GEN nf, GEN pr);
- extern GEN set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch);
- extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx, long v);
- /* Compute t = prod g[i]^e[i] mod pr^n, assuming (t, pr) = 1.
   * Method: modify each g[i] so that it becomes coprime to pr :
   *  x / (p^k u) --> x * (b/p)^v_pr(x) / z^k u, where z = b^e/p^(e-1)
   * b/p = vp^(-1) times something prime to p; both numerator and denominator
-Line 1324  extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, G
+Line 1421  extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, G
 Line 1324  extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, G
 Line 1421  extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, G
   * EX = exponent of (O_K / pr^k)^* used to reduce the product in case the
   * e[i] are large */
  static GEN
- famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prn, GEN EX)
+ famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX)
  {
-   long i,k, l = lg(g), N = degpol(nf[1]);
+   long i, l = lg(g);
-   GEN prnZ,cx,x,u,z, zpow = gzero, p = (GEN)pr[1], b = (GEN)pr[5];
+   GEN prkZ,cx,x,u, zpow = gzero, p = (GEN)pr[1], b = (GEN)pr[5];
-   GEN mul = cgetg(N+1,t_MAT);
+   GEN mul = eltmul_get_table(nf, b);
    GEN newg = cgetg(l+1, t_VEC); /* room for z */
-   prnZ = gcoeff(prn, 1,1);
+   prkZ = gcoeff(prk, 1,1);
-   z = gmod(special_anti_uniformizer(nf, pr), prnZ);
-   for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
    for (i=1; i < l; i++)
    {
      x = (GEN)g[i];
      if (typ(x) != t_COL) x = algtobasis(nf, x);
-     cx = denom(x); x = gmul(x,cx);
+     x = Q_remove_denom(x, &cx);
-     k = pvaluation(cx, p, &u);
+     if (cx)
-     if (!gcmp1(u)) /* could avoid the inversion, but prnZ is small--> cheap */
+     {
-       x = gmul(x, mpinvmod(u, prnZ));
+       long k = pvaluation(cx, p, &u);
-     if (k)
+       if (!gcmp1(u)) /* could avoid the inversion, but prkZ is small--> cheap */
-       zpow = addii(zpow, mulsi(k, (GEN)e[i]));
+         x = gmul(x, mpinvmod(u, prkZ));
-     (void)int_elt_val(nf, x, p, mul, &x, VERYBIGINT);
+       if (k)
-     newg[i] = (long)colreducemodmat(x, prn, NULL);
+         zpow = addii(zpow, mulsi(k, (GEN)e[i]));
+     }
+     (void)int_elt_val(nf, x, p, mul, &x);
+     newg[i] = (long)colreducemodHNF(x, prk, NULL);
    }
    if (zpow == gzero) setlg(newg, l);
    else
    {
-     newg[i] = (long)z;
+     newg[i] = (long)FpV_red(special_anti_uniformizer(nf, pr), prkZ);
      e = concatsp(e, negi(zpow));
    }
    e = gmod(e, EX);
-   return famat_to_nf_modideal_coprime(nf, newg, e, prn);
+   return famat_to_nf_modideal_coprime(nf, newg, e, prk);
  }
  GEN
  famat_ideallog(GEN nf, GEN g, GEN e, GEN bid)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
    GEN vp = gmael(bid, 3,1), ep = gmael(bid, 3,2), arch = gmael(bid,1,2);
    GEN cyc = gmael(bid,2,2), list_set = (GEN)bid[4], U = (GEN)bid[5];
    GEN p1,y0,x,y, psigne;
-   long i;
+   long i, l;
    if (lg(cyc) == 1) return cgetg(1,t_COL);
    y0 = y = cgetg(lg(U), t_COL);
    psigne = zsigne(nf, to_famat(g,e), arch);
-   for (i=1; i<lg(vp); i++)
+   l = lg(vp);
+   for (i=1; i < l; i++)
    {
      GEN pr = (GEN)vp[i], prk;
-     prk = idealpow(nf, pr, (GEN)ep[i]);
+     prk = (l==2)? gmael(bid,1,1): idealpow(nf, pr, (GEN)ep[i]);
      /* TODO: FIX group exponent (should be mod prk, not f !) */
      x = famat_makecoprime(nf, g, e, pr, prk, (GEN)cyc[1]);
      y = zinternallog_pk(nf, x, y, pr, prk, (GEN)list_set[i], &psigne);
-Line 1407  famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
+Line 1506  famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
 Line 1407  famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
 Line 1506  famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
  }
  GEN
+ famat_to_arch(GEN nf, GEN fa, long prec)
+ {
+   GEN g,e, y = NULL;
+   long i,l;
+   if (typ(fa) != t_MAT) return get_arch(nf, fa, prec);
+   if (lg(fa) == 1) return zerovec(lg(nf[6])-1);
+   g = (GEN)fa[1];
+   e = (GEN)fa[2]; l = lg(e);
+   for (i=1; i<l; i++)
+   {
+     GEN t, x = (GEN)g[i];
+     if (typ(x) != t_COL) x = algtobasis(nf,x);
+     x = Q_primpart(x);
+     /* multiplicative arch would be better (save logs), but exponents overflow
+      * [ could keep track of expo separately, but not worth it ] */
+     t = gmul(get_arch(nf,x,prec), (GEN)e[i]);
+     y = y? gadd(y,t): t;
+   }
+   return y;
+ }
+ GEN
  vecmul(GEN x, GEN y)
  {
    long i,lx, tx = typ(x);
-Line 1429  vecinv(GEN x)
+Line 1551  vecinv(GEN x)
 Line 1429  vecinv(GEN x)
 Line 1551  vecinv(GEN x)
  }
  GEN
+ vecpow(GEN x, GEN n)
+ {
+   long i,lx, tx = typ(x);
+   GEN z;
+   if (is_scalar_t(tx)) return powgi(x,n);
+   lx = lg(x); z = cgetg(lx, tx);
+   for (i=1; i<lx; i++) z[i] = (long)vecpow((GEN)x[i], n);
+   return z;
+ }
+ GEN
  vecdiv(GEN x, GEN y) { return vecmul(x, vecinv(y)); }
  /* x,y assumed to be of the same type, either
-Line 1450  GEN
+Line 1583  GEN
 Line 1450  GEN
 Line 1583  GEN
  arch_inv(GEN x) {
    switch (typ(x)) {
      case t_POLMOD: return ginv(x);
+     case t_COL:    return vecinv(x);
      case t_MAT:    return famat_inv(x);
-     default:       return gneg(x); /* t_COL, t_VEC */
+     default:       return gneg(x); /* t_VEC */
    }
  }
-Line 1459  GEN
+Line 1593  GEN
 Line 1459  GEN
 Line 1593  GEN
  arch_pow(GEN x, GEN n) {
    switch (typ(x)) {
      case t_POLMOD: return powgi(x,n);
+     case t_COL:    return vecpow(x,n);
      case t_MAT:    return famat_pow(x,n);
      default:       return gmul(n,x);
    }
-Line 1468  arch_pow(GEN x, GEN n) {
+Line 1603  arch_pow(GEN x, GEN n) {
 Line 1468  arch_pow(GEN x, GEN n) {
 Line 1603  arch_pow(GEN x, GEN n) {
  GEN
  idealmul(GEN nf, GEN x, GEN y)
  {
-   long tx,ty,av,f;
+   gpmem_t av;
+   long tx,ty,f;
    GEN res,ax,ay,p1;
    tx = idealtyp(&x,&ax);
-Line 1489  idealmul(GEN nf, GEN x, GEN y)
+Line 1625  idealmul(GEN nf, GEN x, GEN y)
 Line 1489  idealmul(GEN nf, GEN x, GEN y)
 Line 1625  idealmul(GEN nf, GEN x, GEN y)
            p1 = idealhermite_aux(nf, element_mul(nf,x,y));
            break;
          case id_PRIME:
-           p1 = gmul((GEN)y[1],x);
+         {
-           p1 = two_to_hnf(nf,p1, element_mul(nf,(GEN)y[2],x));
+           GEN mx = eltmul_get_table(nf, x);
+           GEN mpi= eltmul_get_table(nf, (GEN)y[2]);
+           p1 = concatsp(gmul(mx,(GEN)y[1]), gmul(mx,mpi));
+           p1 = idealmat_to_hnf(nf, p1);
            break;
+         }
          default: /* id_MAT */
-           p1 = idealmat_mul(nf,y, principalideal_aux(nf,x));
+           p1 = idealmat_to_hnf(nf, gmul(eltmul_get_table(nf,x), y));
        }break;
      case id_PRIME:
-Line 1518  idealmul(GEN nf, GEN x, GEN y)
+Line 1658  idealmul(GEN nf, GEN x, GEN y)
 Line 1518  idealmul(GEN nf, GEN x, GEN y)
 Line 1658  idealmul(GEN nf, GEN x, GEN y)
  GEN
  idealnorm(GEN nf, GEN x)
  {
-   long av = avma,tetpil;
+   gpmem_t av = avma,tetpil;
    GEN y;
    nf = checknf(nf);
-Line 1536  idealnorm(GEN nf, GEN x)
+Line 1676  idealnorm(GEN nf, GEN x)
 Line 1536  idealnorm(GEN nf, GEN x)
 Line 1676  idealnorm(GEN nf, GEN x)
  }
  /* inverse */
- extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
  /* rewritten from original code by P.M & M.H.
   *
-Line 1547  extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
+Line 1686  extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
 Line 1547  extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
 Line 1686  extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
  static GEN
  hnfideal_inv(GEN nf, GEN I)
  {
-   GEN J, dI = denom(I), IZ,dual;
+   GEN J, dI, IZ,dual;
-   if (gcmp1(dI)) dI = NULL; else I = gmul(I,dI);
+   I = Q_remove_denom(I, &dI);
    if (lg(I)==1) err(talker, "cannot invert zero ideal");
    IZ = gcoeff(I,1,1); /* I \cap Z */
    if (!signe(IZ)) err(talker, "cannot invert zero ideal");
-Line 1576  GEN
+Line 1715  GEN
 Line 1576  GEN
 Line 1715  GEN
  idealinv(GEN nf, GEN x)
  {
    GEN res,ax;
-   long av=avma, tx = idealtyp(&x,&ax);
+   gpmem_t av=avma;
+   long tx = idealtyp(&x,&ax);
    res = ax? cgetg(3,t_VEC): NULL;
    nf=checknf(nf); av=avma;
-Line 1595  idealinv(GEN nf, GEN x)
+Line 1735  idealinv(GEN nf, GEN x)
 Line 1595  idealinv(GEN nf, GEN x)
 Line 1735  idealinv(GEN nf, GEN x)
            case t_COL: x = gmul((GEN)nf[7],x); break;
            case t_POLMOD: x = (GEN)x[2]; break;
          }
-         x = ginvmod(x,(GEN)nf[1]);
+         x = QX_invmod(x,(GEN)nf[1]);
        }
        x = idealhermite_aux(nf,x); break;
      case id_PRIME:
-Line 1617  idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
+Line 1757  idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
 Line 1617  idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
 Line 1757  idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
    if (s < 0) n = negi(n);
    /* now n > 0 */
    x = dummycopy(vp);
-   n1 = dvmdii(n, (GEN)x[3], &r);
+   if (is_pm1(n)) /* n = 1 special cased for efficiency */
-   if (signe(r)) n1 = addis(n1,1); /* n1 = ceil(n/e) */
-   x[1] = (long)powgi((GEN)x[1],n1);
-   if (s < 0)
    {
-     x[2] = ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
+     if (s < 0)
-     *d = (GEN)x[1];
+     {
+       x[2] = x[5];
+       *d = (GEN)x[1];
+     }
+     else
+       *d = NULL;
    }
    else
    {
-     x[2] = (long)element_pow(nf,(GEN)x[2],n);
+     n1 = dvmdii(n, (GEN)x[3], &r);
-     *d = NULL;
+     if (signe(r)) n1 = addis(n1,1); /* n1 = ceil(n/e) */
+     x[1] = (long)powgi((GEN)x[1],n1);
+     if (s < 0)
+     {
+       x[2] = ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
+       *d = (GEN)x[1];
+     }
+     else
+     {
+       x[2] = (long)element_pow(nf,(GEN)x[2],n);
+       *d = NULL;
+     }
    }
    return x;
  }
-Line 1647  idealpowprime(GEN nf, GEN vp, GEN n)
+Line 1800  idealpowprime(GEN nf, GEN vp, GEN n)
 Line 1647  idealpowprime(GEN nf, GEN vp, GEN n)
 Line 1800  idealpowprime(GEN nf, GEN vp, GEN n)
    return x;
  }
- /* x * vp^n */
+ /* x * vp^n. Assume x in HNF (possibly non-integral) */
  GEN
  idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n)
  {
-   GEN denx,y,d;
+   GEN cx,y,dx;
    if (!signe(n)) return x;
    nf = checknf(nf);
-   y = idealpowprime_spec(nf, vp, n, &d);
-   denx = denom(x);
+   /* inert, special cased for efficiency */
-   if (gcmp1(denx)) denx = d; else
+   if (itos((GEN)vp[4]) == degpol(nf[1]))
+     return gmul(x, powgi((GEN)vp[1], n));
+   y = idealpowprime_spec(nf, vp, n, &dx);
+   x = Q_primitive_part(x, &cx);
+   if (cx && dx)
    {
-     x = gmul(denx,x);
+     cx = gdiv(cx, dx);
-     if (d) denx = mulii(d,denx);
+     if (typ(cx) != t_FRAC) dx = NULL;
+     else { dx = (GEN)cx[2]; cx = (GEN)cx[1]; }
+     if (is_pm1(cx)) cx = NULL;
    }
-   x = idealmulspec(nf,x, (GEN)y[1], (GEN)y[2]);
+   x = idealmulspec(nf,x,y);
-   return denx? gdiv(x,denx): x;
+   if (cx) x = gmul(x,cx);
+   if (dx) x = gdiv(x,dx);
+   return x;
  }
+ GEN
+ idealdivpowprime(GEN nf, GEN x, GEN vp, GEN n)
+ {
+   return idealmulpowprime(nf,x,vp, negi(n));
+ }
  /* raise the ideal x to the power n (in Z) */
  GEN
  idealpow(GEN nf, GEN x, GEN n)
  {
-   long tx,N,av,s,i;
+   gpmem_t av;
+   long tx,N,s;
    GEN res,ax,m,cx,n1,a,alpha;
    if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpow");
-Line 1696  idealpow(GEN nf, GEN x, GEN n)
+Line 1864  idealpow(GEN nf, GEN x, GEN n)
 Line 1696  idealpow(GEN nf, GEN x, GEN n)
 Line 1864  idealpow(GEN nf, GEN x, GEN n)
        default:
          n1 = (s<0)? negi(n): n;
-         cx = content(x); if (gcmp1(cx)) cx = NULL; else x = gdiv(x,cx);
+         x = Q_primitive_part(x, &cx);
          a=ideal_two_elt(nf,x); alpha=(GEN)a[2]; a=(GEN)a[1];
-         m = cgetg(N+1,t_MAT); a = powgi(a,n1);
          alpha = element_pow(nf,alpha,n1);
-         for (i=1; i<=N; i++) m[i]=(long)element_mulid(nf,alpha,i);
+         m = eltmul_get_table(nf, alpha);
-         x = hnfmodid(m, a);
+         x = hnfmodid(m, powgi(a,n1));
          if (s<0) x = hnfideal_inv(nf,x);
          if (cx) x = gmul(x, powgi(cx,n));
      }
-Line 1717  idealpow(GEN nf, GEN x, GEN n)
+Line 1884  idealpow(GEN nf, GEN x, GEN n)
 Line 1717  idealpow(GEN nf, GEN x, GEN n)
 Line 1884  idealpow(GEN nf, GEN x, GEN n)
  GEN
  idealpows(GEN nf, GEN ideal, long e)
  {
-   long court[] = {evaltyp(t_INT) | m_evallg(3),0,0};
+   long court[] = {evaltyp(t_INT) | _evallg(3),0,0};
    affsi(e,court); return idealpow(nf,ideal,court);
  }
- GEN
+ static GEN
- init_idele(GEN nf)
+ _idealmulred(GEN nf, GEN x, GEN y, long prec)
  {
-   GEN x = cgetg(3,t_VEC);
+   return ideallllred(nf,idealmul(nf,x,y), NULL, prec);
-   long RU;
-   nf = checknf(nf); RU = lg(nf[6])-1;
-   x[2] = (long)zerovec(RU); return x;
  }
+ typedef struct {
+   GEN nf;
+   long prec;
+ } idealred_muldata;
+ static GEN
+ _mul(void *data, GEN x, GEN y)
+ {
+   idealred_muldata *D = (idealred_muldata *)data;
+   return _idealmulred(D->nf,x,y,D->prec);
+ }
+ static GEN
+ _sqr(void *data, GEN x)
+ {
+   return _mul(data,x,x);
+ }
  /* compute x^n (x ideal, n integer), reducing along the way */
  GEN
  idealpowred(GEN nf, GEN x, GEN n, long prec)
  {
-   long i,j,m,av=avma, s = signe(n);
+   gpmem_t av=avma;
-   GEN y, p1;
+   long s = signe(n);
+   idealred_muldata D;
+   GEN y;
    if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpowred");
-   if (signe(n) == 0) return idealpow(nf,x,n);
+   if (s == 0) return idealpow(nf,x,n);
-   p1 = n+2; m = *p1;
+   D.nf  = nf;
-   y = x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
+   D.prec= prec;
-   for (i=lgefint(n)-2;;)
+   y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul);
-   {
-     for (; j; m<<=1,j--)
-     {
-       y = idealmul(nf,y,y);
-       if (m < 0) y = idealmul(nf,y,x);
-       y = ideallllred(nf,y,NULL,prec);
-     }
-     if (--i == 0) break;
-     m = *++p1; j = BITS_IN_LONG;
-   }
    if (s < 0) y = idealinv(nf,y);
-   if (y == x) y = ideallllred(nf,x,NULL,prec);
+   if (s < 0 || is_pm1(n)) y = ideallllred(nf,y,NULL,prec);
    return gerepileupto(av,y);
  }
  GEN
  idealmulred(GEN nf, GEN x, GEN y, long prec)
  {
-   long av=avma;
+   gpmem_t av = avma;
-   x = idealmul(nf,x,y);
+   return gerepileupto(av, _idealmulred(nf,x,y,prec));
-   return gerepileupto(av, ideallllred(nf,x,NULL,prec));
  }
  long
  isideal(GEN nf,GEN x)
  {
-   long N,av,i,j,k,tx=typ(x),lx;
+   gpmem_t av;
-   GEN p1,minv;
+   long N,i,j,tx=typ(x),lx;
    nf=checknf(nf); lx=lg(x);
    if (tx==t_VEC && lx==3) { x=(GEN)x[1]; tx=typ(x); lx=lg(x); }
-Line 1779  isideal(GEN nf,GEN x)
+Line 1952  isideal(GEN nf,GEN x)
 Line 1779  isideal(GEN nf,GEN x)
 Line 1952  isideal(GEN nf,GEN x)
    if (typ(x)==t_VEC) return (lx==6);
    if (typ(x)!=t_MAT) return 0;
    if (lx == 1) return 1;
-   N=lgef(nf[1])-2; if (lg(x[1]) != N) return 0;
+   N = degpol(nf[1]); if (lg(x[1])-1 != N) return 0;
-   av=avma;
+   av = avma;
-   if (lx != N) x = idealmat_to_hnf(nf,x);
+   if (lx-1 != N) x = idealmat_to_hnf(nf,x);
-   x = gdiv(x,content(x)); minv=ginv(x);
+   x = primpart(x);
-   for (i=1; i<N; i++)
+   for (i=1; i<=N; i++)
-     for (j=1; j<N; j++)
+     for (j=2; j<=N; j++)
-     {
+       if (! hnf_invimage(x, element_mulid(nf,(GEN)x[i],j)))
-       p1=gmul(minv, element_mulid(nf,(GEN)x[i],j));
+       {
-       for (k=1; k<N; k++)
+         avma = av; return 0;
-         if (typ(p1[k])!=t_INT) { avma=av; return 0; }
+       }
-     }
    avma=av; return 1;
  }
  GEN
  idealdiv(GEN nf, GEN x, GEN y)
  {
-   long av=avma,tetpil;
+   gpmem_t av=avma,tetpil;
    GEN z=idealinv(nf,y);
    tetpil=avma; return gerepile(av,tetpil,idealmul(nf,x,z));
-Line 1826  idealdiv(GEN nf, GEN x, GEN y)
+Line 1998  idealdiv(GEN nf, GEN x, GEN y)
 Line 1826  idealdiv(GEN nf, GEN x, GEN y)
 Line 1998  idealdiv(GEN nf, GEN x, GEN y)
  GEN
  idealdivexact(GEN nf, GEN x0, GEN y0)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
    GEN x,y,Nx,Ny,Nz, cy = content(y0);
    nf = checknf(nf);
-Line 1860  idealdivexact(GEN nf, GEN x0, GEN y0)
+Line 2032  idealdivexact(GEN nf, GEN x0, GEN y0)
 Line 1860  idealdivexact(GEN nf, GEN x0, GEN y0)
 Line 2032  idealdivexact(GEN nf, GEN x0, GEN y0)
  GEN
  idealintersect(GEN nf, GEN x, GEN y)
  {
-   long av=avma,lz,i,N;
+   gpmem_t av=avma;
+   long lz,i,N;
    GEN z,dx,dy;
    nf=checknf(nf); N=degpol(nf[1]);
-Line 1877  idealintersect(GEN nf, GEN x, GEN y)
+Line 2050  idealintersect(GEN nf, GEN x, GEN y)
 Line 1877  idealintersect(GEN nf, GEN x, GEN y)
 Line 2050  idealintersect(GEN nf, GEN x, GEN y)
    return gerepileupto(av,z);
  }
- static GEN
+ /*******************************************************************/
- computet2twist(GEN nf, GEN vdir)
+ /*                                                                 */
+ /*                      T2-IDEAL REDUCTION                         */
+ /*                                                                 */
+ /*******************************************************************/
+ GEN
+ computeGtwist(GEN nf, GEN vdir)
  {
-   long j, ru = lg(nf[6]);
+   long i, j, k, l, lG, v, r1, r2;
-   GEN p1,MC, mat = (GEN)nf[5];
+   GEN p1, G = gmael(nf,5,2);
-   if (!vdir) return (GEN)mat[3];
+   if (!vdir) return G;
-   MC=(GEN)mat[2]; p1=cgetg(ru,t_MAT);
+   l = lg(vdir); lG = lg(G);
-   for (j=1; j<ru; j++)
+   p1 = dummycopy(G);
+   nf_get_sign(nf, &r1, &r2);
+   for (i=1; i<l; i++)
    {
-     GEN v = (GEN)vdir[j];
+     v = vdir[i]; if (!v) continue;
-     if (gcmp0(v))
+     if (i <= r1) {
-       p1[j] = MC[j];
+       for (j=1; j<lG; j++) coeff(p1,i,j) = lmul2n(gcoeff(p1,i,j), v);
-     else if (typ(v) == t_INT)
+     } else {
-       p1[j] = lmul2n((GEN)MC[j],itos(v)<<1);
+       k = (i<<1) - r1;
-     else
+       for (j=1; j<lG; j++)
-       p1[j] = lmul((GEN)MC[j],gpui(stoi(4),v,0));
+       {
+         coeff(p1,k-1,j) = lmul2n(gcoeff(p1,k-1,j), v);
+         coeff(p1,k  ,j) = lmul2n(gcoeff(p1,k  ,j), v);
+       }
+     }
    }
-   return mulmat_real(p1,(GEN)mat[1]);
+   return p1;
  }
  static GEN
  chk_vdir(GEN nf, GEN vdir)
  {
+   long l, i, t;
+   GEN v;
    if (!vdir || gcmp0(vdir)) return NULL;
-   if (typ(vdir)!=t_VEC || lg(vdir) != lg(nf[6])) err(idealer5);
+   l = lg(vdir);
-   return vdir;
+   if (l != lg(nf[6])) err(talker, "incorrect vector length in idealred");
+   t = typ(vdir);
+   if (t == t_VECSMALL) return vdir;
+   if (t != t_VEC) err(talker, "not a vector in idealred");
+   v = cgetg(l, t_VECSMALL);
+   for (i=1; i<l; i++) v[i] = itos(gceil((GEN)vdir[i]));
+   return v;
  }
  /* assume I in NxN matrix form (not necessarily HNF) */
  static GEN
- ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *ptprec)
+ idealred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *ptprec)
  {
-   GEN T2, u, y, nf = *ptnf;
+   GEN G, u, y, nf = *ptnf;
    long i, e, prec = *ptprec;
+   e = gexpo(I)>>TWOPOTBITS_IN_LONG;
+   if (e < 0) e = 0;
    for (i=1; ; i++)
    {
-     T2 = computet2twist(nf,vdir);
+     G = computeGtwist(nf,vdir);
-     y = qf_base_change(T2,I,1);
+     y = gmul(G, I);
-     e = 1 + (gexpo(y)>>TWOPOTBITS_IN_LONG);
+     u = lllintern(y, 100, 1, prec);
-     if (e < 0) e = 0;
-     u = lllgramintern(y,100,1, e + prec);
      if (u) break;
      if (i == MAXITERPOL) err(accurer,"ideallllred");
      prec = (prec<<1)-2;
      if (DEBUGLEVEL) err(warnprec,"ideallllred",prec);
-     nf = nfnewprec(nf, (e>>1)+prec);
+     nf = nfnewprec(nf, e + prec);
    }
    *ptprec = prec;
    *ptnf = nf;
-Line 1933  ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *pt
+Line 2126  ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *pt
 Line 1933  ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *pt
 Line 2126  ideallllred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *pt
  }
  GEN
- ideallllred_elt(GEN nf, GEN I)
+ ideallllred_elt(GEN nf, GEN I, GEN vdir)
  {
    long prec = DEFAULTPREC;
-   return ideallllred_elt_i(&nf, I, NULL, &prec);
+   return idealred_elt_i(&nf, I, vdir, &prec);
  }
+ GEN
+ idealred_elt(GEN nf, GEN I)
+ {
+   return ideallllred_elt(nf, I, NULL);
+ }
  GEN
  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
    long N,i,nfprec;
    GEN J,I0,Ired,res,aI,y,x,Nx,b,c1,c,pol;
-Line 1959  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
+Line 2157  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
 Line 1959  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
 Line 2157  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
    if (DEBUGLEVEL>5) msgtimer("entering idealllred");
    I0 = I;
    if (typ(I) != id_MAT || lg(I) != N+1) I = idealhermite_aux(nf,I);
-   c1 = content(I); if (gcmp1(c1)) c1 = NULL; else I = gdiv(I,c1);
+   I = primitive_part(I, &c1);
    if (2 * expi(gcoeff(I,1,1)) >= bit_accuracy(nfprec))
-     Ired = gmul(I, lllintpartial(I));
+     Ired = lllint_ip(I,4);
    else
      Ired = I;
-   y = ideallllred_elt_i(&nf, Ired, chk_vdir(nf,vdir), &prec);
+   y = idealred_elt_i(&nf, Ired, chk_vdir(nf,vdir), &prec);
    if (isnfscalar(y))
    { /* already reduced */
-Line 1974  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
+Line 2172  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
 Line 1974  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
 Line 2172  ideallllred(GEN nf, GEN I, GEN vdir, long prec)
    if (DEBUGLEVEL>5) msgtimer("LLL reduction");
    x = gmul((GEN)nf[7], y); Nx = subres(pol,x);
-   b = gmul(Nx, ginvmod(x,pol));
+   b = gmul(Nx, QX_invmod(x,pol));
-   b = algtobasis_intern(nf,b);
+   b = algtobasis_i(nf,b);
    J = cgetg(N+1,t_MAT); /* = I Nx / x integral */
    for (i=1; i<=N; i++)
      J[i] = (long)element_muli(nf,b,(GEN)Ired[i]);
-   c = content(J); if (!gcmp1(c)) J = gdiv(J,c);
+   J = Q_primitive_part(J, &c);
   /* c = content (I Nx / x) = Nx / den(I/x) --> d = den(I/x) = Nx / c
    * J = (d I / x); I[1,1] = I \cap Z --> d I[1,1] belongs to J and Z */
    if (isnfscalar((GEN)I[1]))
-     b = mulii(gcoeff(I,1,1), divii(Nx, c));
+     b = mulii(gcoeff(I,1,1), c? diviiexact(Nx, c): Nx);
    else
      b = detint(J);
    I = hnfmodid(J,b);
-Line 1994  END:
+Line 2192  END:
 Line 1994  END:
 Line 2192  END:
    switch(typ(aI))
    {
-     case t_POLMOD: case t_MAT: /* compute y, I0 = J y */
+     case t_POLMOD: /* compute y, I0 = J y */
        if (!Nx) y = c1;
        else
        {
-         if (c1) c = gmul(c,c1);
+         c = mul_content(c,c1);
-         y = gmul(x, gdiv(c,Nx));
+         y = c? gmul(x, gdiv(c,Nx)): gdiv(x, Nx);
        }
        break;
+     case t_MAT: /* compute y, I0 = J y */
+       if (!Nx) y = c1;
+       else
+       {
+         c = mul_content(c,c1);
+         y = c? gmul(y, gdiv(c,Nx)): gdiv(y, Nx);
+       }
+       break;
      case t_COL:
        if (y) y = vecinv(gmul(gmael(nf,5,1), y));
        break;
-Line 2020  END:
+Line 2227  END:
 Line 2020  END:
 Line 2227  END:
  GEN
  minideal(GEN nf, GEN x, GEN vdir, long prec)
  {
-   long av = avma, N, tx;
+   gpmem_t av = avma;
-   GEN p1,y;
+   long N, tx;
+   GEN y;
    nf = checknf(nf);
    vdir = chk_vdir(nf,vdir);
-Line 2030  minideal(GEN nf, GEN x, GEN vdir, long prec)
+Line 2238  minideal(GEN nf, GEN x, GEN vdir, long prec)
 Line 2030  minideal(GEN nf, GEN x, GEN vdir, long prec)
 Line 2238  minideal(GEN nf, GEN x, GEN vdir, long prec)
    if (tx == id_PRINCIPAL) return gcopy(x);
    if (tx != id_MAT || lg(x) != N+1) x = idealhermite_aux(nf,x);
-   p1 = computet2twist(nf,vdir);
+   y = gmul(computeGtwist(nf,vdir), x);
-   y = qf_base_change(p1,x,0);
+   y = gmul(x, (GEN)lll(y,prec)[1]);
-   y = gmul(x, (GEN)lllgram(y,prec)[1]);
    return gerepileupto(av, principalidele(nf,y,prec));
  }
+ /*******************************************************************/
+ /*                                                                 */
+ /*                   APPROXIMATION THEOREM                         */
+ /*                                                                 */
+ /*******************************************************************/
+ /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
  static GEN
- appr_reduce(GEN s, GEN y, long N)
+ coprime_part(GEN x, GEN f)
  {
-   GEN p1,u,z = cgetg(N+2,t_MAT);
+   for (;;)
-   long i;
+   {
+     f = mppgcd(x, f); if (is_pm1(f)) break;
-   s=gmod(s,gcoeff(y,1,1)); y=gmul(y,lllint(y));
+     x = diviiexact(x, f);
-   for (i=1; i<=N; i++) z[i]=y[i]; z[N+1]=(long)s;
+   }
-   u=(GEN)ker(z)[1]; p1 = denom(u);
+   return x;
-   if (!gcmp1(p1)) u=gmul(u,p1);
-   p1=(GEN)u[N+1]; setlg(u,N+1);
-   for (i=1; i<=N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
-   return gadd(s, gmul(y,u));
  }
- /* Given a fractional ideal x (if fl=0) or a prime ideal factorization
+ /* x t_INT, f ideal. Write x = x1 x2, sqf(x1) | f, (x2,f) = 1. Return x2 */
-  * with possibly zero or negative exponents (if fl=1), gives a b such that
+ static GEN
-  * v_p(b)=v_p(x) for all prime ideals p dividing x (or in the factorization)
+ nf_coprime_part(GEN nf, GEN x, GEN *listpr)
-  * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
-  * Certainly not the most efficient, but sure.
-  */
- GEN
- idealappr0(GEN nf, GEN x, long fl)
  {
-   long av=avma,tetpil,i,j,k,l,N,r,r2;
+   long v, j, lp = lg(listpr), N = degpol(nf[1]);
-   GEN fact,fact2,list,ep,ep1,ep2,y,z,v,p1,p2,p3,p4,s,pr,alpha,beta,den;
+   GEN x1, x2, ex, p, pr;
-   if (DEBUGLEVEL>4)
+ #if 0 /*1) via many gcds. Expensive ! */
+   GEN f = idealprodprime(nf, (GEN)listpr);
+   f = hnfmodid(f, x); /* first gcd is less expensive since x in Z */
+   x = gscalmat(x, N);
+   for (;;)
    {
-     fprintferr(" entree dans idealappr0() :\n");
+     if (gcmp1(gcoeff(f,1,1))) break;
-     fprintferr(" x = "); outerr(x);
+     x = idealdivexact(nf, x, f);
+     f = hnfmodid(concatsp(f,x), gcoeff(x,1,1)); /* gcd(f,x) */
    }
-   nf=checknf(nf); N=degpol(nf[1]);
+   x2 = x;
-   if (fl)
+ #else /*2) from prime decomposition */
+   x1 = NULL;
+   for (j=1; j<lp; j++)
    {
-     if (typ(x)!=t_MAT || lg(x)!=3)
+     pr = listpr[j]; p = (GEN)pr[1];
-       err(talker,"not a prime ideal factorization in idealappr0");
+     v = ggval(x, p); if (!v) continue;
-     fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
-     if (r==1) return gscalcol_i(gun,N);
+     ex = mulsi(v, (GEN)pr[3]); /* = v_pr(x) > 0 */
-     for (i=1; i<r; i++)
+     x1 = x1? idealmulpowprime(nf, x1, pr, ex)
-       if (signe(ep[i]) < 0) break;
+            : idealpow(nf, pr, ex);
-     if (i < r)
-     {
-       ep1=cgetg(r,t_COL);
-       for (i=1; i<r; i++)
-         ep1[i] = (signe(ep[i])>=0)? zero: lnegi((GEN)ep[i]);
-       fact[2]=(long)ep1; beta=idealappr0(nf,fact,1);
-       fact2=idealfactor(nf,beta);
-       p1=(GEN)fact2[1]; r2=lg(p1);
-       ep2=(GEN)fact2[2]; l=r+r2-1;
-       z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
-       ep1=cgetg(l,t_VEC);
-       for (i=1; i<r; i++)
-         ep1[i] = (signe(ep[i])<=0)? zero: licopy((GEN)ep[i]);
-       j=r-1;
-       for (i=1; i<r2; i++)
-       {
-         p3=(GEN)p1[i]; k=1;
-         while (k<r &&
-           (    !gegal((GEN)p3[1],gmael(list,k,1))
-             || !element_val(nf,(GEN)p3[2],(GEN)list[k]) )) k++;
-         if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
-       }
-       fact=cgetg(3,t_MAT);
-       fact[1]=(long)z; setlg(z,j+1);
-       fact[2]=(long)ep1; setlg(ep1,j+1);
-       alpha=idealappr0(nf,fact,1); tetpil=avma;
-       if (DEBUGLEVEL>2)
-       {
-         fprintferr(" alpha = "); outerr(alpha);
-         fprintferr(" beta = "); outerr(beta);
-       }
-       return gerepile(av,tetpil,element_div(nf,alpha,beta));
-     }
-     y=idmat(N);
-     for (i=1; i<r; i++)
-     {
-       pr=(GEN)list[i];
-       if (signe(ep[i]))
-       {
-         p4=addsi(1,(GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
-         if (cmpis((GEN)pr[4],N))
-         {
-           p2=cgetg(3,t_MAT);
-           p2[1]=(long)gscalcol_i(p1, N);
-           p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
-           y=idealmat_mul(nf,y,p2);
-         }
-         else y=gmul(p1,y);
-       }
-       else y=idealmulprime(nf,y,pr);
-     }
    }
+   x = gscalmat(x, N);
+   x2 = x1? idealdivexact(nf, x, x1): x;
+ #endif
+   return x2;
+ }
+ /* assume L is a list of prime ideals. Return the product */
+ GEN
+ idealprodprime(GEN nf, GEN L)
+ {
+   long l = lg(L), i;
+   GEN z;
+   if (l == 1) return idmat(degpol(nf[1]));
+   z = prime_to_ideal(nf, (GEN)L[1]);
+   for (i=2; i<l; i++) z = idealmulprime(nf,z, (GEN)L[i]);
+   return z;
+ }
+ /* assume L is a list of prime ideals. Return prod L[i]^e[i] */
+ GEN
+ factorbackprime(GEN nf, GEN L, GEN e)
+ {
+   long l = lg(L), i;
+   GEN z;
+   if (l == 1) return idmat(degpol(nf[1]));
+   z = idealpow(nf, (GEN)L[1], (GEN)e[1]);
+   for (i=2; i<l; i++)
+     if (signe(e[i])) z = idealmulpowprime(nf,z, (GEN)L[i],(GEN)e[i]);
+   return z;
+ }
+ /* compute anti-uniformizer for pr, coprime to f outside of pr, integral
+  * outside of p below pr.
+  * sqf = product or primes dividing f, NULL if f a prime power*/
+ GEN
+ anti_unif_mod_f(GEN nf, GEN pr, GEN sqf)
+ {
+   GEN U, V, UV, uv, cx, t, p = (GEN)pr[1];
+   if (!sqf) return gdiv((GEN)pr[5], p);
    else
    {
-     den=denom(x); if (gcmp1(den)) den=NULL; else x=gmul(den,x);
+     GEN d,d1,d2, sqfZ = gcoeff(sqf,1,1); /* = (sqf \cap Z) = (V \cap Z) * p */
-     x=idealhermite_aux(nf,x);
+     U = idealpow(nf,pr,gdeux);
-     fact=idealfactor(nf,x);
+     V = idealdivpowprime(nf,sqf,pr,gun);
-     list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
+     uv = addone_nored(U, V);
-     if (r==1) { avma=av; return gscalcol_i(gun,N); }
+     UV = idealmul(nf,sqf,pr);
-     if (den)
+     t = (GEN)pr[2];
-     {
+     t = makeprimetoideal(nf, UV, uv, t); /* cf unif_mod_f */
-       fact2=idealfactor(nf,den);
-       p1=(GEN)fact2[1]; r2=lg(p1);
+     t = element_inv(nf, t);
-       l=r+r2-1;
+     t = primitive_part(t, &cx); /* p | denom(cx) */
-       z=cgetg(l,t_COL);   for (i=1; i<r; i++) z[i]=list[i];
+     d = denom(cx);
-       ep1=cgetg(l,t_COL); for (i=1; i<r; i++) ep1[i]=ep[i];
+     d2 = coprime_part(d, sqfZ);
-       j=r-1;
+     d1 = diviiexact(d, d2); /* p | d1 */
-       for (i=1; i<r2; i++)
+     cx = gmod(gmul(cx,d1), mulii(d1, gcoeff(V,1,1)));
-       {
+     t = gdiv(colreducemodHNF(gmul(cx,t), gmul(d1,V), NULL), d1);
-         p3=(GEN)p1[i]; k=1;
+     return t; /* v_pr[i](t) = -1, v_pr[j](t) = 0 if i != j */
-         while (k<r && !gegal((GEN)list[k],p3)) k++;
-         if (k==r){ j++; z[j]=(long)p3; ep1[j]=zero; }
-       }
-       fact=cgetg(3,t_MAT);
-       fact[1]=(long)z; setlg(z,j+1);
-       fact[2]=(long)ep1; setlg(ep1,j+1);
-       alpha=idealappr0(nf,fact,1);
-       if (DEBUGLEVEL>2) { fprintferr(" alpha = "); outerr(alpha); }
-       tetpil=avma; return gerepile(av,tetpil,gdiv(alpha,den));
-     }
-     y=x; for (i=1; i<r; i++) y=idealmulprime(nf,y,(GEN)list[i]);
    }
+ }
-   z=cgetg(r,t_VEC);
+ /* compute integral uniformizer for pr, coprime to f outside of pr.
+  * sqf = product or primes dividing f, NULL if f a prime power*/
+ GEN
+ unif_mod_f(GEN nf, GEN pr, GEN sqf)
+ {
+   GEN U, V, UV, uv, t = (GEN)pr[2];
+   if (!sqf) return t;
+   else
+   {
+     U = idealpow(nf,pr,gdeux);
+     V = idealdivpowprime(nf,sqf,pr,gun);
+     uv = addone_nored(U, V);
+     UV = idealmulprime(nf,sqf,pr);
+     return makeprimetoideal(nf, UV, uv, t);
+   }
+ }
+ static GEN
+ appr_reduce(GEN s, GEN y)
+ {
+   long i, N = lg(y)-1;
+   GEN d, u, z = cgetg(N+2,t_MAT);
+   s = gmod(s, gcoeff(y,1,1));
+   y = lllint_ip(y,100);
+   for (i=1; i<=N; i++) z[i] = y[i];
+   z[N+1] = (long)s;
+   u = (GEN)ker(z)[1];
+   u = Q_remove_denom(u, NULL);
+   d = (GEN)u[N+1]; setlg(u,N+1);
+   for (i=1; i<=N; i++) u[i] = (long)diviiround((GEN)u[i],d);
+   return gadd(s, gmul(y,u));
+ }
+ /* L0 in K^*.
+  * If ptd1 == NULL, assume (L0,f) = 1
+  *   return L integral, L0 = L mod f
+  *
+  * Otherwise, assume v_pr(L0) <= 0 for all pr | f and set *ptd1 = d1
+  *   return L integral, L0 = L/d1 mod f, and such that
+  *   if (L*I,f) = 1 for some integral I, then d1 | L*I  */
+ GEN
+ make_integral(GEN nf, GEN L0, GEN f, GEN *listpr, GEN *ptd1)
+ {
+   GEN fZ, t, L, D2, d1, d2, d;
+   if (ptd1) *ptd1 = NULL;
+   L = Q_remove_denom(L0, &d);
+   if (!d) return L0;
+   /* L0 = L / d, L integral */
+   fZ = gcoeff(f,1,1);
+   /* Kill denom part coprime to fZ */
+   d2 = coprime_part(d, fZ);
+   t = mpinvmod(d2, fZ); if (!is_pm1(t)) L = gmul(L,t);
+   if (egalii(d, d2)) return L;
+   d1 = diviiexact(d, d2);
+   /* L0 = (L / d1) mod f. d1 not coprime to f
+    * write (d1) = D1 D2, D2 minimal, (D2,f) = 1. */
+   D2 = nf_coprime_part(nf, d1, listpr);
+   t = idealaddtoone_i(nf, D2, f); /* in D2, 1 mod f */
+   L = element_muli(nf,t,L);
+   /* if (L0, f) = 1, then L in D1 ==> in D1 D2 = (d1) */
+   if (!ptd1) return Q_div_to_int(L, d1); /* exact division */
+   *ptd1 = d1; return L;
+ }
+ void
+ check_listpr(GEN x)
+ {
+   long l = lg(x), i;
+   for (i=1; i<l; i++) checkprimeid((GEN)x[i]);
+ }
+ /* Given a prime ideal factorization with possibly zero or negative
+  * exponents, gives b such that v_p(b) = v_p(x) for all prime ideals pr | x
+  * and v_pr(b)> = 0 for all other pr.
+  * For optimal performance, all [anti-]uniformizers should be precomputed,
+  * but no support for this yet.
+  * No garbage collecting */
+ GEN
+ idealapprfact_i(GEN nf, GEN x)
+ {
+   gpmem_t av = avma;
+   GEN tau, pi, z, d, sqf, sqfsafe, list, e, e2;
+   long s, flag, i, r, N;
+   nf = checknf(nf);
+   if (typ(x) != t_MAT || lg(x) != 3)
+     err(talker,"not a factorization in idealapprfact");
+   check_listpr((GEN)x[1]);
+   if (DEBUGLEVEL>4) {
+     fprintferr(" entering idealapprfact():\n");
+     fprintferr(" x = %Z\n", x);
+   }
+   list= (GEN)x[1];
+   e   = (GEN)x[2]; r = lg(e); N = degpol(nf[1]);
+   if (r==1) return gscalcol_i(gun, N);
+   sqf = (r == 2)? NULL: idealprodprime(nf, list);
+   sqfsafe = NULL;
+   z = gun; flag = 0;
    for (i=1; i<r; i++)
    {
-     pr=(GEN)list[i]; p4=addsi(1, (GEN)ep[i]); p1=powgi((GEN)pr[1],p4);
+     s = signe(e[i]);
-     if (cmpis((GEN)pr[4],N))
+     if (s < 0)
      {
-       p2=cgetg(3,t_MAT);
+       flag = 1;
-       p2[1]=(long)gscalcol_i(p1,N);
+       if (!sqfsafe) sqfsafe = sqf? sqf: prime_to_ideal(nf, (GEN)list[1]);
-       p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
+       tau = anti_unif_mod_f(nf, (GEN)list[i], sqf);
-       z[i]=ldiv(idealmat_mul(nf,y,p2),p1);
+       tau = make_integral(nf, tau, sqfsafe, (GEN*)list, &d);
+       tau = gdiv(tau, d);
+       z = element_mul(nf, z, element_pow(nf, tau, negi((GEN)e[i])));
      }
-     else z[i]=ldiv(y,p1);
+     else if (s > 0)
+     {
+       pi = unif_mod_f(nf, (GEN)list[i], sqf);
+       z = element_mul(nf, z, element_pow(nf, pi, (GEN)e[i]));
+     }
    }
-   v=idealaddmultoone(nf,z);
+   if (z == gun) { avma = av; return gscalcol_i(gun, N); }
-   s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
+   e2 = cgetg(r, t_VEC);
-   for (i=1; i<r; i++)
+   for (i=1; i<r; i++) e2[i] = laddis((GEN)e[i], 1);
+   x = factorbackprime(nf, list,e2);
+   if (flag) /* denominator */
    {
-     pr=(GEN)list[i];
+     z = Q_remove_denom(z, &d);
-     if (signe(ep[i]))
+     x = Q_muli_to_int(x, d);
-       s=gadd(s,element_mul(nf,(GEN)v[i],element_pow(nf,(GEN)pr[2],(GEN)ep[i])));
-     else s=gadd(s,(GEN)v[i]);
    }
-   p3 = appr_reduce(s,y,N);
+   else d = NULL;
-   if (DEBUGLEVEL>2)
+   z = appr_reduce(z, x);
-     { fprintferr(" sortie de idealappr0 p3 = "); outerr(p3); }
+   if (d) z = gdiv(z, d);
-   return gerepileupto(av,p3);
+   return z;
  }
+ GEN
+ idealapprfact(GEN nf, GEN x) {
+   gpmem_t av = avma;
+   return gerepileupto(av, idealapprfact_i(nf, x));
+ }
+ GEN
+ idealappr(GEN nf, GEN x) {
+   gpmem_t av = avma;
+   return gerepileupto(av, idealapprfact_i(nf, idealfactor(nf, x)));
+ }
+ GEN
+ idealappr0(GEN nf, GEN x, long fl) {
+   return fl? idealapprfact(nf, x): idealappr(nf, x);
+ }
  /* Given a prime ideal factorization x with possibly zero or negative exponents,
-  * and a vector y of elements of nf, gives a b such that
+  * and a vector w of elements of nf, gives a b such that
-  * v_p(b-y_p)>=v_p(x) for all prime ideals p in the ideal factorization
+  * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
   * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
-  * Certainly not the most efficient, but sure.
+  * Certainly not the most efficient, but sure. */
-  */
  GEN
- idealchinese(GEN nf, GEN x, GEN y)
+ idealchinese(GEN nf, GEN x, GEN w)
  {
-   long ty=typ(y),av=avma,i,j,k,l,N,r,r2;
+   gpmem_t av = avma;
-   GEN fact,fact2,list,ep,ep1,ep2,z,t,v,p1,p2,p3,p4,s,pr,den;
+   long ty = typ(w), i,N,r;
+   GEN L,e,z,t,y,v,s,den;
-   if (DEBUGLEVEL>4)
-   {
-     fprintferr(" entree dans idealchinese() :\n");
-     fprintferr(" x = "); outerr(x);
-     fprintferr(" y = "); outerr(y);
-   }
    nf=checknf(nf); N=degpol(nf[1]);
-   if (typ(x)!=t_MAT ||(lg(x)!=3))
+   if (typ(x) != t_MAT || lg(x) != 3)
      err(talker,"not a prime ideal factorization in idealchinese");
-   fact=x; list=(GEN)fact[1]; ep=(GEN)fact[2]; r=lg(list);
+   L = (GEN)x[1]; r = lg(L);
-   if (!is_vec_t(ty) || lg(y)!=r)
+   e = (GEN)x[2];
+   if (!is_vec_t(ty) || lg(w) != r)
      err(talker,"not a suitable vector of elements in idealchinese");
    if (r==1) return gscalcol_i(gun,N);
-   den=denom(y);
+   den = denom(w);
    if (!gcmp1(den))
    {
-     fact2=idealfactor(nf,den);
+     GEN fa = famat_reduce(famat_mul(x, idealfactor(nf,den)));
-     p1=(GEN)fact2[1]; r2=lg(p1);
+     L = (GEN)fa[1]; r = lg(L);
-     ep2=(GEN)fact2[2]; l=r+r2-1;
+     e = (GEN)fa[2];
-     z=cgetg(l,t_VEC); for (i=1; i<r; i++) z[i]=list[i];
-     ep1=cgetg(l,t_VEC); for (i=1; i<r; i++) ep1[i]=ep[i];
-     j=r-1;
-     for (i=1; i<r2; i++)
-     {
-       p3=(GEN)p1[i]; k=1;
-       while (k<r && !gegal((GEN)list[k],p3)) k++;
-       if (k==r) { j++; z[j]=(long)p3; ep1[j]=ep2[i]; }
-       else ep1[k]=ladd((GEN)ep1[k],(GEN)ep2[i]);
-     }
-     r=j+1; setlg(z,r); setlg(ep1,r); list=z; ep=ep1;
    }
    for (i=1; i<r; i++)
-     if (signe(ep[i])<0) ep[i] = zero;
+     if (signe(e[i]) < 0) e[i] = zero;
-   t=idmat(N);
+   y = factorbackprime(nf, L, e);
+   z = cgetg(r,t_VEC);
    for (i=1; i<r; i++)
+     z[i] = (long)idealdivpowprime(nf,y, (GEN)L[i], (GEN)e[i]);
+   v = idealaddmultoone(nf,z); s = NULL;
+   for (i=1; i<r; i++)
    {
-     pr=(GEN)list[i]; p4=(GEN)ep[i];
+     t = element_mul(nf, (GEN)v[i], (GEN)w[i]);
-     if (signe(p4))
+     s = s? gadd(s,t): t;
-     {
-       if (cmpis((GEN)pr[4],N))
-       {
-         p2=cgetg(3,t_MAT);
-         p2[1]=(long)gscalcol_i(powgi((GEN)pr[1],p4), N);
-         p2[2]=(long)element_pow(nf,(GEN)pr[2],p4);
-         t=idealmat_mul(nf,t,p2);
-       }
-       else t=gmul(powgi((GEN)pr[1],p4),t);
-     }
    }
-   z=cgetg(r,t_VEC);
+   return gerepileupto(av, appr_reduce(s,y));
+ }
+ static GEN
+ mat_ideal_two_elt2(GEN nf, GEN x, GEN a)
+ {
+   GEN L, e, fact = idealfactor(nf,a);
+   long i, v, r;
+   L = (GEN)fact[1]; r = lg(L);
+   e = (GEN)fact[2];
    for (i=1; i<r; i++)
    {
-     pr=(GEN)list[i]; p4=(GEN)ep[i];
+     v = idealval(nf,x,(GEN)L[i]);
-     if (cmpis((GEN)pr[4],N))
+     e[i] = lstoi(v);
-     {
-       p2=cgetg(3,t_MAT); p1=powgi((GEN)pr[1],p4);
-       p2[1]=(long)gscalcol_i(p1,N);
-       p2[2]=(long)element_pow(nf,(GEN)pr[5],p4);
-       z[i]=ldiv(idealmat_mul(nf,t,p2),p1);
-     }
-     else z[i]=ldiv(t,powgi((GEN)pr[1],p4));
    }
-   v=idealaddmultoone(nf,z);
+   return centermod(idealapprfact_i(nf,fact), gcoeff(x,1,1));
-   s=cgetg(N+1,t_COL); for (i=1; i<=N; i++) s[i]=zero;
-   for (i=1; i<r; i++)
-     s = gadd(s,element_mul(nf,(GEN)v[i],(GEN)y[i]));
-   p3 = appr_reduce(s,t,N);
-   if (DEBUGLEVEL>2)
-     { fprintferr(" sortie de idealchinese() : p3 = "); outerr(p3); }
-   return gerepileupto(av,p3);
  }
- GEN
- idealappr(GEN nf, GEN x) { return idealappr0(nf,x,0); }
- GEN
- idealapprfact(GEN nf, GEN x) { return idealappr0(nf,x,1); }
  /* Given an integral ideal x and a in x, gives a b such that
-  * x=aZ_K+bZ_K using a different algorithm than ideal_two_elt
+  * x = aZ_K + bZ_K using the approximation theorem */
-  */
  GEN
  ideal_two_elt2(GEN nf, GEN x, GEN a)
  {
-   long ta=typ(a), av=avma,tetpil,i,r;
+   gpmem_t av = avma;
-   GEN con,ep,b,list,fact;
+   GEN cx, b;
    nf = checknf(nf);
-   if (!is_extscalar_t(ta) && typ(a)!=t_COL)
+   if (typ(a) != t_COL) a = algtobasis(nf, a);
-     err(typeer,"ideal_two_elt2");
    x = idealhermite_aux(nf,x);
    if (gcmp0(x))
    {
      if (!gcmp0(a)) err(talker,"element not in ideal in ideal_two_elt2");
      avma=av; return gcopy(a);
    }
-   con = content(x);
+   x = Q_primitive_part(x, &cx);
-   if (gcmp1(con)) con = NULL; else { x = gdiv(x,con); a = gdiv(a,con); }
+   if (cx) a = gdiv(a, cx);
-   a = principalideal(nf,a);
+   if (!hnf_invimage(x, a))
-   if (!gcmp1(denom(gauss(x,a))))
      err(talker,"element does not belong to ideal in ideal_two_elt2");
-   fact=idealfactor(nf,a); list=(GEN)fact[1];
+   b = mat_ideal_two_elt2(nf, x, a);
-   r=lg(list); ep = (GEN)fact[2];
+   b = cx? gmul(b,cx): gcopy(b);
-   for (i=1; i<r; i++) ep[i]=lstoi(idealval(nf,x,(GEN)list[i]));
+   return gerepileupto(av, b);
-   b = centermod(idealappr0(nf,fact,1), gcoeff(x,1,1));
-   tetpil=avma; b = con? gmul(b,con): gcopy(b);
-   return gerepile(av,tetpil,b);
  }
- /* Given 2 integral ideals x and y in a number field nf gives a beta
+ /* Given 2 integral ideals x and y in nf, returns a beta in nf such that
-  * belonging to nf such that beta.x is an integral ideal coprime to y
+  * beta * x is an integral ideal coprime to y */
-  */
  GEN
  idealcoprime(GEN nf, GEN x, GEN y)
  {
-   long av=avma,tetpil,i,r;
+   gpmem_t av = avma;
-   GEN fact,list,p2,ep;
+   long i, r;
+   GEN list, ep, fact = idealfactor(nf,y);
-   if (DEBUGLEVEL>4)
+   list = (GEN)fact[1];
-   {
+   ep   = (GEN)fact[2]; r = lg(ep);
-     fprintferr(" entree dans idealcoprime() :\n");
+   for (i=1; i<r; i++) ep[i] = lstoi(-idealval(nf,x,(GEN)list[i]));
-     fprintferr(" x = "); outerr(x);
+   return gerepileupto(av, idealapprfact_i(nf,fact));
-     fprintferr(" y = "); outerr(y);
-   }
-   fact=idealfactor(nf,y); list=(GEN)fact[1];
-   r=lg(list); ep = (GEN)fact[2];
-   for (i=1; i<r; i++) ep[i]=lstoi(-idealval(nf,x,(GEN)list[i]));
-   tetpil=avma; p2=idealappr0(nf,fact,1);
-   if (DEBUGLEVEL>4)
-     { fprintferr(" sortie de idealcoprime() : p2 = "); outerr(p2); }
-   return gerepile(av,tetpil,p2);
  }
+ /*******************************************************************/
+ /*                                                                 */
+ /*                  LINEAR ALGEBRA OVER Z_K  (HNF,SNF)             */
+ /*                                                                 */
+ /*******************************************************************/
  /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
  long
  isinvector(GEN v, GEN x, long n)
-Line 2347  isinvector(GEN v, GEN x, long n)
+Line 2625  isinvector(GEN v, GEN x, long n)
 Line 2347  isinvector(GEN v, GEN x, long n)
 Line 2625  isinvector(GEN v, GEN x, long n)
  }
  GEN
- elt_mul_get_table(GEN nf, GEN x)
- {
-   long i,lx = lg(x);
-   GEN mul=cgetg(lx,t_MAT);
-   /* assume w_1 = 1 */
-   mul[1]=(long)x;
-   for (i=2; i<lx; i++) mul[i] = (long)element_mulid(nf,x,i);
-   return mul;
- }
- GEN
- elt_mul_table(GEN mul, GEN z)
- {
-   long av = avma, i, lx = lg(mul);
-   GEN p1 = gmul((GEN)z[1], (GEN)mul[1]);
-   for (i=2; i<lx; i++)
-     if (!gcmp0((GEN)z[i])) p1 = gadd(p1, gmul((GEN)z[i], (GEN)mul[i]));
-   return gerepileupto(av, p1);
- }
- GEN
  element_mulvec(GEN nf, GEN x, GEN v)
  {
-   long lv=lg(v),i;
+   long lv = lg(v), i;
    GEN y = cgetg(lv,t_COL);
    if (typ(x) == t_COL)
    {
-     GEN mul = elt_mul_get_table(nf,x);
+     GEN mul = eltmul_get_table(nf,x);
-     for (i=1; i<lv; i++)
+     for (i=1; i<lv; i++) y[i] = lmul(mul,(GEN)v[i]);
-       y[i] = (long)elt_mul_table(mul,(GEN)v[i]);
    }
    else
    { /* scalar */
-     for (i=1; i<lv; i++)
+     for (i=1; i<lv; i++) y[i] = lmul(x, (GEN)v[i]);
-       y[i] = lmul(x, (GEN)v[i]);
    }
    return y;
  }
-Line 2393  static GEN
+Line 2646  static GEN
 Line 2393  static GEN
 Line 2646  static GEN
  element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
  {
    long lv,j;
-   GEN y, mul = elt_mul_get_table(nf,x);
+   GEN y, mul = eltmul_get_table(nf,x);
    lv=min(lg(m),lim+1); y=cgetg(lv,t_VEC);
-   for (j=1; j<lv; j++)
+   for (j=1; j<lv; j++) y[j] = lmul(mul,gcoeff(m,i,j));
-     y[j] = (long)elt_mul_table(mul,gcoeff(m,i,j));
    return y;
  }
  /* Given an element x and an ideal in matrix form (not necessarily HNF),
-  * gives an r such that x-r is in ideal and r is small. No checks
+  * gives an r such that x-r is in ideal and r is small. No checks */
-  */
  GEN
  element_reduce(GEN nf, GEN x, GEN ideal)
  {
-   long tx=typ(x),av=avma,tetpil,N,i;
+   long tx = typ(x), N, i;
-   GEN p1,u;
+   gpmem_t av=avma;
+   GEN u,d;
    if (is_extscalar_t(tx))
-     x = algtobasis_intern(checknf(nf), x);
+     x = algtobasis_i(checknf(nf), x);
    N = lg(x);
    if (typ(ideal) != t_MAT || lg(ideal) != N) err(typeer,"element_reduce");
-   p1=cgetg(N+1,t_MAT);
+   u = ker( concatsp(ideal,x) ); /* do NOT use deplin. Much, much slower -- KB */
-   for (i=1; i<N; i++) p1[i]=ideal[i];
+   u = (GEN)u[1];
-   p1[N]=(long)x; u=(GEN)ker(p1)[1];
+   d = (GEN)u[N]; setlg(u,N);
-   p1=(GEN)u[N]; setlg(u,N);
+   for (i=1; i<N; i++) u[i] = (long)gdivround((GEN)u[i],d);
-   for (i=1; i<N; i++) u[i]=lround(gdiv((GEN)u[i],p1));
+   return gerepileupto(av, gadd(x, gmul(ideal,u)));
-   u=gmul(ideal,u); tetpil=avma;
-   return gerepile(av,tetpil,gadd(u,x));
  }
  /* A torsion-free module M over Z_K will be given by a row vector [A,I] with
-Line 2430  element_reduce(GEN nf, GEN x, GEN ideal)
+Line 2680  element_reduce(GEN nf, GEN x, GEN ideal)
 Line 2430  element_reduce(GEN nf, GEN x, GEN ideal)
 Line 2680  element_reduce(GEN nf, GEN x, GEN ideal)
   * that [A,I] is a pseudo-basis if k=n
   */
- /* Given a torsion-free module x as above outputs a pseudo-basis for x in
+ #define swap(x,y) { long _t=x; x=y; y=_t; }
-  * Hermite Normal Form
-  */
+ /* Given a torsion-free module x as above outputs a pseudo-basis for x in HNF */
  GEN
  nfhermite(GEN nf, GEN x)
  {
-   long av0 = avma, av,lim,i,j,def,k,m;
+   long i, j, def, k, m;
+   gpmem_t av0 = avma, av, lim;
    GEN p1,p2,y,A,I,J;
-   nf=checknf(nf);
+   nf = checknf(nf);
    if (typ(x)!=t_VEC || lg(x)!=3) err(talker,"not a module in nfhermite");
-   A=(GEN)x[1]; I=(GEN)x[2]; k=lg(A)-1;
+   A = (GEN)x[1];
+   I = (GEN)x[2]; k = lg(A)-1;
    if (typ(A)!=t_MAT) err(talker,"not a matrix in nfhermite");
    if (typ(I)!=t_VEC || lg(I)!=k+1)
      err(talker,"not a correct ideal list in nfhermite");
    if (!k) err(talker,"not a matrix of maximal rank in nfhermite");
-   m=lg(A[1])-1;
+   m = lg(A[1])-1;
-   if (k<m) err(talker,"not a matrix of maximal rank in nfhermite");
+   if (k < m) err(talker,"not a matrix of maximal rank in nfhermite");
    av = avma; lim = stack_lim(av, 1);
-   p1 = cgetg(k+1,t_MAT); for (j=1; j<=k; j++) p1[j]=A[j];
+   A = matalgtobasis(nf,A);
-   A = p1; I = dummycopy(I);
+   I = dummycopy(I);
-   J = cgetg(k+1,t_VEC);
+   J = zerovec(k); def = k+1;
-   for (j=1; j<=k; j++)
-   {
-     if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
-     J[j] = zero;
-   }
-   def = k+1;
    for (i=m; i>=1; i--)
    {
      GEN den,p4,p5,p6,u,v,newid, invnewid = NULL;
      def--; j=def; while (j>=1 && gcmp0(gcoeff(A,i,j))) j--;
      if (!j) err(talker,"not a matrix of maximal rank in nfhermite");
-     if (j==def) j--;
+     if (j==def) j--; else { swap(A[j], A[def]); swap(I[j], I[def]); }
-     else
+     p1 = gcoeff(A,i,def);
-     {
+     p2 = element_inv(nf, p1);
-       p1=(GEN)A[j]; A[j]=A[def]; A[def]=(long)p1;
+     A[def] = (long)element_mulvec(nf,p2,(GEN)A[def]);
-       p1=(GEN)I[j]; I[j]=I[def]; I[def]=(long)p1;
+     I[def] = (long)idealmul(nf,p1,(GEN)I[def]);
-     }
-     p1=gcoeff(A,i,def); p2=element_inv(nf,p1);
-     A[def]=(long)element_mulvec(nf,p2,(GEN)A[def]);
-     I[def]=(long)idealmul(nf,p1,(GEN)I[def]);
      for (  ; j; j--)
      {
-       p1=gcoeff(A,i,j);
+       p1 = gcoeff(A,i,j);
        if (gcmp0(p1)) continue;
-       p2=idealmul(nf,p1,(GEN)I[j]);
+       p2 = idealmul(nf,p1,(GEN)I[j]);
        newid = idealadd(nf,p2,(GEN)I[def]);
        invnewid = hnfideal_inv(nf,newid);
        p4 = idealmul(nf, p2,        invnewid);
        p5 = idealmul(nf,(GEN)I[def],invnewid);
        y = idealaddtoone(nf,p4,p5);
-       u=element_div(nf,(GEN)y[1],p1); v=(GEN)y[2];
+       u = (GEN)y[1]; u = element_div(nf, u, p1);
-       p6=gsub((GEN)A[j],element_mulvec(nf,p1,(GEN)A[def]));
+       v = (GEN)y[2];
-       A[def]=ladd(element_mulvec(nf,u,(GEN)A[j]),
+       p6 = gsub((GEN)A[j], element_mulvec(nf, p1,(GEN)A[def]));
-                   element_mulvec(nf,v,(GEN)A[def]));
+       A[def] = ladd(element_mulvec(nf, u, (GEN)A[j]),
-       A[j]=(long)p6;
+                     element_mulvec(nf, v, (GEN)A[def]));
-       I[j]=(long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
+       A[j] = (long)p6;
-       I[def]=(long)newid; den=denom((GEN)I[j]);
+       I[j] = (long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
+       I[def] = (long)newid; den = denom((GEN)I[j]);
        if (!gcmp1(den))
        {
-         I[j]=lmul(den,(GEN)I[j]);
+         I[j] = lmul(den,(GEN)I[j]);
-         A[j]=ldiv((GEN)A[j],den);
+         A[j] = ldiv((GEN)A[j],den);
        }
      }
      if (!invnewid) invnewid = hnfideal_inv(nf,(GEN)I[def]);
-     p1=(GEN)I[def]; J[def]=(long)invnewid;
+     J[def] = (long)invnewid;
+     p1 = (GEN)I[def];
      for (j=def+1; j<=k; j++)
      {
-       p2 = gsub(element_reduce(nf,gcoeff(A,i,j),idealmul(nf,p1,(GEN)J[j])),
+       GEN c = gcoeff(A,i,j);
-                 gcoeff(A,i,j));
+       p2 = gsub(element_reduce(nf,c, idealmul(nf,p1,(GEN)J[j])), c);
-       A[j] = ladd((GEN)A[j],element_mulvec(nf,p2,(GEN)A[def]));
+       A[j] = ladd((GEN)A[j], element_mulvec(nf, p2,(GEN)A[def]));
      }
      if (low_stack(lim, stack_lim(av1,1)))
      {
-       GEN *gptr[3];
+       GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
        if(DEBUGMEM>1) err(warnmem,"nfhermite, i = %ld", i);
-       gptr[0]=&A; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+       gerepilemany(av,gptr,3);
      }
    }
    y = cgetg(3,t_VEC);
-   p1 = cgetg(m+1,t_MAT); y[1] = (long)p1;
+   A += k-m; A[0] = evaltyp(t_MAT)|evallg(m+1); y[1] = (long)A;
-   p2 = cgetg(m+1,t_VEC); y[2] = (long)p2;
+   I += k-m; I[0] = evaltyp(t_VEC)|evallg(m+1); y[2] = (long)I;
-   for (j=1; j<=m; j++) p1[j] = lcopy((GEN)A[j+k-m]);
+   return gerepilecopy(av0, y);
-   for (j=1; j<=m; j++) p2[j] = lcopy((GEN)I[j+k-m]);
-   return gerepileupto(av0, y);
  }
+ static GEN
+ idealmulelt(GEN nf, GEN elt, GEN x)
+ {
+   long t = typ(elt);
+   GEN z;
+   if (t == t_POL || t == t_POLMOD) elt = algtobasis(nf,elt);
+   if (isnfscalar(elt)) elt = (GEN)elt[1];
+   z = element_mulvec(nf, elt, x);
+   settyp(z, t_MAT); return z;
+ }
+ static GEN
+ zero_nfbezout(GEN nf,GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
+ {
+   gpmem_t av, tetpil;
+   GEN pab,d;
+   d=idealmulelt(nf,b,B); *di=idealinv(nf,idealmat_to_hnf(nf,d));
+   av=avma; pab=idealmul(nf,A,B); tetpil=avma;
+   *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
+   *v=element_inv(nf,b);
+   *u=gzero; return d;
+ }
+ /* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
+  * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in
+  * B.di. Assume A, B non-zero, but a or b can be zero (not both)
+  */
+ static GEN
+ nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *u,GEN *v,GEN *w,GEN *di)
+ {
+   GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
+   gpmem_t av, tetpil;
+   if (gcmp0(a))
+   {
+     if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
+     return zero_nfbezout(nf,b,A,B,u,v,w,di);
+   }
+   if (gcmp0(b))
+     return zero_nfbezout(nf,a,B,A,v,u,w,di);
+   av = avma;
+   pa = idealmulelt(nf,a,A);
+   pb = idealmulelt(nf,b,B);
+   d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
+   pa1=idealmullll(nf,pa,dinv);
+   pb1=idealmullll(nf,pb,dinv);
+   uv=idealaddtoone(nf,pa1,pb1);
+   pab=idealmul(nf,A,B); tetpil=avma;
+   pu=element_div(nf,(GEN)uv[1],a);
+   pv=element_div(nf,(GEN)uv[2],b);
+   d=gcopy(d); dinv=gcopy(dinv);
+   pw=idealmul(nf,pab,dinv);
+   *u=pu; *v=pv; *w=pw; *di=dinv;
+   gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
+   gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
+   return d;
+ }
  /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
   * three components. I=[b_1,...,b_n] is a row vector of k fractional ideals
   * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
-Line 2536  nfhermite(GEN nf, GEN x)
+Line 2841  nfhermite(GEN nf, GEN x)
 Line 2536  nfhermite(GEN nf, GEN x)
 Line 2841  nfhermite(GEN nf, GEN x)
  GEN
  nfsmith(GEN nf, GEN x)
  {
-   long av,tetpil,i,j,k,l,lim,c,n,m,N;
+   long i, j, k, l, c, n, m, N;
+   gpmem_t av, tetpil, lim;
    GEN p1,p2,p3,p4,z,b,u,v,w,d,dinv,unnf,A,I,J;
    nf=checknf(nf); N=degpol(nf[1]);
-Line 2562  nfsmith(GEN nf, GEN x)
+Line 2868  nfsmith(GEN nf, GEN x)
 Line 2562  nfsmith(GEN nf, GEN x)
 Line 2868  nfsmith(GEN nf, GEN x)
    {
      do
      {
-       c=0;
+       c = 0;
        for (j=i-1; j>=1; j--)
        {
-         p1=gcoeff(A,i,j);
+         p1=gcoeff(A,i,j); if (gcmp0(p1)) continue;
-         if (!gcmp0(p1))
-         {
+         p2 = gcoeff(A,i,i);
-           p2=gcoeff(A,i,i);
+         d = nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
-           d=nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
+         if (gcmp0(u)) b = element_mulvec(nf,v,(GEN)A[j]);
-           if (!gcmp0(u))
+         else
-           {
+         {
-             if (!gcmp0(v))
+           if (gcmp0(v)) b = element_mulvec(nf,u,(GEN)A[i]);
-               b=gadd(element_mulvec(nf,u,(GEN)A[i]),
+           else
-                      element_mulvec(nf,v,(GEN)A[j]));
+             b = gadd(element_mulvec(nf,u,(GEN)A[i]),
-             else b=element_mulvec(nf,u,(GEN)A[i]);
+                      element_mulvec(nf,v,(GEN)A[j]));
-           }
+         }
-           else b=element_mulvec(nf,v,(GEN)A[j]);
+         A[j] = lsub(element_mulvec(nf,p2,(GEN)A[j]),
-           A[j]=lsub(element_mulvec(nf,p2,(GEN)A[j]),
+                     element_mulvec(nf,p1,(GEN)A[i]));
-                     element_mulvec(nf,p1,(GEN)A[i]));
+         A[i] = (long)b; J[j] = (long)w; J[i] = (long)d;
-           A[i]=(long)b; J[j]=(long)w; J[i]=(long)d;
-         }
        }
        for (j=i-1; j>=1; j--)
        {
-         p1=gcoeff(A,j,i);
+         p1 = gcoeff(A,j,i); if (gcmp0(p1)) continue;
-         if (!gcmp0(p1))
-         {
+         p2 = gcoeff(A,i,i);
-           p2=gcoeff(A,i,i);
+         d = nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
-           d=nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
+         if (gcmp0(u)) b = element_mulvecrow(nf,v,A,j,i);
-           if (gcmp0(u))
+         else
-             b=element_mulvecrow(nf,v,A,j,i);
+         {
-           else
+           if (gcmp0(v))
-           {
+             b = element_mulvecrow(nf,u,A,i,i);
-             if (gcmp0(v))
+           else
-               b=element_mulvecrow(nf,u,A,i,i);
+             b = gadd(element_mulvecrow(nf,u,A,i,i),
-             else
+                      element_mulvecrow(nf,v,A,j,i));
-               b=gadd(element_mulvecrow(nf,u,A,i,i),
+         }
-                      element_mulvecrow(nf,v,A,j,i));
+         p3 = gsub(element_mulvecrow(nf,p2,A,j,i),
-           }
+                   element_mulvecrow(nf,p1,A,i,i));
-           p3=gsub(element_mulvecrow(nf,p2,A,j,i),
+         for (k=1; k<=i; k++) { coeff(A,j,k) = p3[k]; coeff(A,i,k) = b[k]; }
-                   element_mulvecrow(nf,p1,A,i,i));
+         I[j] = (long)w; I[i] = (long)d; c = 1;
-           for (k=1; k<=i; k++) { coeff(A,j,k)=p3[k]; coeff(A,i,k)=b[k]; }
-           I[j]=(long)w; I[i]=(long)d; c++;
-         }
        }
-       if (!c)
+       if (c) continue;
-       {
-         b=gcoeff(A,i,i); if (gcmp0(b)) break;
-         b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
+       b=gcoeff(A,i,i); if (gcmp0(b)) break;
-         for (k=1; k<i; k++)
-           for (l=1; l<i; l++)
-           {
-             p3 = gcoeff(A,k,l);
-             if (!gcmp0(p3))
-             {
-               p4 = idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]));
-               if (!gegal(idealadd(nf,b,p4), b))
-               {
-                 b=idealdiv(nf,(GEN)I[k],(GEN)I[i]);
-                 p4=gauss(idealdiv(nf,(GEN)J[i],idealmul(nf,p3,(GEN)J[l])),b);
-                 l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
-                 if (l>N) err(talker,"bug2 in nfsmith");
-                 p3=element_mulvecrow(nf,(GEN)b[l],A,k,i);
-                 for (l=1; l<=i; l++)
-                   coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
-                 k = l = i; c = 1;
+       b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
-               }
+       for (k=1; k<i; k++)
-             }
+         for (l=1; l<i; l++)
-           }
+         {
-       }
+           p3 = gcoeff(A,k,l);
+           if (gcmp0(p3)) continue;
+           p4 = idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]));
+           if (hnfdivide(b, p4)) continue;
+           b = idealdiv(nf,(GEN)I[k],(GEN)I[i]);
+           p1 = idealdiv(nf,(GEN)J[i], idealmul(nf,p3,(GEN)J[l]));
+           p4 = gauss(p4, b);
+           l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
+           if (l>N) err(talker,"bug2 in nfsmith");
+           p3 = element_mulvecrow(nf,(GEN)b[l],A,k,i);
+           for (l=1; l<=i; l++)
+             coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
+           k = l = i; c = 1;
+         }
        if (low_stack(lim, stack_lim(av,1)))
        {
-         GEN *gptr[3];
+         GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
          if(DEBUGMEM>1) err(warnmem,"nfsmith");
-         gptr[0]=&A; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+         gerepilemany(av,gptr,3);
        }
      }
      while (c);
-Line 2652  nfsmith(GEN nf, GEN x)
+Line 2951  nfsmith(GEN nf, GEN x)
 Line 2652  nfsmith(GEN nf, GEN x)
 Line 2951  nfsmith(GEN nf, GEN x)
    return gerepile(av,tetpil,z);
  }
- /*******************************************************************/
- /*                                                                 */
- /*          ALGEBRE LINEAIRE DANS LES CORPS DE NOMBRES             */
- /*                                                                 */
- /*******************************************************************/
  #define trivlift(x) ((typ(x)==t_POLMOD)? (GEN)x[2]: lift_intern(x))
  GEN
- element_mulmodpr2(GEN nf, GEN x, GEN y, GEN prhall)
+ element_mulmodpr(GEN nf, GEN x, GEN y, GEN modpr)
  {
-   long av=avma;
+   gpmem_t av=avma;
    GEN p1;
-   nf=checknf(nf); checkprhall(prhall);
+   nf=checknf(nf); checkmodpr(modpr);
    p1 = element_mul(nf,x,y);
-   return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+   return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
  }
- /* On ne peut PAS definir ca comme les autres par
-  * #define element_divmodpr() nfreducemodpr(element_div())
-  * car le element_div ne marche pas en general
-  */
  GEN
- element_divmodpr(GEN nf, GEN x, GEN y, GEN prhall)
+ element_divmodpr(GEN nf, GEN x, GEN y, GEN modpr)
  {
-   long av=avma;
+   gpmem_t av=avma;
    GEN p1;
-   nf=checknf(nf); checkprhall(prhall);
+   nf=checknf(nf); checkmodpr(modpr);
    p1=lift_intern(gdiv(gmodulcp(gmul((GEN)nf[7],trivlift(x)), (GEN)nf[1]),
                        gmodulcp(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1])));
-   p1=algtobasis_intern(nf,p1);
+   p1=algtobasis_i(nf,p1);
-   return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+   return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
  }
  GEN
- element_invmodpr(GEN nf, GEN y, GEN prhall)
+ element_invmodpr(GEN nf, GEN y, GEN modpr)
  {
-   long av=avma;
+   gpmem_t av=avma;
    GEN p1;
-   p1=ginvmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
+   p1=QX_invmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
-   p1=algtobasis_intern(nf,p1);
+   p1=algtobasis_i(nf,p1);
-   return gerepileupto(av,nfreducemodpr(nf,p1,prhall));
+   return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
  }
  GEN
- element_powmodpr(GEN nf,GEN x,GEN k,GEN prhall)
+ element_powmodpr(GEN nf,GEN x,GEN k,GEN pr)
  {
-   long av=avma,N,s;
+   gpmem_t av=avma;
-   GEN y,z;
+   GEN z,T,p,modpr;
-   nf=checknf(nf); checkprhall(prhall);
+   nf = checknf(nf);
-   N=degpol(nf[1]);
+   modpr = nf_to_ff_init(nf,&pr,&T,&p);
-   s=signe(k); k=(s>=0)?k:negi(k);
+   z = nf_to_ff(nf,x,modpr);
-   z=x; y = gscalcol_i(gun,N);
+   z = FpXQ_pow(z,k,T,p);
-   for(;;)
+   z = ff_to_nf(z,modpr);
-   {
+   if (typ(z) != t_COL) z = algtobasis(nf, z);
-     if (mpodd(k)) y=element_mulmodpr(nf,z,y,prhall);
+   return gerepilecopy(av, z);
-     k=shifti(k,-1);
-     if (signe(k)) z=element_sqrmodpr(nf,z,prhall);
-     else
-     {
-       cgiv(k); if (s<0) y = element_invmodpr(nf,y,prhall);
-       return gerepileupto(av,y);
-     }
-   }
  }
- /* x est une matrice dont les coefficients sont des vecteurs dans la base
-  * d'entiers modulo un ideal premier prhall, sous forme reduite modulo prhall.
-  */
  GEN
- nfkermodpr(GEN nf, GEN x, GEN prhall)
+ nfkermodpr(GEN nf, GEN x, GEN pr)
  {
-   ulong av0,av,av1,lim;
+   gpmem_t av = avma;
-   long i,j,k,r,t,n,m,N;
+   GEN T,p,modpr;
-   GEN c,d,y,unnf,munnf,zeromodp,zeronf,p,pp,prh;
-   nf=checknf(nf); checkprhall(prhall);
+   nf = checknf(nf);
    if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
-   n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
+   modpr = nf_to_ff_init(nf, &pr,&T,&p);
-   prh=(GEN)prhall[1]; av0=avma;
+   x = modprM(lift(x), nf, modpr);
-   N=degpol(nf[1]); pp=gcoeff(prh,1,1);
+   return gerepilecopy(av, modprM_lift(FqM_ker(x,T,p), modpr));
-   zeromodp=gmodulsg(0,pp);
-   unnf=cgetg(N+1,t_COL); unnf[1]=(long)gmodulsg(1,pp);
-   zeronf=cgetg(N+1,t_COL); zeronf[1]=(long)zeromodp;
-   av=avma; munnf=cgetg(N+1,t_COL); munnf[1]=(long)gmodulsg(-1,pp);
-   for (i=2; i<=N; i++)
-     zeronf[i] = munnf[i] = unnf[i] = (long)zeromodp;
-   m=lg(x[1])-1; x=dummycopy(x); r=0;
-   c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
-   d=new_chunk(n+1); av1=avma; lim=stack_lim(av1,1);
-   for (k=1; k<=n; k++)
-   {
-     j=1;
-     while (j<=m && (c[j] || gcmp0(gcoeff(x,j,k)))) j++;
-     if (j > m) { r++; d[k]=0; continue; }
-       p=element_divmodpr(nf,munnf,gcoeff(x,j,k),prhall);
-       c[j]=k; d[k]=j; coeff(x,j,k)=(long)munnf;
-       for (i=k+1; i<=n; i++)
-         coeff(x,j,i)=(long)element_mulmodpr(nf,p,gcoeff(x,j,i),prhall);
-       for (t=1; t<=m; t++)
-         if (t!=j)
-         {
-         p = gcoeff(x,t,k); if (gcmp0(p)) continue;
-         coeff(x,t,k) = (long)zeronf;
-           for (i=k+1; i<=n; i++)
-             coeff(x,t,i)=ladd(gcoeff(x,t,i),
-                               element_mulmodpr(nf,p,gcoeff(x,j,i),prhall));
-           if (low_stack(lim, stack_lim(av1,1)))
-           {
-             if (DEBUGMEM>1) err(warnmem,"nfkermodpr, k = %ld / %ld",k,n);
-             x=gerepilecopy(av1,x);
-           }
-         }
-   }
-   if (!r) { avma=av0; return cgetg(1,t_MAT); }
-   av1=avma; y=cgetg(r+1,t_MAT);
-   for (j=k=1; j<=r; j++,k++)
-   {
-     p=cgetg(n+1,t_COL); y[j]=(long)p; while (d[k]) k++;
-     for (i=1; i<k; i++) p[i]=d[i]? lcopy(gcoeff(x,d[i],k)): (long)zeronf;
-     p[k]=(long)unnf; for (i=k+1; i<=n; i++) p[i]=(long)zeronf;
-   }
-   return gerepile(av,av1,y);
  }
  /* a.x=b ou b est un vecteur */
  GEN
- nfsolvemodpr(GEN nf, GEN a, GEN b, GEN prhall)
+ nfsolvemodpr(GEN nf, GEN a, GEN b, GEN pr)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long nbli,nbco,i,j,k;
+   GEN T,p,modpr;
-   GEN aa,x,p,m,u;
-   nf=checknf(nf); checkprhall(prhall);
+   nf = checknf(nf);
    if (typ(a)!=t_MAT) err(typeer,"nfsolvemodpr");
-   nbco=lg(a)-1; nbli=lg(a[1])-1;
+   modpr = nf_to_ff_init(nf, &pr,&T,&p);
-   if (typ(b)!=t_COL) err(typeer,"nfsolvemodpr");
+   a = modprM(lift(a), nf, modpr);
-   if (lg(b)!=nbco+1) err(mattype1,"nfsolvemodpr");
+   b = modprM(lift(b), nf, modpr);
-   x=cgetg(nbli+1,t_COL);
+   return gerepilecopy(av, modprM_lift(FqM_gauss(a,b,T,p), modpr));
-   for (j=1; j<=nbco; j++) x[j]=b[j];
-   aa=cgetg(nbco+1,t_MAT);
-   for (j=1; j<=nbco; j++)
-   {
-     aa[j]=lgetg(nbli+1,t_COL);
-     for (i=1; i<=nbli; i++) coeff(aa,i,j)=coeff(a,i,j);
-   }
-   for (i=1; i<nbli; i++)
-   {
-     p=gcoeff(aa,i,i); k=i;
-     if (gcmp0(p))
-     {
-       k=i+1; while (k<=nbli && gcmp0(gcoeff(aa,k,i))) k++;
-       if (k>nbco) err(matinv1);
-       for (j=i; j<=nbco; j++)
-       {
-         u=gcoeff(aa,i,j); coeff(aa,i,j)=coeff(aa,k,j);
-         coeff(aa,k,j)=(long)u;
-       }
-       u=(GEN)x[i]; x[i]=x[k]; x[k]=(long)u;
-       p=gcoeff(aa,i,i);
-     }
-     for (k=i+1; k<=nbli; k++)
-     {
-       m=gcoeff(aa,k,i);
-       if (!gcmp0(m))
-       {
-         m=element_divmodpr(nf,m,p,prhall);
-         for (j=i+1; j<=nbco; j++)
-           coeff(aa,k,j)=lsub(gcoeff(aa,k,j),
-                              element_mulmodpr(nf,m,gcoeff(aa,i,j),prhall));
-         x[k]=lsub((GEN)x[k],element_mulmodpr(nf,m,(GEN)x[i],prhall));
-       }
-     }
-   }
-   /* Resolution systeme triangularise */
-   p=gcoeff(aa,nbli,nbco); if (gcmp0(p)) err(matinv1);
-   x[nbli]=(long)element_divmodpr(nf,(GEN)x[nbli],p,prhall);
-   for (i=nbli-1; i>0; i--)
-   {
-     m=(GEN)x[i];
-     for (j=i+1; j<=nbco; j++)
-       m=gsub(m,element_mulmodpr(nf,gcoeff(aa,i,j),(GEN)x[j],prhall));
-     x[i]=(long)element_divmodpr(nf,m,gcoeff(aa,i,i),prhall);
-   }
-   return gerepilecopy(av,x);
  }
+ #if 0
  GEN
- nfsuppl(GEN nf, GEN x, long n, GEN prhall)
+ nfsuppl(GEN nf, GEN x, GEN pr)
  {
-   long av=avma,av2,k,s,t,N,lx=lg(x);
+   gpmem_t av = avma;
-   GEN y,p1,p2,p,unmodp,zeromodp,unnf,zeronf,prh;
+   GEN T,p,modpr;
-   stackzone *zone;
-   k=lx-1; if (k>n) err(suppler2);
+   nf = checknf(nf);
-   if (k && lg(x[1])!=n+1) err(talker,"incorrect dimension in nfsupl");
+   if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
-   N=degpol(nf[1]); prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
+   modpr = nf_to_ff_init(nf, &pr,&T,&p);
+   x = modprM(lift(x), nf, modpr);
-   zone  = switch_stack(NULL, 2*(3 + 2*lg(p) + N+1) + (n+3)*(n+1));
+   return gerepilecopy(av, modprM_lift(FqM_suppl(x,T,p), modpr));
-   switch_stack(zone,1);
-   unmodp=gmodulsg(1,p); zeromodp=gmodulsg(0,p);
-   unnf=gscalcol_proto(unmodp,zeromodp,N);
-   zeronf=gscalcol_proto(zeromodp,zeromodp,N);
-   y = idmat_intern(n,unnf,zeronf);
-   switch_stack(zone,0); av2=avma;
-   for (s=1; s<=k; s++)
-   {
-     p1=nfsolvemodpr(nf,y,(GEN)x[s],prhall); t=s;
-     while (t<=n && gcmp0((GEN)p1[t])) t++;
-     avma=av2; if (t>n) err(suppler2);
-     p2=(GEN)y[s]; y[s]=x[s]; if (s!=t) y[t]=(long)p2;
-   }
-   avma=av; y=gcopy(y);
-   free(zone); return y;
  }
+ #endif
- /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
-    t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
- GEN
- nfidealdet1(GEN nf, GEN a, GEN b)
- {
-   long av=avma;
-   GEN x,p1,res,da,db;
-   a = idealinv(nf,a);
-   da = denom(a); if (!gcmp1(da)) a = gmul(da,a);
-   db = denom(b); if (!gcmp1(db)) b = gmul(db,b);
-   x = idealcoprime(nf,a,b);
-   p1 = idealaddtoone(nf, idealmul(nf,x,a), b);
-   res = cgetg(5,t_VEC);
-   res[1] = lmul(x,da);
-   res[2] = ldiv((GEN)p1[2],db);
-   res[3] = lnegi(db);
-   res[4] = (long) element_div(nf,(GEN)p1[1],(GEN)res[1]);
-   return gerepileupto(av,res);
- }
  /* Given a pseudo-basis pseudo, outputs a multiple of its ideal determinant */
  GEN
  nfdetint(GEN nf,GEN pseudo)
  {
    GEN pass,c,v,det1,piv,pivprec,vi,p1,x,I,unnf,zeronf,id,idprod;
-   long i,j,k,rg,n,n1,m,m1,av=avma,av1,tetpil,lim,cm=0,N;
+   long i, j, k, rg, n, n1, m, m1, cm=0, N;
+   gpmem_t av=avma, av1, tetpil, lim;
    nf=checknf(nf); N=degpol(nf[1]);
    if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
-Line 2999  nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
+Line 3146  nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
 Line 2999  nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
 Line 3146  nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
      x[i]=(long)element_reduce(nf,(GEN)x[i],detmat);
  }
- static GEN
- zero_nfbezout(GEN nf,GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
- {
-   long av, tetpil;
-   GEN pab,d;
-   d=idealmulelt(nf,b,B); *di=idealinv(nf,idealmat_to_hnf(nf,d));
-   av=avma; pab=idealmul(nf,A,B); tetpil=avma;
-   *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
-   *v=element_inv(nf,b);
-   *u=gzero; return d;
- }
- /* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
-  * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in
-  * B.di. Assume A, B non-zero, but a or b can be zero (not both)
-  */
- static GEN
- nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *u,GEN *v,GEN *w,GEN *di)
- {
-   GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
-   long av,tetpil;
-   if (gcmp0(a))
-   {
-     if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
-     return zero_nfbezout(nf,b,A,B,u,v,w,di);
-   }
-   if (gcmp0(b))
-     return zero_nfbezout(nf,a,B,A,v,u,w,di);
-   av = avma;
-   pa=idealmulelt(nf,a,A);
-   pb=idealmulelt(nf,b,B);
-   d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
-   pa1=idealmullll(nf,pa,dinv);
-   pb1=idealmullll(nf,pb,dinv);
-   uv=idealaddtoone(nf,pa1,pb1);
-   pab=idealmul(nf,A,B); tetpil=avma;
-   pu=element_div(nf,(GEN)uv[1],a);
-   pv=element_div(nf,(GEN)uv[2],b);
-   d=gcopy(d); dinv=gcopy(dinv);
-   pw=idealmul(nf,pab,dinv);
-   *u=pu; *v=pv; *w=pw; *di=dinv;
-   gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
-   gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
-   return d;
- }
  /* A usage interne. Pas de verifs ni gestion de pile */
  GEN
  idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
-Line 3062  idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
+Line 3157  idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
 Line 3062  idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
 Line 3157  idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
    return den? gdiv(z,den): z;
  }
- /* A usage interne. Pas de verifs ni gestion de pile */
  GEN
- idealmulelt(GEN nf, GEN elt, GEN x)
- {
-   long t = typ(elt);
-   GEN z;
-   if (t == t_POL || t == t_POLMOD) elt = algtobasis(nf,elt);
-   if (isnfscalar(elt)) elt = (GEN)elt[1];
-   z = element_mulvec(nf, elt, x);
-   settyp(z, t_MAT); return z;
- }
- GEN
  nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
  {
-   long av0=avma,li,co,av,tetpil,i,j,jm1,def,ldef,lim,N;
+   long li, co, i, j, jm1, def, ldef, N;
-   GEN b,q,w,p1,p2,d,u,v,den,x,I,J,dinv,unnf,wh;
+   gpmem_t av0=avma, av, lim;
+   GEN b,q,w,p1,d,u,v,den,x,I,J,dinv,wh,unnf;
    nf=checknf(nf); N=degpol(nf[1]);
    if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
      err(talker,"not a module in nfhermitemod");
-   x=(GEN)pseudo[1]; I=(GEN)pseudo[2];
+   x = (GEN)pseudo[1];
-   if (typ(x)!=t_MAT) err(talker,"not a matrix in nfhermitemod");
+   I = (GEN)pseudo[2];
-   co=lg(x);
+   if (typ(x) != t_MAT) err(talker,"not a matrix in nfhermitemod");
-   if (typ(I)!=t_VEC || lg(I)!=co)
+   co = lg(x);
+   if (typ(I) != t_VEC || lg(I) != co)
      err(talker,"not a correct ideal list in nfhermitemod");
    if (co==1) return cgetg(1,t_MAT);
-   li=lg(x[1]); x=dummycopy(x); I=dummycopy(I);
+   li= lg(x[1]);
-   unnf=gscalcol_i(gun,N);
+   x = matalgtobasis(nf, x);
-   for (j=1; j<co; j++)
+   I = dummycopy(I);
-     if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
+   unnf = gscalcol(gun,N);
-   den=denom(detmat); if (!gcmp1(den)) detmat=gmul(den,detmat);
+   den = denom(detmat); if (!gcmp1(den)) detmat = gmul(den,detmat);
-   detmat=gmul(detmat,lllintpartial(detmat));
+   detmat = lllint_ip(detmat,100);
-   av=avma; lim=stack_lim(av,1);
+   av = avma; lim = stack_lim(av,1);
-   def=co; ldef=(li>co)?li-co+1:1;
+   def = co; ldef = (li>co)? li-co+1: 1;
    for (i=li-1; i>=ldef; i--)
    {
      def--; j=def-1; while (j && gcmp0(gcoeff(x,i,j))) j--;
-Line 3124  nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
+Line 3209  nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
 Line 3124  nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
 Line 3209  nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-       GEN *gptr[2];
+       GEN *gptr[2]; gptr[0]=&x; gptr[1]=&I;
        if(DEBUGMEM>1) err(warnmem,"[1]: nfhermitemod");
-       gptr[0]=&x; gptr[1]=&I; gerepilemany(av,gptr,2);
+       gerepilemany(av,gptr,2);
      }
    }
-   b=detmat; wh=cgetg(li,t_MAT); def--;
+   b = detmat; wh = cgetg(li,t_MAT); def--;
    for (i=li-1; i>=1; i--)
    {
      d = nfbezout(nf,gcoeff(x,i,i+def),unnf,(GEN)I[i+def],b,&u,&v,&w,&dinv);
      p1 = element_mulvec(nf,u,(GEN)x[i+def]);
      nfcleanmod(nf,p1,i,idealmullll(nf,b,dinv));
-     wh[i]=(long)p1; coeff(wh,i,i)=(long)unnf; I[i+def]=(long)d;
+     wh[i] = (long)p1; coeff(wh,i,i) = (long)unnf;
-     if (i>1) b=idealmul(nf,b,dinv);
+     I[i+def] = (long)d;
+     if (i>1) b = idealmul(nf,b,dinv);
    }
    J=cgetg(li,t_VEC); J[1]=zero;
-   for (j=2; j<li; j++) J[j]=(long)idealinv(nf,(GEN)I[j+def]);
+   for (j=2; j<li; j++) J[j] = (long)idealinv(nf,(GEN)I[j+def]);
    for (i=li-2; i>=1; i--)
    {
      for (j=i+1; j<li; j++)
      {
-       q=idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
+       q = idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
-       p1=gsub(element_reduce(nf,gcoeff(wh,i,j),q),gcoeff(wh,i,j));
+       p1 = gsub(element_reduce(nf,gcoeff(wh,i,j),q), gcoeff(wh,i,j));
-       wh[j]=(long)gadd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
+       wh[j] = ladd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-       GEN *gptr[3];
+       GEN *gptr[3]; gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J;
        if(DEBUGMEM>1) err(warnmem,"[2]: nfhermitemod");
-       gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J; gerepilemany(av,gptr,3);
+       gerepilemany(av,gptr,3);
      }
    }
-   tetpil=avma; p1=cgetg(3,t_VEC); p1[1]=lcopy(wh);
+   I += def; I[0] = evaltyp(t_VEC)|evallg(li);
-   p2=cgetg(li,t_VEC); p1[2]=(long)p2;
+   p1=cgetg(3,t_VEC);
-   for (j=1; j<li; j++) p2[j]=lcopy((GEN)I[j+def]);
+   p1[1] = (long)wh;
-   return gerepile(av0,tetpil,p1);
+   p1[2] = (long)I; return gerepilecopy(av0, p1);
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>