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

Return to base2.c CVS log

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

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

-version 1.1, 2001/10/02 11:17:02
+version 1.2, 2002/09/11 07:26:49
 Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  /*                                                                 */
  /*******************************************************************/
  #include "pari.h"
+ #include "parinf.h"
- extern GEN caractducos(GEN p, GEN x, int v);
+ #define RXQX_rem(x,y,T) RXQX_divrem((x),(y),(T),ONLY_REM)
- extern GEN element_muli(GEN nf, GEN x, GEN y);
+ #define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM)
- extern GEN element_mulid(GEN nf, GEN x, long i);
+ extern GEN addone_aux2(GEN x, GEN y);
- extern GEN eleval(GEN f,GEN h,GEN a);
+ extern GEN addshiftw(GEN x, GEN y, long d);
- extern GEN ideal_better_basis(GEN nf, GEN x, GEN M);
+ extern GEN gmul_mat_smallvec(GEN x, GEN y);
- extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t, long v);
+ extern GEN hnf_invimage(GEN A, GEN b);
+ extern GEN norm_by_embed(long r1, GEN x);
+ extern GEN ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound);
+ extern GEN polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N);
+ extern GEN FqX_gcd(GEN P, GEN Q, GEN T, GEN p);
+ extern GEN FqX_factor(GEN x, GEN T, GEN p);
+ extern GEN DDF2(GEN x, long hint);
+ extern GEN eltabstorel(GEN x, GEN T, GEN pol, GEN k);
+ extern GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
+ extern GEN FpVQX_red(GEN z, GEN T, GEN p);
+ extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
+ extern GEN RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr);
+ extern GEN RXQX_mul(GEN x, GEN y, GEN T);
+ extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
+ extern GEN _ei(long n, long i);
+ extern GEN col_to_ff(GEN x, long v);
+ extern GEN element_mulidid(GEN nf, long i, long j);
+ extern GEN eltmulid_get_table(GEN nf, long i);
+ extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
  extern GEN mat_to_vecpol(GEN x, long v);
- extern GEN nfidealdet1(GEN nf, GEN a, GEN b);
+ extern GEN merge_factor_i(GEN f, GEN g);
- extern GEN nfsuppl(GEN nf, GEN x, long n, GEN prhall);
+ extern GEN mulmat_pol(GEN A, GEN x);
+ extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
+ extern GEN pidealprimeinv(GEN nf, GEN x);
  extern GEN pol_to_monic(GEN pol, GEN *lead);
- extern GEN pol_to_vec(GEN x, long N);
  extern GEN quicktrace(GEN x, GEN sym);
- extern GEN respm(GEN f1,GEN f2,GEN pm);
+ extern GEN sqr_by_tab(GEN tab, GEN x);
+ extern GEN to_polmod(GEN x, GEN mod);
+ extern GEN unnf_minus_x(GEN x);
+ extern int isrational(GEN x);
+ extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t);
+ /* FIXME: backward compatibility. Should use the proper nf_* equivalents */
+ #define compat_PARTIAL 1
+ #define compat_ROUND2  2
  static void
- allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1, GEN *ptw2)
+ allbase_check_args(GEN f, long flag, GEN *dx, GEN *ptw)
  {
-   GEN w;
+   GEN w = *ptw;
+   if (DEBUGLEVEL) (void)timer2();
    if (typ(f)!=t_POL) err(notpoler,"allbase");
    if (degpol(f) <= 0) err(constpoler,"allbase");
-   if (DEBUGLEVEL) timer2();
-   switch(code)
+   *dx = w? factorback(w, NULL): ZX_disc(f);
-   {
+   if (!signe(*dx)) err(talker,"reducible polynomial in allbase");
-     case 0: case 1:
+   if (!w) *ptw = auxdecomp(absi(*dx), (flag & nf_PARTIALFACT)? 0: 1);
-       *y = ZX_disc(f);
-       if (!signe(*y)) err(talker,"reducible polynomial in allbase");
-       w = auxdecomp(absi(*y),1-code);
-       break;
-     default:
-       w = (GEN)code;
-       *y = factorback(w, NULL);
-   }
    if (DEBUGLEVEL) msgtimer("disc. factorisation");
-   *ptw1 = (GEN)w[1];
-   *ptw2 = (GEN)w[2];
  }
  /*******************************************************************/
-Line 62  allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1
+Line 82  allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1
 Line 62  allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1
 Line 82  allbase_check_args(GEN f, long code, GEN *y, GEN *ptw1
  /*                            ROUND 2                              */
  /*                                                                 */
  /*******************************************************************/
- /*  Normalized quotient and remainder ( -1/2 |y| < r = x-q*y <= 1/2 |y| ) */
- static GEN
- rquot(GEN x, GEN y)
- {
-   long av=avma,av1;
-   GEN u,v,w,p;
-   u=absi(y); v=shifti(x,1); w=shifti(y,1);
-   if (cmpii(u,v)>0) p=subii(v,u);
-   else p=addsi(-1,addii(u,v));
-   av1=avma; return gerepile(av,av1,divii(p,w));
- }
- /* space needed lx + 2*ly */
- static GEN
- rrmdr(GEN x, GEN y)
- {
-   long av=avma,tetpil,k;
-   GEN r,ys2;
-   if (!signe(x)) return gzero;
-   r = resii(x,y); tetpil = avma;
-   ys2 = shifti(y,-1);
-   k = absi_cmp(r, ys2);
-   if (k>0 || (k==0 && signe(r)>0))
-   {
-     avma = tetpil;
-     if (signe(y) == signe(r)) r = subii(r,y); else r = addii(r,y);
-     return gerepile(av,tetpil,r);
-   }
-   avma = tetpil; return r;
- }
  /* companion matrix of unitary polynomial x */
  static GEN
  companion(GEN x) /* cf assmat */
-Line 117  companion(GEN x) /* cf assmat */
+Line 104  companion(GEN x) /* cf assmat */
 Line 117  companion(GEN x) /* cf assmat */
 Line 104  companion(GEN x) /* cf assmat */
  static GEN
  mulmati(GEN x, GEN y)
  {
-   long n = lg(x),i,j,k,av;
+   gpmem_t av;
+   long n = lg(x),i,j,k;
    GEN z = cgetg(n,t_MAT),p1,p2;
    for (j=1; j<n; j++)
-Line 138  mulmati(GEN x, GEN y)
+Line 126  mulmati(GEN x, GEN y)
 Line 138  mulmati(GEN x, GEN y)
 Line 126  mulmati(GEN x, GEN y)
  }
  static GEN
- powmati(GEN x, long m)
+ _mulmati(void *data /*ignored*/, GEN x, GEN y) {
- {
+   (void)data; return mulmati(x,y);
-   long av=avma,j;
+ }
-   GEN y = x;
+ static GEN
+ _sqrmati(void *data /*ignored*/, GEN x) {
+   (void)data; return mulmati(x,x);
+ }
-   j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
+ static GEN
-   for (; j; m<<=1,j--)
+ powmati(GEN x, GEN n)
-   {
+ {
-     y=mulmati(y,y);
+   gpmem_t av = avma;
-     if (m<0) y=mulmati(y,x);
+   GEN y = leftright_pow(x, n, NULL, &_sqrmati, &_mulmati);
-   }
    return gerepileupto(av,y);
  }
  static GEN
  rtran(GEN v, GEN w, GEN q)
  {
-   long av,tetpil;
+   gpmem_t av,tetpil;
    GEN p1;
    if (signe(q))
-Line 169  rtran(GEN v, GEN w, GEN q)
+Line 159  rtran(GEN v, GEN w, GEN q)
 Line 169  rtran(GEN v, GEN w, GEN q)
 Line 159  rtran(GEN v, GEN w, GEN q)
  /* return (v - qw) mod m (only compute entries k0,..,n)
   * v and w are expected to have entries smaller than m */
  static GEN
- mtran(GEN v, GEN w, GEN q, GEN m, long k0)
+ mtran(GEN v, GEN w, GEN q, GEN m, GEN mo2, long k0)
  {
-   long k,l;
+   long k;
    GEN p1;
    if (signe(q))
-   {
+     for (k=lg(v)-1; k >= k0; k--)
-     l = lgefint(m) << 2;
-     for (k=lg(v)-1; k>= k0; k--)
      {
-       long av = avma; (void)new_chunk(l);
+       gpmem_t av = avma;
        p1 = subii((GEN)v[k], mulii(q,(GEN)w[k]));
-       avma = av; v[k]=(long)rrmdr(p1, m);
+       p1 = centermodii(p1, m, mo2);
+       v[k] = lpileuptoint(av, p1);
      }
-   }
    return v;
  }
-Line 236  static void
+Line 224  static void
 Line 236  static void
 Line 224  static void
  rowred(GEN a, GEN rmod)
  {
    long j,k,pro, c = lg(a), r = lg(a[1]);
-   long av=avma, lim=stack_lim(av,1);
+   gpmem_t av=avma, lim=stack_lim(av,1);
-   GEN q;
+   GEN q, rmodo2 = shifti(rmod,-1);
    for (j=1; j<r; j++)
    {
      for (k=j+1; k<c; k++)
        while (signe(gcoeff(a,j,k)))
        {
-         q=rquot(gcoeff(a,j,j),gcoeff(a,j,k));
+         q=diviiround(gcoeff(a,j,j),gcoeff(a,j,k));
-         pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod, j);
+         pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod,rmodo2, j);
          a[j]=a[k]; a[k]=pro;
        }
      if (signe(gcoeff(a,j,j)) < 0)
        for (k=j; k<r; k++) coeff(a,k,j)=lnegi(gcoeff(a,k,j));
      for (k=1; k<j; k++)
      {
-       q=rquot(gcoeff(a,j,k),gcoeff(a,j,j));
+       q=diviiround(gcoeff(a,j,k),gcoeff(a,j,j));
-       a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod, k);
+       a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod,rmodo2, k);
      }
      if (low_stack(lim, stack_lim(av,1)))
      {
-Line 268  rowred(GEN a, GEN rmod)
+Line 256  rowred(GEN a, GEN rmod)
 Line 268  rowred(GEN a, GEN rmod)
 Line 256  rowred(GEN a, GEN rmod)
  }
  /* Calcule d/x  ou  d est entier et x matrice triangulaire inferieure
-  * entiere dont les coeff diagonaux divisent d (resultat entier).
+  * entiere dont les coeff diagonaux divisent d (resultat entier). */
-  */
  static GEN
- matinv(GEN x, GEN d, long n)
+ matinv(GEN x, GEN d)
  {
-   long i,j,k,av,av1;
+   gpmem_t av,av1;
+   long i,j,k, n = lg(x[1]); /* Warning: lg(x) from ordmax is bogus */
    GEN y,h;
-   y=idmat(n);
+   y = idmat(n-1);
-   for (i=1; i<=n; i++)
+   for (i=1; i<n; i++)
-     coeff(y,i,i)=ldivii(d,gcoeff(x,i,i));
+     coeff(y,i,i) = (long)diviiexact(d,gcoeff(x,i,i));
    av=avma;
-   for (i=2; i<=n; i++)
+   for (i=2; i<n; i++)
      for (j=i-1; j; j--)
      {
        for (h=gzero,k=j+1; k<=i; k++)
-Line 298  matinv(GEN x, GEN d, long n)
+Line 286  matinv(GEN x, GEN d, long n)
 Line 298  matinv(GEN x, GEN d, long n)
 Line 286  matinv(GEN x, GEN d, long n)
  static GEN
  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
  {
-   long sp,hard_case_exponent,i,n=lg(cf)-1,av=avma, av2,limit;
+   long sp,i,n=lg(cf)-1;
-   GEN T,T2,Tn,m,v,delta, *w;
+   gpmem_t av=avma, av2,limit;
+   GEN T,T2,Tn,m,v,delta,hard_case_exponent, *w;
    const GEN pp = sqri(p);
+   const GEN ppo2 = shifti(pp,-1);
    const long pps = (2*expi(pp)+2<BITS_IN_LONG)? pp[2]: 0;
    if (cmpis(p,n) > 0)
    {
-     hard_case_exponent = 0;
+     hard_case_exponent = NULL;
      sp = 0; /* gcc -Wall */
    }
    else
-Line 313  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 303  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 313  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 303  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
      long k;
      k = sp = itos(p);
      i=1; while (k < n) { k *= sp; i++; }
-     hard_case_exponent = i;
+     hard_case_exponent = stoi(i);
    }
    T=cgetg(n+1,t_MAT); for (i=1; i<=n; i++) T[i]=lgetg(n+1,t_COL);
    T2=cgetg(2*n+1,t_MAT); for (i=1; i<=2*n; i++) T2[i]=lgetg(n+1,t_COL);
    Tn=cgetg(n*n+1,t_MAT); for (i=1; i<=n*n; i++) Tn[i]=lgetg(n+1,t_COL);
    v = new_chunk(n+1);
-   w =  (GEN*)new_chunk(n+1);
+   w = (GEN*)new_chunk(n+1);
    av2 = avma; limit = stack_lim(av2,1);
    delta=gun; m=idmat(n);
    for(;;)
    {
-     long j,k,h, av0 = avma;
+     long j, k, h;
+     gpmem_t av0 = avma;
      GEN t,b,jp,hh,index,p1, dd = sqri(delta), ppdd = mulii(dd,pp);
+     GEN ppddo2 = shifti(ppdd,-1);
      if (DEBUGLEVEL > 3)
        fprintferr("ROUND2: epsilon = %ld\tavma = %ld\n",epsilon,avma);
-     b=matinv(m,delta,n);
+     b=matinv(m,delta);
      for (i=1; i<=n; i++)
      {
        for (j=1; j<=n; j++)
-Line 344  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 336  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 344  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 336  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
              GEN p2 = mulii(gcoeff(m,i,h),gcoeff(cf[h],j,k));
              if (p2!=gzero) p1 = addii(p1,p2);
            }
-           coeff(T,j,k) = (long)rrmdr(p1, ppdd);
+           coeff(T,j,k) = (long)centermodii(p1, ppdd, ppddo2);
          }
        p1 = mulmati(m, mulmati(T,b));
        for (j=1; j<=n; j++)
          for (k=1; k<=n; k++)
-           coeff(p1,j,k)=(long)rrmdr(divii(gcoeff(p1,j,k),dd),pp);
+           coeff(p1,j,k)=(long)centermodii(divii(gcoeff(p1,j,k),dd),pp,ppo2);
        w[i] = p1;
      }
-Line 381  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 373  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 381  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 373  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
        for (i=1; i<=n; i++)
          for (j=1; j<=n; j++)
          {
-           long av1 = avma;
+           gpmem_t av1 = avma;
            p1 = gzero;
            for (k=1; k<=n; k++)
              for (h=1; h<=n; h++)
              {
                const GEN r=modii(gcoeff(w[i],k,h),p);
                const GEN s=modii(gcoeff(w[j],h,k),p);
                const GEN p2 = mulii(r,s);
                if (p2!=gzero) p1 = addii(p1,p2);
              }
            coeff(T,i,j) = lpileupto(av1,p1);
-Line 405  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 397  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 405  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 397  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
            coeff(T2,j,i)=(i==j)? ps: 0;
            coeff(T2,j,n+i)=smodis(gcoeff(t,i,j),ps);
          }
        rowred_long(T2,pps);
      }
      else
      {
-Line 417  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 409  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 417  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 409  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
          }
        rowred(T2,pp);
      }
-     jp=matinv(T2,p,n);
+     jp=matinv(T2,p);
      if (pps)
      {
        for (k=1; k<=n; k++)
        {
-         long av1=avma;
+         gpmem_t av1=avma;
          t = mulmati(mulmati(jp,w[k]), T2);
          for (h=i=1; i<=n; i++)
            for (j=1; j<=n; j++)
-Line 447  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 439  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 447  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 439  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
        index = mulii(index,gcoeff(Tn,i,i));
      if (gcmp1(index)) break;
-     m = mulmati(matinv(Tn,index,n), m);
+     m = mulmati(matinv(Tn,index), m);
      hh = delta = mulii(index,delta);
      for (i=1; i<=n; i++)
        for (j=1; j<=n; j++)
-Line 486  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
+Line 478  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 486  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
 Line 478  ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
   *
   *  2) discriminant of K (in *y).
   */
- GEN
+ static GEN
- allbase(GEN f, long code, GEN *y)
+ allbase2(GEN f, int flag, GEN *dx, GEN *dK, GEN *ptw)
  {
-   GEN w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
+   GEN w,w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
-   long av=avma,tetpil,n,h,j,i,k,r,s,t,v,mf;
+   gpmem_t av=avma,tetpil;
+   long n,h,j,i,k,r,s,t,v,mf;
-   allbase_check_args(f,code,y, &w1,&w2);
+   w = ptw? *ptw: NULL;
+   allbase_check_args(f,flag,dx, &w);
+   w1 = (GEN)w[1];
+   w2 = (GEN)w[2];
    v = varn(f); n = degpol(f); h = lg(w1)-1;
    cf = (GEN*)cgetg(n+1,t_VEC);
    cf[2]=companion(f);
-Line 501  allbase(GEN f, long code, GEN *y)
+Line 497  allbase(GEN f, long code, GEN *y)
 Line 501  allbase(GEN f, long code, GEN *y)
 Line 497  allbase(GEN f, long code, GEN *y)
    a=idmat(n); da=gun;
    for (i=1; i<=h; i++)
    {
-     long av1 = avma;
+     gpmem_t av1 = avma;
      mf=itos((GEN)w2[i]); if (mf==1) continue;
      if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
-Line 513  allbase(GEN f, long code, GEN *y)
+Line 509  allbase(GEN f, long code, GEN *y)
 Line 513  allbase(GEN f, long code, GEN *y)
 Line 509  allbase(GEN f, long code, GEN *y)
        for (s=r; s; s--)
          while (signe(gcoeff(bt,s,r)))
          {
-           q=rquot(gcoeff(at,s,s),gcoeff(bt,s,r));
+           q=diviiround(gcoeff(at,s,s),gcoeff(bt,s,r));
            pro=rtran((GEN)at[s],(GEN)bt[r],q);
            for (t=s-1; t; t--)
            {
-             q=rquot(gcoeff(at,t,s),gcoeff(at,t,t));
+             q=diviiround(gcoeff(at,t,s),gcoeff(at,t,t));
              pro=rtran(pro,(GEN)at[t],q);
            }
            at[s]=bt[r]; bt[r]=(long)pro;
-Line 528  allbase(GEN f, long code, GEN *y)
+Line 524  allbase(GEN f, long code, GEN *y)
 Line 528  allbase(GEN f, long code, GEN *y)
 Line 524  allbase(GEN f, long code, GEN *y)
        {
          while (signe(gcoeff(at,j,k)))
          {
-           q=rquot(gcoeff(at,j,j),gcoeff(at,j,k));
+           q=diviiround(gcoeff(at,j,j),gcoeff(at,j,k));
            pro=rtran((GEN)at[j],(GEN)at[k],q);
            at[j]=at[k]; at[k]=(long)pro;
          }
-Line 537  allbase(GEN f, long code, GEN *y)
+Line 533  allbase(GEN f, long code, GEN *y)
 Line 537  allbase(GEN f, long code, GEN *y)
 Line 533  allbase(GEN f, long code, GEN *y)
          for (k=1; k<=j; k++) coeff(at,k,j)=lnegi(gcoeff(at,k,j));
        for (k=j+1; k<=n; k++)
        {
-         q=rquot(gcoeff(at,j,k),gcoeff(at,j,j));
+         q=diviiround(gcoeff(at,j,k),gcoeff(at,j,j));
          at[k]=(long)rtran((GEN)at[k],(GEN)at[j],q);
        }
      }
-Line 549  allbase(GEN f, long code, GEN *y)
+Line 545  allbase(GEN f, long code, GEN *y)
 Line 549  allbase(GEN f, long code, GEN *y)
 Line 545  allbase(GEN f, long code, GEN *y)
        }
      tetpil=avma; a=gtrans(at);
      {
        GEN *gptr[2];
        da = icopy(da); gptr[0]=&a; gptr[1]=&da;
        gerepilemanysp(av1,tetpil,gptr,2);
      }
    }
+   *dK = *dx;
    for (j=1; j<=n; j++)
-     *y = divii(mulii(*y,sqri(gcoeff(a,j,j))), sqri(da));
+     *dK = diviiexact(mulii(*dK,sqri(gcoeff(a,j,j))), sqri(da));
-   tetpil=avma; *y=icopy(*y);
+   tetpil=avma; *dK = icopy(*dK);
    at=cgetg(n+1,t_VEC); v=varn(f);
    for (k=1; k<=n; k++)
    {
      q=cgetg(k+2,t_POL); at[k]=(long)q;
      q[1] = evalsigne(1) | evallgef(2+k) | evalvarn(v);
-     for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,k,j),da);
+     for (j=1; j<=k; j++) q[j+1] = ldiv(gcoeff(a,k,j),da);
    }
-   gptr[0]=&at; gptr[1]=y;
+   gptr[0] = &at; gptr[1] = dK;
    gerepilemanysp(av,tetpil,gptr,2);
    return at;
  }
  GEN
- base2(GEN x, GEN *y)
+ base2(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_ROUND2, NULL); }
- {
-   return allbase(x,0,y);
- }
  GEN
- discf2(GEN x)
+ discf2(GEN x) { return nfdiscf0(x, compat_ROUND2, NULL); }
- {
-   GEN y;
-   long av=avma,tetpil;
-   allbase(x,0,&y); tetpil=avma;
-   return gerepile(av,tetpil,icopy(y));
- }
  /*******************************************************************/
  /*                                                                 */
  /*                            ROUND 4                              */
-Line 595  GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
+Line 582  GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
 Line 595  GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
 Line 582  GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
  GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
  static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
  static GEN maxord(GEN p,GEN f,long mf);
- static GEN nbasis(GEN ibas,GEN pd);
  static GEN testb2(GEN p,GEN fa,long Fa,GEN theta,GEN pmf,long Ft,GEN ns);
  static GEN testc2(GEN p,GEN fa,GEN pmr,GEN pmf,GEN alph2,
                    long Ea,GEN thet2,long Et,GEN ns);
-Line 609  fnz(GEN x,long j)
+Line 595  fnz(GEN x,long j)
 Line 609  fnz(GEN x,long j)
 Line 595  fnz(GEN x,long j)
    return 1;
  }
- /* retourne la base, dans y le discf et dans ptw la factorisation (peut
+ /* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
-  etre partielle) de discf */
+ static GEN
+ get_coprimes(GEN a, GEN b)
+ {
+   long i, k = 1;
+   GEN u = cgetg(3, t_VEC);
+   u[1] = (long)a;
+   u[2] = (long)b;
+   /* u1,..., uk 2 by 2 coprime */
+   while (k+1 < lg(u))
+   {
+     GEN d, c = (GEN)u[k+1];
+     if (is_pm1(c)) { k++; continue; }
+     for (i=1; i<=k; i++)
+     {
+       if (is_pm1(u[i])) continue;
+       d = mppgcd(c, (GEN)u[i]);
+       if (d == gun) continue;
+       c = diviiexact(c, d);
+       u[i] = (long)diviiexact((GEN)u[i], d);
+       u = concatsp(u, d);
+     }
+     u[++k] = (long)c;
+   }
+   for (i = k = 1; i < lg(u); i++)
+     if (!is_pm1(u[i])) u[k++] = u[i];
+   setlg(u, k); return u;
+ }
+ /* return integer basis. Set dK = disc(K), dx = disc(f), w (possibly partial)
+  * factorization of dK. *ptw can be set by the caller, in which case it is
+  * taken to be the factorization of disc(f), then overwritten
+  * [no consistency check] */
  GEN
- allbase4(GEN f,long code, GEN *y, GEN *ptw)
+ allbase(GEN f, int flag, GEN *dx, GEN *dK, GEN *index, GEN *ptw)
  {
-   GEN w,w1,w2,a,da,b,db,bas,q,p1,*gptr[3];
+   GEN w, w1, w2, a, da, p1, ordmax;
-   long v,n,mf,h,lfa,i,j,k,l,tetpil,av = avma;
+   long n, lw, i, j, k, l;
-   allbase_check_args(f,code,y, &w1,&w2);
+   if (flag & nf_ROUND2) return allbase2(f,flag,dx,dK,ptw);
-   v = varn(f); n = degpol(f); h = lg(w1)-1;
+   w = ptw? *ptw: NULL;
-   a = NULL; /* gcc -Wall */
+   allbase_check_args(f, flag, dx, &w);
-   da= NULL;
+   w1 = (GEN)w[1];
-   for (i=1; i<=h; i++)
+   w2 = (GEN)w[2];
+   n = degpol(f); lw = lg(w1);
+   ordmax = cgetg(1, t_VEC);
+   /* "complete" factorization first */
+   for (i=1; i<lw; i++)
    {
-     mf=itos((GEN)w2[i]); if (mf == 1) continue;
+     long mf = itos((GEN)w2[i]);
-     if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
+     if (mf == 1) { ordmax = concatsp(ordmax, gun); continue; }
-     b = maxord((GEN)w1[i],f,mf); db = gun;
+     CATCH(invmoder) { /* caught false prime, update factorization */
+       GEN x = (GEN)global_err_data;
+       GEN p = mppgcd((GEN)x[1], (GEN)x[2]);
+       GEN N, u;
+       if (DEBUGLEVEL) err(warner,"impossible inverse: %Z", x);
+       u = get_coprimes(p, diviiexact((GEN)x[1],p));
+       l = lg(u);
+       /* no small factors, but often a prime power */
+       for (k = 1; k < l; k++) u[k] = coeff(auxdecomp((GEN)u[k], 2),1,1);
+       w1[i] = u[1];
+       w1 = concatsp(w1, vecextract_i(u, 2, l-1));
+       N = *dx;
+       w2[i] = lstoi(pvaluation(N, (GEN)w1[i], &N));
+       k  = lw;
+       lw = lg(w1);
+       for ( ; k < lw; k++) w2[k] = lstoi(pvaluation(N, (GEN)w1[k], &N));
+     } RETRY {
+       if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
+       ordmax = concatsp(ordmax, _vec( maxord((GEN)w1[i],f,mf) ));
+     } ENDCATCH;
+   }
+   a = NULL; /* gcc -Wall */
+   da = NULL;
+   for (i=1; i<lw; i++)
+   {
+     GEN db, b = (GEN)ordmax[i];
+     if (b == gun) continue;
+     db = gun;
      for (j=1; j<=n; j++)
      {
        p1 = denom(gcoeff(b,j,j));
        if (cmpii(p1,db) > 0) db = p1;
      }
-     if (db != gun)
+     if (db == gun) continue;
-     { /* db = denom(diag(b)), (da,db) = 1 */
-       b = gmul(b,db);
+     /* db = denom(diag(b)), (da,db) = 1 */
-       if (!da) { da=db; a=b; }
+     b = Q_muli_to_int(b,db);
-       else
+     if (!da) { da = db; a = b; }
+     else
+     {
+       j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
+       b = gmul(da,b);
+       a = gmul(db,a); da = mulii(da,db);
+       k = j-1; p1 = cgetg(2*n-k+1,t_MAT);
+       for (j=1; j<=k; j++)
        {
-         j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
+         p1[j] = a[j];
-         b = gmul(da,b); a = gmul(db,a);
+         coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
-         k=j-1; p1=cgetg(2*n-k+1,t_MAT);
-         for (j=1; j<=k; j++)
-         {
-           p1[j] = a[j];
-           coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
-         }
-         for (  ; j<=n;     j++) p1[j] = a[j];
-         for (  ; j<=2*n-k; j++) p1[j] = b[j+k-n];
-         da = mulii(da,db); a = hnfmodid(p1, da);
        }
+       for (  ; j<=n;     j++) p1[j] = a[j];
+       for (  ; j<=2*n-k; j++) p1[j] = b[j+k-n];
+       a = hnfmodid(p1, da);
      }
-     if (DEBUGLEVEL>5)
+     if (DEBUGLEVEL>5) fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
-       fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
    }
+   *dK = *dx;
    if (da)
    {
-     for (j=1; j<=n; j++)
+     *index = diviiexact(da, gcoeff(a,1,1));
-       *y = mulii(divii(*y,sqri(da)),sqri(gcoeff(a,j,j)));
+     for (j=2; j<=n; j++) *index = mulii(*index, diviiexact(da, gcoeff(a,j,j)));
+     *dK = diviiexact(*dK, sqri(*index));
      for (j=n-1; j; j--)
        if (cmpis(gcoeff(a,j,j),2) > 0)
        {
-         p1=shifti(gcoeff(a,j,j),-1);
+         p1 = shifti(gcoeff(a,j,j),-1);
          for (k=j+1; k<=n; k++)
            if (cmpii(gcoeff(a,j,k),p1) > 0)
              for (l=1; l<=j; l++)
-               coeff(a,l,k)=lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
+               coeff(a,l,k) = lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
        }
+     a = gdiv(a, da);
    }
-   lfa = 0;
+   else
-   if (ptw)
    {
-     for (j=1; j<=h; j++)
+     *index = gun;
-     {
+     a = idmat(n);
-       k=ggval(*y,(GEN)w1[j]);
-       if (k) { lfa++; w1[lfa]=w1[j]; w2[lfa]=k; }
-     }
    }
-   tetpil=avma; *y=icopy(*y);
-   bas=cgetg(n+1,t_VEC); v=varn(f);
-   for (k=1; k<=n; k++)
-   {
-     q=cgetg(k+2,t_POL); bas[k]=(long)q;
-     q[1] = evalsigne(1) | evallgef(k+2) | evalvarn(v);
-     if (da)
-       for (j=1; j<=k; j++) q[j+1]=ldiv(gcoeff(a,j,k),da);
-     else
-     {
-       for (j=2; j<=k; j++) q[j]=zero;
-       q[j]=un;
-     }
-   }
    if (ptw)
    {
-     *ptw=w=cgetg(3,t_MAT);
+     long lfa = 1;
-     w[1]=lgetg(lfa+1,t_COL);
+     GEN W1, W2, D = *dK;
-     w[2]=lgetg(lfa+1,t_COL);
-     for (j=1; j<=lfa; j++)
+     w = cgetg(3,t_MAT);
+     W1 = cgetg(lw, t_COL); w[1] = (long)W1;
+     W2 = cgetg(lw, t_COL); w[2] = (long)W2;
+     for (j=1; j<lw; j++)
      {
-       coeff(w,j,1)=(long)icopy((GEN)w1[j]);
+       k = pvaluation(D, (GEN)w1[j], &D);
-       coeff(w,j,2)=lstoi(w2[j]);
+       if (k) { W1[lfa] = w1[j]; W2[lfa] = lstoi(k); lfa++; }
      }
-     gptr[2]=ptw;
+     setlg(W1, lfa);
+     setlg(W2, lfa); *ptw = w;
    }
-   gptr[0]=&bas; gptr[1]=y;
+   return mat_to_vecpol(a, varn(f));
-   gerepilemanysp(av,tetpil,gptr, ptw?3:2);
-   return bas;
  }
- extern GEN merge_factor_i(GEN f, GEN g);
  static GEN
  update_fact(GEN x, GEN f)
  {
-Line 723  update_fact(GEN x, GEN f)
+Line 760  update_fact(GEN x, GEN f)
 Line 723  update_fact(GEN x, GEN f)
 Line 760  update_fact(GEN x, GEN f)
    for (i=1; i<l; i++)
    {
      k = pvaluation(d, (GEN)p[i], &d);
      if (k) { q[iq] = p[i]; e[iq] = lstoi(k); iq++; }
    }
    setlg(q,iq); setlg(e,iq);
    return merge_factor_i(decomp(d), g);
  }
- /* if y is non-NULL, it receives the discriminant
-  * return basis if (ret_basis != 0), discriminant otherwise
-  */
  static GEN
- nfbasis00(GEN x0, long flag, GEN p, long ret_basis, GEN *y)
+ unscale_vecpol(GEN v, GEN h)
  {
-   GEN x, disc, basis, lead;
+   long i, l;
-   GEN *gptr[2];
+   GEN w;
-   long k, tetpil, av = avma, l = lgef(x0), smll;
+   if (!h) return v;
+   l = lg(v); w = cgetg(l, typ(v));
+   for (i=1; i<l; i++) w[i] = (long)unscale_pol((GEN)v[i], h);
+   return w;
+ }
-   if (typ(x0)!=t_POL) err(typeer,"nfbasis00");
+ /* FIXME: have to deal with compatibility flags */
-   if (l<=3) err(zeropoler,"nfbasis00");
+ static void
-   for (k=2; k<l; k++)
+ _nfbasis(GEN x0, long flag, GEN fa, GEN *pbas, GEN *pdK)
-     if (typ(x0[k])!=t_INT) err(talker,"polynomial not in Z[X] in nfbasis");
+ {
+   GEN x, dx, dK, basis, lead, index;
+   long fl = 0;
-   x = pol_to_monic(x0,&lead);
+   if (typ(x0)!=t_POL) err(typeer,"nfbasis");
+   if (!degpol(x0)) err(zeropoler,"nfbasis");
+   check_pol_int(x0, "nfbasis");
-   if (!p || gcmp0(p))
+   x = pol_to_monic(x0, &lead);
-     smll = (flag & 1); /* small basis */
+   if (fa && gcmp0(fa)) fa = NULL; /* compatibility. NULL is the proper arg */
-   else
+   if (fa && lead) fa = update_fact(x, fa);
-   {
+   if (flag & compat_PARTIAL) fl |= nf_PARTIALFACT;
-     if (lead) p = update_fact(x,p);
+   if (flag & compat_ROUND2)  fl |= nf_ROUND2;
-     smll = (long) p;   /* factored basis */
+   basis = allbase(x, fl, &dx, &dK, &index, &fa);
-   }
+   if (pbas) *pbas = unscale_vecpol(basis, lead);
+   if (pdK)  *pdK = dK;
-   if (flag & 2)
-     basis = allbase(x,smll,&disc); /* round 2 */
-   else
-     basis = allbase4(x,smll,&disc,NULL); /* round 4 */
-   if (!ret_basis) return gerepilecopy(av,disc);
-   tetpil=avma;
-   if (!lead) basis = gcopy(basis);
-   else
-   {
-     long v = varn(x);
-     GEN pol = gmul(polx[v],lead);
-     tetpil = avma; basis = gsubst(basis,v,pol);
-   }
-   if (!y)
-     return gerepile(av,tetpil,basis);
-   *y = gcopy(disc);
-   gptr[0]=&basis; gptr[1]=y;
-   gerepilemanysp(av,tetpil,gptr,2);
-   return basis;
  }
  GEN
- nfbasis(GEN x, GEN *y, long flag, GEN p)
+ nfbasis(GEN x, GEN *pdK, long flag, GEN fa)
  {
-   return nfbasis00(x,flag,p,1,y);
+   gpmem_t av = avma;
+   GEN bas; _nfbasis(x, flag, fa, &bas, pdK);
+   gerepileall(av, 2, &bas, pdK); return bas;
  }
  GEN
- nfbasis0(GEN x, long flag, GEN p)
+ nfbasis0(GEN x, long flag, GEN fa)
  {
-   return nfbasis00(x,flag,p,1,NULL);
+   gpmem_t av = avma;
+   GEN bas; _nfbasis(x, flag, fa, &bas, NULL);
+   return gerepilecopy(av, bas);
  }
  GEN
- nfdiscf0(GEN x, long flag, GEN p)
+ nfdiscf0(GEN x, long flag, GEN fa)
  {
-   return nfbasis00(x,flag,p,0,&p);
+   gpmem_t av = avma;
+   GEN dK; _nfbasis(x, flag, fa, NULL, &dK);
+   return gerepilecopy(av, dK);
  }
  GEN
- base(GEN x, GEN *y)
+ base(GEN x, GEN *pdK) { return nfbasis(x, pdK, 0, NULL); }
- {
-   return allbase4(x,0,y,NULL);
- }
  GEN
- smallbase(GEN x, GEN *y)
+ smallbase(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_PARTIAL, NULL); }
- {
-   return allbase4(x,1,y,NULL);
- }
  GEN
- factoredbase(GEN x, GEN p, GEN *y)
+ factoredbase(GEN x, GEN fa, GEN *pdK) { return nfbasis(x, pdK, 0, fa); }
- {
-   return allbase4(x,(long)p,y,NULL);
- }
  GEN
- discf(GEN x)
+ discf(GEN x) { return nfdiscf0(x, 0, NULL); }
- {
-   GEN y;
-   long av=avma,tetpil;
-   allbase4(x,0,&y,NULL); tetpil=avma;
-   return gerepile(av,tetpil,icopy(y));
- }
  GEN
- smalldiscf(GEN x)
+ smalldiscf(GEN x) { return nfdiscf0(x, nf_PARTIALFACT, NULL); }
- {
-   GEN y;
-   long av=avma,tetpil;
-   allbase4(x,1,&y,NULL); tetpil=avma;
-   return gerepile(av,tetpil,icopy(y));
- }
  GEN
- factoreddiscf(GEN x, GEN p)
+ factoreddiscf(GEN x, GEN fa) { return nfdiscf0(x, 0, fa); }
- {
-   GEN y;
-   long av=avma,tetpil;
-   allbase4(x,(long)p,&y,NULL); tetpil=avma;
-   return gerepile(av,tetpil,icopy(y));
- }
  /* return U if Z[alpha] is not maximal or 2*dU < m-1; else return NULL */
  static GEN
-Line 852  dedek(GEN f, long mf, GEN p,GEN g)
+Line 848  dedek(GEN f, long mf, GEN p,GEN g)
 Line 852  dedek(GEN f, long mf, GEN p,GEN g)
 Line 848  dedek(GEN f, long mf, GEN p,GEN g)
    GEN k,h;
    long dk;
-   if (DEBUGLEVEL>=3)
-   {
-     fprintferr("  entering dedek ");
-     if (DEBUGLEVEL>5)
-       fprintferr("with parameters p=%Z,\n  f=%Z",p,f);
-     fprintferr("\n");
-   }
    h = FpX_div(f,g,p);
    k = gdivexact(gadd(f, gneg_i(gmul(g,h))), p);
    k = FpX_gcd(k, FpX_gcd(g,h, p), p);
    dk = degpol(k);
-   if (DEBUGLEVEL>=3) fprintferr("  gcd has degree %ld\n", dk);
+   if (DEBUGLEVEL>2)
+   {
+     fprintferr("  dedek: gcd has degree %ld\n", dk);
+     if (DEBUGLEVEL>5) fprintferr("initial parameters p=%Z,\n  f=%Z\n",p,f);
+   }
    if (2*dk >= mf-1) return FpX_div(f,k,p);
    return dk? (GEN)NULL: f;
  }
-Line 873  dedek(GEN f, long mf, GEN p,GEN g)
+Line 866  dedek(GEN f, long mf, GEN p,GEN g)
 Line 873  dedek(GEN f, long mf, GEN p,GEN g)
 Line 866  dedek(GEN f, long mf, GEN p,GEN g)
  static GEN
  maxord(GEN p,GEN f,long mf)
  {
-   long j,r, av = avma, flw = (cmpsi(degpol(f),p) < 0);
+   const gpmem_t av = avma;
+   long j,r, flw = (cmpsi(degpol(f),p) < 0);
    GEN w,g,h,res;
    if (flw)
-Line 902  maxord(GEN p,GEN f,long mf)
+Line 896  maxord(GEN p,GEN f,long mf)
 Line 902  maxord(GEN p,GEN f,long mf)
 Line 896  maxord(GEN p,GEN f,long mf)
  static GEN
  polmodiaux(GEN x, GEN y, GEN ys2)
  {
-   if (typ(x)!=t_INT)
+   if (typ(x)!=t_INT) x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
-     x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
+   return centermodii(x,y,ys2);
-   x = modii(x,y);
-   if (cmpii(x,ys2) > 0) x = subii(x,y);
-   return x;
  }
  /* x polynomial with integer or rational coeff. Reduce them mod y IN PLACE */
-Line 931  polmodi_keep(GEN x, GEN y)
+Line 922  polmodi_keep(GEN x, GEN y)
 Line 931  polmodi_keep(GEN x, GEN y)
 Line 922  polmodi_keep(GEN x, GEN y)
  }
  static GEN
+ shiftpol(GEN x, long v)
+ {
+   x = addshiftw(x, zeropol(v), 1);
+   setvarn(x,v); return x;
+ }
+ /* Sylvester's matrix, mod p^m (assumes f1 monic) */
+ static GEN
+ sylpm(GEN f1,GEN f2,GEN pm)
+ {
+   long n,j,v=varn(f1);
+   GEN a,h;
+   n=degpol(f1); a=cgetg(n+1,t_MAT);
+   h = FpX_res(f2,f1,pm);
+   for (j=1;; j++)
+   {
+     a[j] = (long)pol_to_vec(h,n);
+     if (j == n) break;
+     h = FpX_res(shiftpol(h,v),f1,pm);
+   }
+   return hnfmodid(a,pm);
+ }
+ /* polynomial gcd mod p^m (assumes f1 monic) */
+ GEN
+ gcdpm(GEN f1,GEN f2,GEN pm)
+ {
+   gpmem_t av=avma,tetpil;
+   long n,c,v=varn(f1);
+   GEN a,col;
+   n=degpol(f1); a=sylpm(f1,f2,pm);
+   for (c=1; c<=n; c++)
+     if (signe(resii(gcoeff(a,c,c),pm))) break;
+   if (c > n) { avma=av; return zeropol(v); }
+   col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
+   return gerepile(av,tetpil, gtopolyrev(col,v));
+ }
+ /* reduced resultant mod p^m (assumes x monic) */
+ GEN
+ respm(GEN x,GEN y,GEN pm)
+ {
+   gpmem_t av = avma;
+   GEN p1 = sylpm(x,y,pm);
+   p1 = gcoeff(p1,1,1);
+   if (egalii(p1,pm)) { avma = av; return gzero; }
+   return gerepileuptoint(av, icopy(p1));
+ }
+ static GEN
  dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U)
  {
-   long n=degpol(f),dU,c;
+   long n=degpol(f),dU,i;
    GEN b,ha,pd,pdp;
    if (n == 1) return gscalmat(gun, 1);
-   if (DEBUGLEVEL>=3)
+   if (DEBUGLEVEL>5)
    {
-     fprintferr("  entering Dedekind Basis ");
+     fprintferr("  entering Dedekind Basis with parameters p=%Z\n",p);
-     if (DEBUGLEVEL>5)
+     fprintferr("  f = %Z,\n  alpha = %Z\n",f,alpha);
-     {
-       fprintferr("with parameters p=%Z\n",p);
-       fprintferr("  f = %Z,\n  alpha = %Z",f,alpha);
-     }
-     fprintferr("\n");
    }
-   ha = pd = gpuigs(p,mf/2); pdp = mulii(pd,p);
+   ha = pd = gpowgs(p,mf/2); pdp = mulii(pd,p);
    dU = typ(U)==t_POL? degpol(U): 0;
    b = cgetg(n,t_MAT); /* Z[a] + U/p Z[a] is maximal */
    /* skip first column = gscalcol(pd,n) */
-   for (c=1; c<n; c++)
+   for (i=1; i<n; i++)
    {
-     if (c == dU)
+     if (i == dU)
      {
-       ha = gdiv(gmul(pd,eleval(f,U,alpha)),p);
+       ha = gdiv(gmul(pd,RX_RXQ_compo(U,alpha,f)),p);
        ha = polmodi(ha,pdp);
      }
      else
      {
-       GEN p2, mod;
+       GEN d, mod;
        ha = gmul(ha,alpha);
-       p2 = content(ha); /* to cancel denominator */
+       ha = Q_remove_denom(ha, &d);
-       if (gcmp1(p2)) { p2 = NULL; mod = pdp; }
+       mod = d? mulii(pdp,d): pdp;
-       else
-       {
-         ha = gdiv(ha,p2);
-         if (typ(p2)==t_INT)
-           mod = divii(pdp, mppgcd(pdp,p2));
-         else
-           mod = mulii(pdp, (GEN)p2[2]); /* p2 = a / p^e */
-       }
        ha = FpX_res(ha, f, mod);
-       if (p2) ha = gmul(ha,p2);
+       if (d) ha = gdivexact(ha,d);
      }
-     b[c] = (long)pol_to_vec(ha,n);
+     b[i] = (long)pol_to_vec(ha,n);
    }
    b = hnfmodid(b,pd);
    if (DEBUGLEVEL>5) fprintferr("  new order: %Z\n",b);
-   return gdiv(b,pd);
+   return gdiv(b, pd);
  }
  static GEN
  get_partial_order_as_pols(GEN p, GEN f)
  {
-   long i,j, n = degpol(f), vf = varn(f);
+   GEN b = maxord(p,f, ggval(ZX_disc(f),p));
-   GEN b,ib,h,col;
+   return mat_to_vecpol(b, varn(f));
+ }
-   b = maxord(p,f, ggval(ZX_disc(f),p));
+ long
-   ib = cgetg(n+1,t_VEC);
+ FpX_val(GEN x0, GEN t, GEN p, GEN *py)
-   for (i=1; i<=n; i++)
+ {
+   long k;
+   GEN r, y, x = x0;
+   for (k=0; ; k++)
    {
-     h=cgetg(i+2,t_POL); ib[i]=(long)h; col=(GEN)b[i];
+     y = FpX_divres(x, t, p, &r);
-     h[1]=evalsigne(1)|evallgef(i+2)|evalvarn(vf);
+     if (signe(r)) break;
-     for (j=1;j<=i;j++) h[j+1]=col[j];
+     x = y;
    }
-   return ib;
+   *py = x; return k;
  }
+ /* e in Qp, f i Zp. Return r = e mod (f, pk) */
+ static GEN
+ QpX_mod(GEN e, GEN f, GEN pk)
+ {
+   GEN mod, d;
+   e = Q_remove_denom(e, &d);
+   mod = d? mulii(pk,d): pk;
+   e = FpX_res(e, centermod(f, mod), mod);
+   e = FpX_center(e, mod);
+   if (d) e = gdiv(e, d);
+   return e;
+ }
  /* if flag != 0, factorization to precision r (maximal order otherwise) */
  GEN
  Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long flag)
  {
-   GEN res,pr,pk,ph,pdr,unmodp,b1,b2,b3,a1,e,f1,f2;
+   GEN fred,res,pr,pk,ph,pdr,b1,b2,a,e,f1,f2;
    if (DEBUGLEVEL>2)
    {
-Line 1016  Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,lo
+Line 1066  Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,lo
 Line 1016  Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,lo
 Line 1066  Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,lo
      }
      fprintferr("\n");
    }
-   unmodp = gmodulsg(1,p);
-   b1=lift_intern(gmul(chi,unmodp));
+   (void)FpX_val(chi, nu, p, &b1); /* nu irreducible mod p */
-   a1=gun; b2=gun;
+   b2 = FpX_div(chi, b1, p);
-   b3=lift_intern(gmul(nu,unmodp));
+   a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
-   while (degpol(b3) > 0)
+   pdr = respm(f, derivpol(f), gpowgs(p,mf+1));
-   {
+   e = RX_RXQ_compo(a, theta, f);
-     GEN p1;
-     b1 = FpX_div(b1,b3, p);
-     b2 = FpX_red(gmul(b2,b3), p);
-     b3 = FpX_extgcd(b2,b1, p, &a1,&p1); /* p1 = junk */
-     p1 = leading_term(b3);
-     if (!gcmp1(p1))
-     { /* FpX_extgcd does not return normalized gcd */
-       p1 = mpinvmod(p1,p);
-       b3 = gmul(b3,p1);
-       a1 = gmul(a1,p1);
-     }
-   }
-   pdr = respm(f,derivpol(f),gpuigs(p,mf+1));
-   e = eleval(f,FpX_red(gmul(a1,b2), p),theta);
    e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,p)),pdr);
    pr = flag? gpowgs(p,flag): mulii(p,sqri(pdr));
-   pk=p; ph=mulii(pdr,pr);
+   pk = p; ph = mulii(pdr,pr);
    /* E(t) - e(t) belongs to p^k Op, which is contained in p^(k-df)*Zp[xi] */
    while (cmpii(pk,ph) < 0)
-   {
+   { /* e <-- e^2(3-2e) mod p^2k */
+     pk = sqri(pk);
      e = gmul(gsqr(e), gsubsg(3,gmul2n(e,1)));
-     e = gres(e,f); pk = sqri(pk);
+     e = QpX_mod(e, f, pk);
-     e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,pk)), pdr);
    }
-   f1 = gcdpm(f,gmul(pdr,gsubsg(1,e)), ph);
+   fred = centermod(f, ph);
-   f1 = FpX_res(f1,f, pr);
+   f1 = gcdpm(fred, gmul(pdr,gsubsg(1,e)), ph);
-   f2 = FpX_res(FpX_div(f,f1, pr), f, pr);
+   fred = centermod(fred, pr); /* pr | ph */
+   f1 = centermod(f1, pr);
+   f2 = FpX_div(fred,f1, pr);
+   f2 = FpX_center(f2, pr);
-   if (DEBUGLEVEL>2)
+   if (DEBUGLEVEL>5)
    {
      fprintferr("  leaving Decomp");
-     if (DEBUGLEVEL>5)
+     fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n\n", f1,f2,e);
-       fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n", f1,f2,e);
-     fprintferr("\n");
    }
    if (flag)
    {
-     b1=factorpadic4(f1,p,flag);
+     b1 = factorpadic4(f1,p,flag);
-     b2=factorpadic4(f2,p,flag); res=cgetg(3,t_MAT);
+     b2 = factorpadic4(f2,p,flag); res = cgetg(3,t_MAT);
-     res[1]=lconcat((GEN)b1[1],(GEN)b2[1]);
+     res[1] = (long)concatsp((GEN)b1[1],(GEN)b2[1]);
-     res[2]=lconcat((GEN)b1[2],(GEN)b2[2]); return res;
+     res[2] = (long)concatsp((GEN)b1[2],(GEN)b2[2]); return res;
    }
    else
    {
-     GEN ib1,ib2;
+     GEN ib1, ib2;
-     long n1,n2,i;
+     long n, n1, n2, i;
-     ib1 = get_partial_order_as_pols(p,f1); n1=lg(ib1)-1;
+     ib1 = get_partial_order_as_pols(p,f1); n1 = lg(ib1)-1;
-     ib2 = get_partial_order_as_pols(p,f2); n2=lg(ib2)-1;
+     ib2 = get_partial_order_as_pols(p,f2); n2 = lg(ib2)-1;
-     res=cgetg(n1+n2+1,t_VEC);
+     n = n1+n2;
+     res = cgetg(n+1, t_VEC);
      for (i=1; i<=n1; i++)
-       res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib1[i]),e),f), pdr);
+       res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib1[i]),e), f, pdr);
-     e=gsubsg(1,e); ib2 -= n1;
+     e = gsubsg(1,e); ib2 -= n1;
-     for (   ; i<=n1+n2; i++)
+     for (   ; i<=n; i++)
-       res[i]=(long)polmodi(gmod(gmul(gmul(pdr,(GEN)ib2[i]),e),f), pdr);
+       res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib2[i]),e), f, pdr);
-     return nbasis(res,pdr);
+     res = vecpol_to_mat(res, n);
+     return gdiv(hnfmodid(res,pdr), pdr); /* normalized integral basis */
    }
  }
  /* minimum extension valuation: res[0]/res[1] (both are longs) */
- static long *
+ static void
- vstar(GEN p,GEN h)
+ vstar(GEN p,GEN h, long *L, long *E)
  {
-   static long res[2];
    long m,first,j,k,v,w;
    m=degpol(h); first=1; k=1; v=0;
-Line 1098  vstar(GEN p,GEN h)
+Line 1136  vstar(GEN p,GEN h)
 Line 1098  vstar(GEN p,GEN h)
 Line 1136  vstar(GEN p,GEN h)
        first=0;
      }
    m = cgcd(v,k);
-   res[0]=v/m; res[1]=k/m; return res;
+   *L = v/m;
+   *E = k/m;
  }
  /* reduce the element elt modulo rd, taking first care of the denominators */
-Line 1107  redelt(GEN elt, GEN rd, GEN pd)
+Line 1146  redelt(GEN elt, GEN rd, GEN pd)
 Line 1107  redelt(GEN elt, GEN rd, GEN pd)
 Line 1146  redelt(GEN elt, GEN rd, GEN pd)
  {
    GEN den, nelt, nrd, relt;
-   den  = ggcd(denom(content(elt)), pd);
+   den  = ggcd(Q_denom(elt), pd);
    nelt = gmul(den, elt);
    nrd  = gmul(den, rd);
-Line 1123  redelt(GEN elt, GEN rd, GEN pd)
+Line 1162  redelt(GEN elt, GEN rd, GEN pd)
 Line 1123  redelt(GEN elt, GEN rd, GEN pd)
 Line 1162  redelt(GEN elt, GEN rd, GEN pd)
  GEN
  polsymmodpp(GEN g, GEN pp)
  {
-   long av1, av2, d = degpol(g), i, k;
+   gpmem_t av1, av2;
+   long d = degpol(g), i, k;
    GEN s , y;
    y = cgetg(d + 1, t_COL);
    y[1] = lstoi(d);
    for (k = 1; k < d; k++)
    {
      av1 = avma;
      s = centermod(gmulsg(k, polcoeff0(g,d-k,-1)), pp);
      for (i = 1; i < k; i++)
        s = gadd(s, gmul((GEN)y[k-i+1], polcoeff0(g,d-i,-1)));
      av2 = avma;
      y[k + 1] = lpile(av1, av2, centermod(gneg(s), pp));
    }
-Line 1152  manage_cache(GEN chi, GEN pp, GEN ns)
+Line 1192  manage_cache(GEN chi, GEN pp, GEN ns)
 Line 1152  manage_cache(GEN chi, GEN pp, GEN ns)
 Line 1192  manage_cache(GEN chi, GEN pp, GEN ns)
    {
      if (DEBUGLEVEL > 4)
        fprintferr("newtonsums: result too large to fit in cache\n");
      return polsymmodpp(chi, pp);
    }
    if (!signe((GEN)ns[1]))
    {
      ns2 = polsymmodpp(chi, pp);
      for (j = 1; j <= n; j++)
        affii((GEN)ns2[j], (GEN)ns[j]);
    }
    return ns;
  }
  /* compute the Newton sums modulo pp of the characteristic
     polynomial of a(x) mod g(x) */
  static GEN
  newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
  {
    GEN va, pa, s, ns2;
-   long j, k, n = degpol(chi), av2, lim;
+   long j, k, n = degpol(chi);
+   gpmem_t av2, lim;
    ns2 = manage_cache(chi, pp, ns);
-Line 1184  newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
+Line 1225  newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
 Line 1184  newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
 Line 1225  newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
    for (j = 1; j <= n; j++)
    {
      pa = gmul(pa, a);
-     if (pp) pa = polmodi(pa, pp);
+     pa = polmodi(pa, pp);
      pa = gmod(pa, chi);
-     if (pp) pa = polmodi(pa, pp);
+     pa = polmodi(pa, pp);
      s  = gzero;
      for (k = 0; k <= n-1; k++)
        s = addii(s, mulii(polcoeff0(pa, k, -1), (GEN)ns2[k+1]));
-     if (pp) va[j] = (long)centermod(s, pp);
+     va[j] = (long)centermod(s, pp);
      if (low_stack(lim, stack_lim(av2, 1)))
      {
        GEN *gptr[2];
-Line 1213  static GEN
+Line 1254  static GEN
 Line 1213  static GEN
 Line 1254  static GEN
  newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
  {
    GEN v, c, s, t;
-   long n = degpol(chi), j, k, vn = varn(chi), av = avma, av2, lim;
+   long n = degpol(chi), j, k, vn = varn(chi);
+   gpmem_t av = avma, av2, lim;
    v = newtonsums(a, chi, pp, ns);
    av2 = avma;
-Line 1240  newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
+Line 1282  newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
 Line 1240  newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
 Line 1282  newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
        c = gerepilecopy(av2, c);
      }
    }
    k = (n%2)? 1: 2;
    for (  ; k <= n+1; k += 2)
      c[k] = lneg((GEN)c[k]);
    return gerepileupto(av, gtopoly(c, vn));
  }
- static GEN
+ /* return v_p(n!) */
+ static long
+ val_fact(long n, GEN pp)
+ {
+   long p, q, v;
+   if (is_bigint(pp)) return 0;
+   q = p = itos(pp); v = 0;
+   do { v += n/q; q *= p; } while (n >= q);
+   return v;
+ }
+ static GEN
  mycaract(GEN f, GEN beta, GEN p, GEN pp, GEN ns)
  {
-   GEN p1, p2, chi, chi2, npp;
+   GEN p1, chi, npp;
-   long j, a, v = varn(f), n = degpol(f);
+   long v = varn(f), n = degpol(f);
    if (gcmp0(beta)) return zeropol(v);
-   p1 = content(beta);
+   beta = Q_primitive_part(beta,&p1);
-   if (gcmp1(p1)) p1 = NULL;
-   else beta = gdiv(beta, p1);
    if (!pp)
-     chi = caractducos(f, beta, v);
+     chi = ZX_caract(f, beta, v);
    else
    {
-     a = 0;
+     npp = mulii(pp, gpowgs(p, val_fact(n, p)));
-     for (j = 1; j <= n; j++) /* compute the extra precision needed */
+     if (p1) npp = mulii(npp, gpowgs(denom(p1), n));
-       a += ggval(stoi(j), p);
-     npp = mulii(pp, gpowgs(p, a));
-     if (p1) npp = gmul(npp, gpowgs(denom(p1), n));
      chi = newtoncharpoly(beta, f, npp, ns);
    }
-   if (p1)
+   if (p1) chi = rescale_pol(chi,p1);
-   {
-     chi2 = cgetg(n+3, t_POL);
-     chi2[1] = chi[1];
-     p2 = gun;
-     for (j = n+2; j >= 2; j--)
-     {
-       chi2[j] = lmul((GEN)chi[j], p2);
-       p2 = gmul(p2, p1);
-     }
-   }
-   else
-     chi2 = chi;
-   if (!pp) return chi2;
+   if (!pp) return chi;
    /* this can happen only if gamma is incorrect (not an integer) */
-   if (divise(denom(content(chi2)), p)) return NULL;
+   if (divise(Q_denom(chi), p)) return NULL;
-   return redelt(chi2, pp, pp);
+   return redelt(chi, pp, pp);
  }
  /* Factor characteristic polynomial of beta mod p */
  static GEN
  factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
  {
-   long av,l;
+   gpmem_t av;
    GEN chi,nu, b = cgetg(4,t_VEC);
+   long l;
    chi=mycaract(f,beta,p,pp,ns);
    av=avma; nu=(GEN)factmod(chi,p)[1]; l=lg(nu)-1;
-Line 1314  factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
+Line 1350  factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
 Line 1314  factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
 Line 1350  factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
  static GEN
  getprime(GEN p, GEN chi, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep)
  {
    long v = varn(chi), L, E, r, s;
-   GEN chin, pip, pp, vn;
+   GEN chin, pip, pp;
    if (gegal(nup, polx[v]))
      chin = chip;
    else
      chin = mycaract(chip, nup, p, NULL, NULL);
-   vn = vstar(p, chin);
-   L  = vn[0];
-   E  = vn[1];
-   cbezout(L, -E, &r, &s);
-   if (r <= 0)
+   vstar(p, chin, &L, &E);
-   {
+   (void)cbezout(L, -E, &r, &s);
-     long q = (-r) / E;
+   if (r <= 0)
-     q++;
+   {
-     r += q*E;
+     long q = 1 + ((-r) / E);
-     s += q*L;
+     r += q*E;
+     s += q*L;
    }
-   pip = eleval(chi, nup, phi);
+   pip = RX_RXQ_compo(nup, phi, chi);
-   pip = lift_intern(gpuigs(gmodulcp(pip, chi), r));
+   pip = lift_intern(gpowgs(gmodulcp(pip, chi), r));
-   pp  = gpuigs(p, s);
+   pp  = gpowgs(p, s);
    *Lp = L;
-   *Ep = E;
+   *Ep = E; return gdiv(pip, pp);
-   return gdiv(pip, pp);
  }
  static GEN
  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr, GEN pmf, long mf,
               GEN ns)
  {
    long l, v = varn(fx);
-Line 1366  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
+Line 1396  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1366  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1396  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
    for (;;)
    {
      if (signe(pdr)) break;
      if (!nalph) nalph = gadd(alph, gmul(p, polx[v]));
-     /* nchi was too reduced at this point */
+     /* nchi was too reduced at this point; try a larger precision */
-     nchi = mycaract(fx, nalph, NULL, NULL, ns);
+     pmf  = sqri(pmf);
-     pdr = respm(nchi, derivpol(nchi), pmf);
+     nchi = mycaract(fx, nalph, p, pmf, ns);
+     pdr  = respm(nchi, derivpol(nchi), pmf);
      if (signe(pdr)) break;
      if (DEBUGLEVEL >= 6)
        fprintferr("  non separable polynomial in update_alpha!\n");
      /* at this point, we assume that chi is not square-free */
      nalph = gadd(nalph, gmul(p, polx[v]));
      w = factcp(p, fx, nalph, NULL, ns);
      nchi = (GEN)w[1];
      nnu  = (GEN)w[2];
      l    = itos((GEN)w[3]);
      if (l > 1) return Decomp(p, fx, mf, nalph, nchi, nnu, 0);
      pdr = respm(nchi, derivpol(nchi), pmr);
-Line 1389  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
+Line 1420  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1389  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1420  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
    {
      npmr  = mulii(sqri(pdr), p);
      nchi  = polmodi(nchi, npmr);
-     if (!nalph) nalph = redelt(alph, npmr, pmf);
+     nalph = redelt(nalph? nalph: alph, npmr, pmf);
-     else nalph = redelt(nalph, npmr, pmf);
    }
    affii(gzero, (GEN)ns[1]); /* kill cache again (contains data for fx) */
-Line 1404  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
+Line 1434  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1404  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
 Line 1434  update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr
    return rep;
  }
- extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
+ /* flag != 0 iff we're looking for the p-adic factorization,
- /* flag != 0 iff we're looking for the p-adic factorization,
     in which case it is the p-adic precision we want */
  GEN
  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
  {
-   long Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
+   long L, E, Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
    GEN p1, alph, chi, nu, w, phi, pmf, pdr, pmr, kapp, pie, chib, ns;
-   GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, vb, opa;
+   GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, opa;
-   if (DEBUGLEVEL >= 3)
+   if (DEBUGLEVEL>2)
    {
-     if (flag)
+     fprintferr("  entering Nilord");
-       fprintferr("  entering Nilord2 (factorization)");
+     if (DEBUGLEVEL>4)
-     else
-       fprintferr("  entering Nilord2 (basis/discriminant)");
-     if (DEBUGLEVEL >= 5)
      {
        fprintferr(" with parameters: p = %Z, expo = %ld\n", p, mf);
        fprintferr("  fx = %Z, gx = %Z", fx, gx);
      }
      fprintferr("\n");
    }
    pmf = gpowgs(p, mf + 1);
    pdr = respm(fx, derivpol(fx), pmf);
    pmr = mulii(sqri(pdr), p);
-Line 1436  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1461  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1436  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1461  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
    chi = polmodi_keep(fx, pmr);
    alph = polx[v];
    nu = gx;
    N  = degpol(fx);
    oE = 0;
    opa = NULL;
-Line 1452  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1477  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1452  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1477  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
    for(;;)
    {
-     /* kappa need to be recomputed */
+     /* kappa needs to be recomputed */
      kapp = NULL;
      Fa   = degpol(nu);
      /* the prime element in Zp[alpha] */
-Line 1471  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1496  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1471  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1496  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
        pia  = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
        pia  = redelt(pia, pmr, pmf);
      }
-     oE = Ea; opa = eleval(fx, pia, alph);
+     oE = Ea; opa = RX_RXQ_compo(pia, alph, fx);
      if (DEBUGLEVEL >= 5)
        fprintferr("  Fa = %ld and Ea = %ld \n", Fa, Ea);
      /* we change alpha such that nu = pia */
      if (La > 1)
      {
-       alph = gadd(alph, eleval(fx, pia, alph));
+       alph = gadd(alph, RX_RXQ_compo(pia, alph, fx));
        w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
        alph = (GEN)w[1];
-Line 1490  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1515  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1490  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1515  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
      }
      /* if Ea*Fa == N then O = Zp[alpha] */
      if (Ea*Fa == N)
      {
        if (flag) return NULL;
-Line 1515  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1540  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1515  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1540  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          eq = (long)(vn / N);
          er = (long)(vn*Ea/N - eq*Ea);
        }
        else
        {
          chib = mycaract(chi, beta, NULL, NULL, ns);
-         vb = vstar(p, chib);
+         vstar(p, chib, &L, &E);
-         eq = (long)(vb[0] / vb[1]);
+         eq = (long)(L / E);
-         er = (long)(vb[0]*Ea / vb[1] - eq*Ea);
+         er = (long)(L*Ea / E - eq*Ea);
        }
        /* eq and er are such that gamma = beta.p^-eq.nu^-er is a unit */
        if (eq) gamm = gdiv(beta, gpowgs(p, eq));
        else gamm = beta;
-Line 1532  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1557  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1532  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1557  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          /* kappa = nu^-1 in Zp[alpha] */
          if (!kapp)
          {
-           kapp = ginvmod(nu, chi);
+           kapp = QX_invmod(nu, chi);
            kapp = redelt(kapp, pmr, pmf);
            kapp = gmodulcp(kapp, chi);
          }
-Line 1544  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1569  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1544  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1569  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          fprintferr("  gamma = %Z\n", gamm);
        if (er || !chib)
          chig = mycaract(chi, gamm, p, pmf, ns);
        else
        {
          chig = poleval(chib, gmul(polx[v], gpowgs(p, eq)));
-Line 1552  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1577  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1552  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1577  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          chig = polmodi(chig, pmf);
        }
-       if (!chig || !gcmp1(denom(content(chig))))
+       if (!chig || !gcmp1(Q_denom(chig)))
        {
-         /* the valuation of beta was wrong... This also means
+         /* Valuation of beta was wrong. This means that
-            that chi_gamma has more than one factor modulo p   */
+            gamma fails the v*-test */
-         if (!chig) chig = mycaract(chi, gamm, p, NULL, NULL);
+         chib = mycaract(chi, beta, p, NULL, ns);
+         vstar(p, chib, &L, &E);
+         eq = (long)(-L / E);
+         er = (long)(-L*Ea / E - eq*Ea);
-         vb = vstar(p, chig);
+         if (eq) gamm = gmul(beta, gpowgs(p, eq));
-         eq = (long)(-vb[0] / vb[1]);
+         else gamm = beta;
-         er = (long)(-vb[0]*Ea / vb[1] - eq*Ea);
+         if (er)
-         if (eq) gamm = gmul(gamm, gpowgs(p, eq));
-         if (er)
          {
            gamm = gmul(gamm, gpowgs(nu, er));
            gamm = gmod(gamm, chi);
            gamm = redelt(gamm, p, pmr);
          }
-         if (eq || er) chig = mycaract(chi, gamm, p, pmf, ns);
+         chig = mycaract(chi, gamm, p, pmf, ns);
        }
        nug  = (GEN)factmod(chig, p)[1];
-Line 1578  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1604  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1578  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1604  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
        if (l > 1)
        {
          /* there are at least 2 factors mod. p => chi can be split */
-         phi  = eleval(fx, gamm, alph);
+         phi  = RX_RXQ_compo(gamm, alph, fx);
          phi  = redelt(phi, p, pmf);
          if (flag) mf += 3;
          return Decomp(p, fx, mf, phi, chig, nug, flag);
-Line 1594  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1620  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1594  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1620  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          if (gcmp1((GEN)w[1]))
          {
            /* there are at least 2 factors mod. p => chi can be split */
-           phi = eleval(fx, (GEN)w[2], alph);
+           phi = RX_RXQ_compo((GEN)w[2], alph, fx);
            phi = redelt(phi, p, pmf);
            if (flag) mf += 3;
            return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
-Line 1613  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1639  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1613  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1639  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
        for (i = 1;; i++)
        {
          if (i >= lg(w))
            err(talker, "bug in nilord (no suitable root), is p = %Z a prime?",
                (long)p);
          delt = gneg_i(gsubst(gcoeff(w, 2, i), nv, polx[v]));
          eta  = gsub(gamm, delt);
          if (typ(delt) == t_INT)
          {
            chie = poleval(chig, gadd(polx[v], delt));
-Line 1625  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1651  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1625  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1651  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
            nue  = lift((GEN)nue[l]);
          }
          else
          {
            p1   = factcp(p, chi, eta, pmf, ns);
            chie = (GEN)p1[1];
            nue  = (GEN)p1[2];
            l    = itos((GEN)p1[3]);
          }
          if (l > 1)
          {
            /* there are at least 2 factors mod. p => chi can be split */
-           delete_var();
+           (void)delete_var();
-           phi = eleval(fx, eta, alph);
+           phi = RX_RXQ_compo(eta, alph, fx);
            phi = redelt(phi, p, pmf);
            if (flag) mf += 3;
            return Decomp(p, fx, mf, phi, chie, nue, flag);
          }
          /* if vp(eta) = vp(gamma - delta) > 0 */
          if (gegal(nue, polx[v])) break;
        }
-       delete_var();
+       (void)delete_var();
-       if (!signe(modii((GEN)chie[2], pmr)))
+       if (!signe(modii(constant_term(chie), pmr)))
          chie = mycaract(chi, eta, p, pmf, ns);
        pie = getprime(p, chi, eta, chie, nue, &Le, &Ee);
-Line 1660  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1686  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1660  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1686  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
          if (gcmp1((GEN)w[1]))
          {
            /* there are at least 2 factors mod. p => chi can be split */
-           phi = eleval(fx, (GEN)w[2], alph);
+           phi = RX_RXQ_compo((GEN)w[2], alph, fx);
            phi = redelt(phi, p, pmf);
            if (flag) mf += 3;
            return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
          }
          break;
        }
        if (eq) delt = gmul(delt, gpowgs(p, eq));
        if (er) delt = gmul(delt, gpowgs(nu, er));
        beta = gsub(beta, delt);
      }
      /* we replace alpha by a new alpha with a larger F or E */
-     alph = eleval(fx, (GEN)w[2], alph);
+     alph = RX_RXQ_compo((GEN)w[2], alph, fx);
      chi  = (GEN)w[3];
      nu   = (GEN)w[4];
-Line 1702  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
+Line 1728  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1702  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
 Line 1728  nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
  }
  /* Returns [1,phi,chi,nu] if phi non-primary
   *         [2,phi,chi,nu] with F_phi = lcm (F_alpha, F_theta)
   *                         and E_phi = E_alpha
   */
  static GEN
  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, long Ft, GEN ns)
-Line 1711  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
+Line 1737  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1711  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1737  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
    long m, Dat, t, v = varn(fa);
    GEN b, w, phi, h;
    Dat = clcm(Fa, Ft) + 3;
    b = cgetg(5, t_VEC);
    m = p[2];
    if (degpol(p) > 0 || m < 0) m = 0;
    for (t = 1;; t++)
-Line 1725  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
+Line 1751  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1725  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1751  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
      if (h[2] > 1) { b[1] = un; break; }
      if (lgef(w[2]) == Dat) { b[1] = deux; break; }
    }
    b[2] = (long)phi;
    b[3] = w[1];
    b[4] = w[2];
    return b;
  }
-Line 1736  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
+Line 1762  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1736  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
 Line 1762  testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, lon
   *         [2, phi, chi, nu] if E_phi = lcm (E_alpha, E_theta)
   */
  static GEN
  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, long Ea, GEN thet2,
         long Et, GEN ns)
  {
    GEN b, c1, c2, c3, psi, phi, w, h;
-Line 1744  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
+Line 1770  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
 Line 1744  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
 Line 1770  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
    b=cgetg(5, t_VEC);
-   cbezout(Ea, Et, &r, &s); t = 0;
+   (void)cbezout(Ea, Et, &r, &s); t = 0;
    while (r < 0) { r = r + Et; t++; }
    while (s < 0) { s = s + Ea; t++; }
-   c1 = lift_intern(gpuigs(gmodulcp(alph2, fa), s));
+   c1 = lift_intern(gpowgs(gmodulcp(alph2, fa), s));
-   c2 = lift_intern(gpuigs(gmodulcp(thet2, fa), r));
+   c2 = lift_intern(gpowgs(gmodulcp(thet2, fa), r));
-   c3 = gdiv(gmod(gmul(c1, c2), fa), gpuigs(p, t));
+   c3 = gdiv(gmod(gmul(c1, c2), fa), gpowgs(p, t));
    psi = redelt(c3, pmr, pmr);
    phi = gadd(polx[v], psi);
-Line 1758  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
+Line 1784  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
 Line 1758  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
 Line 1784  testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, lon
    w = factcp(p, fa, phi, pmf, ns);
    h = (GEN)w[3];
    b[1] = (h[2] > 1)? un: deux;
    b[2] = (long)phi;
    b[3] = w[1];
    b[4] = w[2];
    return b;
  }
- /* evaluate h(a) mod f */
+ /*******************************************************************/
- GEN
+ /*                                                                 */
- eleval(GEN f,GEN h,GEN a)
+ /*    2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index)       */
+ /*                                                                 */
+ /*******************************************************************/
+ /* beta = generators of prime ideal (vectors). Return "small" random elt in P */
+ static GEN
+ random_elt_in_P(GEN beta, long small)
  {
-   long n,av,tetpil;
+   long z, i, j, la, lbeta = lg(beta);
-   GEN y;
+   GEN a = NULL, B;
-   if (typ(h) != t_POL) return gcopy(h);
+   /* compute sum random(-7..8) * beta[i] */
-   av = tetpil = avma;
+   if (small)
-   n=lgef(h)-1; y=(GEN)h[n];
-   for (n--; n>=2; n--)
    {
-     y = gadd(gmul(y,a),(GEN)h[n]);
+     la = lg(beta[1]);
-     tetpil=avma; y = gmod(y,f);
+     if (small > 7) small = 0;
+     for (i=1; i<lbeta; i++)
+     { /* Warning: change definition of 'small' if you change this */
+       z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
+       if (z >= 9) z -= 7;
+       if (small) z %= small;
+       if (!z) continue;
+       B = (GEN)beta[i];
+       if (a)
+         for (j = 1; j < la; j++) a[j] += z * B[j];
+       else {
+         a = cgetg(la, t_VECSMALL);
+         for (j = 1; j <la; j++) a[j] = z * B[j];
+       }
+     }
    }
-   return gerepile(av,tetpil,y);
+   else
+     for (i=1; i<lbeta; i++)
+     {
+       z = mymyrand() >> (BITS_IN_RANDOM-5);
+       if (z >= 9) z -= 7;
+       if (!z) continue;
+       B = gmulsg(z, (GEN)beta[i]);
+       a = a? gadd(a, B): B;
+     }
+   return a;
  }
- GEN addshiftw(GEN x, GEN y, long d);
+ /* routines to check uniformizer algebraically (as t_POL) */
- static GEN
+ /* a t_POL, D t_INT (or NULL if 1), q = p^(f+1). a/D in pr | p, norm(pr) = pf.
- shiftpol(GEN x, long v)
+  * Return 1 if (a/D,p) = pr, and 0 otherwise */
+ static int
+ is_uniformizer(GEN D, GEN a, GEN T, GEN q)
  {
-   x = addshiftw(x, zeropol(v), 1);
+   GEN N = ZX_resultant_all(T, a, D, 0); /* norm(a) */
-   setvarn(x,v); return x;
+   return (resii(N, q) != gzero);
  }
- /* Sylvester's matrix, mod p^m (assumes f1 monic) */
+ /* a/D or a/D + p uniformizer ? */
  static GEN
- sylpm(GEN f1,GEN f2,GEN pm)
+ prime_check_elt(GEN D, GEN Dp, GEN a, GEN T, GEN q, int ramif)
  {
-   long n,j,v=varn(f1);
+   if (is_uniformizer(D,a,T,q)) return a;
-   GEN a,h;
+   /* can't succeed if e > 1 */
+   if (!ramif && is_uniformizer(D,gadd(a,Dp),T,q)) return a;
+   return NULL;
+ }
-   n=degpol(f1); a=cgetg(n+1,t_MAT);
+ /* P given by an Fp basis */
-   h = FpX_res(f2,f1,pm);
+ static GEN
-   for (j=1;; j++)
+ compute_pr2(GEN nf, GEN P, GEN p, int *ramif)
+ {
+   long i, j, m = lg(P)-1;
+   GEN mul, P2 = gmul(p, P);
+   for (i=1; i<=m; i++)
    {
-     a[j] = (long)pol_to_vec(h,n);
+     mul = eltmul_get_table(nf, (GEN)P[i]);
-     if (j == n) break;
+     for (j=1; j<=i; j++)
-     h = FpX_res(shiftpol(h,v),f1,pm);
+       P2 = concatsp(P2, gmul(mul, (GEN)P[j]));
    }
-   return hnfmodid(a,pm);
+   P2 = hnfmodid(P2, sqri(p)); /* HNF of pr^2 */
+   *ramif = egalii(gcoeff(P2,1,1), p); /* pr/p ramified ? */
+   return P2;
  }
- /* polynomial gcd mod p^m (assumes f1 monic) */
+ static GEN
- GEN
+ random_unif_loop_pol(GEN nf, GEN P, GEN D, GEN Dp, GEN beta, GEN pol,
- gcdpm(GEN f1,GEN f2,GEN pm)
+                      GEN p, GEN q)
  {
-   long n,c,v=varn(f1),av=avma,tetpil;
+   long small, i, c = 0, m = lg(P)-1, N = degpol(nf[1]), keep = getrand();
-   GEN a,col;
+   int ramif;
+   GEN a, P2;
+   gpmem_t av;
-   n=degpol(f1); a=sylpm(f1,f2,pm);
+   for(i=1; i<=m; i++)
-   for (c=1; c<=n; c++)
+     if ((a = prime_check_elt(D,Dp,(GEN)beta[i],pol,q,0))) return a;
-     if (signe(resii(gcoeff(a,c,c),pm))) break;
-   if (c > n) { avma=av; return zeropol(v); }
-   col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
+   (void)setrand(1);
-   return gerepile(av,tetpil, gtopolyrev(col,v));
+   if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
+   P2 = compute_pr2(nf, P, p, &ramif);
+   a = mulis(Dp, 8*degpol(nf[1]));
+   if (is_bigint(a)) small = 0;
+   else
+   {
+     ulong mod = itou(Dp);
+     small = itos(p);
+     for (i=1; i<=m; i++)
+       beta[i] = (long)u_Fp_FpV(pol_to_vec((GEN)beta[i], N), mod);
+   }
+   for(av = avma; ;avma = av)
+   {
+     if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
+     a = random_elt_in_P(beta, small);
+     if (!a) continue;
+     if (small) a = vec_to_pol(small_to_vec(a), varn(pol));
+     a = centermod(a, Dp);
+     if ((a = prime_check_elt(D,Dp,a,pol,q,ramif)))
+     {
+       if (DEBUGLEVEL) fprintferr("\n");
+       (void)setrand(keep); return a;
+     }
+   }
  }
- /* reduced resultant mod p^m (assumes x monic) */
+ /* routines to check uniformizer numerically (as t_COL) */
- GEN
- respm(GEN x,GEN y,GEN pm)
- {
-   long av = avma;
-   GEN p1 = sylpm(x,y,pm);
-   p1 = gcoeff(p1,1,1);
+ static int
-   if (egalii(p1,pm)) { avma = av; return gzero; }
+ vec_is_uniformizer(GEN a, GEN q, long r1)
-   return gerepileuptoint(av, icopy(p1));
+ {
+   long e;
+   GEN N = grndtoi( norm_by_embed(r1, a), &e );
+   if (e > -5) err(precer, "vec_is_uniformizer");
+   return (resii(N, q) != gzero);
  }
- /* Normalized integral basis */
  static GEN
- nbasis(GEN ibas,GEN pd)
+ prime_check_eltvec(GEN a, GEN M, GEN p, GEN q, long r1, long small, int ramif)
  {
-   long k, n = lg(ibas)-1;
+   GEN ar;
-   GEN a = cgetg(n+1,t_MAT);
+   long i,l;
-   for (k=1; k<=n; k++) a[k] = (long)pol_to_vec((GEN)ibas[k],n);
+   a = centermod(a, p);
-   return gdiv(hnfmodid(a,pd), pd);
+   if (small) ar = gmul_mat_smallvec(M, a);
+   else       ar = gmul(M, a);
+   if (vec_is_uniformizer(ar, q, r1)) return small? small_to_col(a): a;
+   /* can't succeed if e > 1 */
+   if (ramif) return NULL;
+   l = lg(ar);
+   for (i=1; i<l; i++) ar[i] = ladd((GEN)ar[i], p);
+   if (!vec_is_uniformizer(ar, q, r1)) return NULL;
+   if (small) a = small_to_col(a);
+   a[1] = laddii((GEN)a[1],p); return a;
  }
- /*******************************************************************/
- /*                                                                 */
- /*                   BUCHMANN-LENSTRA ALGORITHM                    */
- /*                                                                 */
- /*******************************************************************/
- static GEN lens(GEN nf,GEN p,GEN a);
- GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
- /* return a Z basis of Z_K's p-radical, modfrob = x--> x^p-x */
  static GEN
- pradical(GEN nf, GEN p, GEN *modfrob)
+ random_unif_loop_vec(GEN nf, GEN P, GEN p, GEN q)
  {
-   long i,N = degpol(nf[1]);
+   long small, r1, i, c = 0, m = lg(P)-1, keep = getrand();
-   GEN p1,m,frob,rad;
+   int ramif;
+   GEN a, P2, beta, M = gmael(nf,5,1);
+   gpmem_t av;
-   frob = cgetg(N+1,t_MAT);
+   r1 = nf_get_r1(nf);
-   for (i=1; i<=N; i++)
+   a = mulis(p, 8*degpol(nf[1]));
-     frob[i] = (long) element_powid_mod_p(nf,i,p, p);
+   if (is_bigint(a)) { small = 0; beta = P; }
-   /* p1 = smallest power of p st p^k >= N */
-   p1=p; while (cmpis(p1,N)<0) p1=mulii(p1,p);
-   if (p1==p) m = frob;
    else
    {
-     m=cgetg(N+1,t_MAT); p1 = divii(p1,p);
+     small = itos(p); beta = cgetg(m+1, t_MAT);
-     for (i=1; i<=N; i++)
+     for (i=1; i<=m; i++)
-       m[i]=(long)element_pow_mod_p(nf,(GEN)frob[i],p1, p);
+       beta[i] = (long)u_Fp_FpV((GEN)P[i], itou(p));
    }
-   rad = FpM_ker(m, p);
+   for(i=1; i<=m; i++)
-   for (i=1; i<=N; i++)
+     if ((a = prime_check_eltvec((GEN)beta[i],M,p,q,r1,small,0))) return a;
-     coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
-   *modfrob = frob; return rad;
- }
- static GEN
+   (void)setrand(1);
- project(GEN algebre, GEN x, long k, long kbar)
+   if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
- {
+   P2 = compute_pr2(nf, P, p, &ramif);
-   x = inverseimage(algebre,x);
-   x += k; x[0] = evaltyp(t_COL) | evallg(kbar+1);
-   return x;
- }
- /* Calcule le polynome minimal de alpha dans algebre (coeffs dans Z) */
+   for(av = avma; ;avma = av)
- static GEN
- pol_min(GEN alpha,GEN nf,GEN algebre,long kbar,GEN p)
- {
-   long av=avma,tetpil,i,N,k;
-   GEN p1,puiss;
-   N = lg(nf[1])-3; puiss=cgetg(N+2,t_MAT);
-   k = N-kbar; p1=alpha;
-   puiss[1] = (long)gscalcol_i(gun,kbar);
-   for (i=2; i<=N+1; i++)
    {
-     if (i>2) p1 = element_mul(nf,p1,alpha);
+     if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
-     puiss[i] = (long) project(algebre,p1,k,kbar);
+     a = random_elt_in_P(beta, small);
+     if (!a) continue;
+     if ((a = prime_check_eltvec(a,M,p,q,r1,small,0)))
+     {
+       if (DEBUGLEVEL) fprintferr("\n");
+       (void)setrand(keep); return a;
+     }
    }
-   puiss = lift_intern(puiss);
-   p1 = (GEN)FpM_ker(puiss, p)[1]; tetpil=avma;
-   return gerepile(av,tetpil,gtopolyrev(p1,0));
  }
- /* Evalue le polynome pol en alpha,element de nf */
+ /* L = Fp-bases for prime ideals of O_K/p. Return a p-uniformizer for
+  * L[i] using the approximation theorem */
  static GEN
- eval_pol(GEN nf,GEN pol,GEN alpha,GEN algebre,GEN algebre1)
+ uniformizer_appr(GEN nf, GEN L, long i, GEN p)
  {
-   long av=avma,tetpil,i,kbar,k, lx = lgef(pol)-1, N = degpol(nf[1]);
+   long j, l;
-   GEN res;
+   GEN inter = NULL, P = (GEN)L[i], P2, Q, u, v, pi;
+   int ramif;
-   kbar = lg(algebre1)-1; k = N-kbar;
+   l = lg(L)-1;
-   res = gscalcol_i((GEN)pol[lx], N);
+   for (j=l; j; j--)
-   for (i=2; i<lx; i++)
    {
-     res = element_mul(nf,alpha,res);
+     if (j == i) continue;
-     res[1] = ladd((GEN)res[1],(GEN)pol[i]);
+     Q = (GEN)L[j]; if (typ(Q) != t_MAT) break;
+     inter = inter? FpM_intersect(inter, Q, p): Q;
    }
-   res = project(algebre,res,k,kbar); tetpil=avma;
+   if (!inter)
-   return gerepile(av,tetpil,gmul(algebre1,res));
+   {
+     setlg(L, l);
+     inter = idealprodprime(nf, L);
+   }
+   else
+   {
+     inter = hnfmodid(inter, p);
+     for (   ; j; j--)
+     {
+       Q = (GEN)L[j];
+       inter = idealmul(nf, inter, Q);
+     }
+   }
+   /* inter = prod L[j], j != i */
+   P2 = compute_pr2(nf, P, p, &ramif);
+   u = addone_aux2(P2, inter);
+   v = unnf_minus_x(u);
+   if (!ramif) pi = gmul(p, v);
+   else
+   {
+     l = lg(P)-1;
+     for (j=l; j; j--)
+       if (!hnf_invimage(P2, (GEN)P[j])) break;
+     pi = element_muli(nf,(GEN)P[j],v);
+   }
+   pi = gadd(pi, u);
+   pi =  centermod(pi, p);
+   if (!ramif && hnf_invimage(P2, pi)) pi[1] = laddii((GEN)pi[1], p);
+   return pi;
  }
+ /* Input: an ideal mod p, P != Z_K (Fp-basis, in matrix form)
+  * Output: x such that P = (p,x) */
  static GEN
- kerlens2(GEN x, GEN p)
+ uniformizer(GEN nf, GEN P, GEN p)
  {
-   long i,j,k,t,nbc,nbl,av,av1;
+   GEN q, D, Dp, w, beta, a, T = (GEN)nf[1];
-   GEN a,c,l,d,y,q;
+   long i, f, N = degpol(T), m = lg(P)-1;
-   av=avma; a=gmul(x,gmodulsg(1,p));
+   if (!m) return gscalcol_i(p,N);
-   nbl=nbc=lg(x)-1;
-   c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
+   /* we want v_p(Norm(beta)) = p^f, f = N-m */
-   l=new_chunk(nbc+1);
+   f = N-m;
-   d=new_chunk(nbc+1);
+   P = centermod(P, p);
-   k = t = 1;
+   q = mulii(gpowgs(p,f), p);
-   while (t<=nbl && k<=nbc)
+   if (typ(nf[5]) == t_VEC) /* dummy nf from padicff */
    {
-     for (j=1; j<k; j++)
+     GEN M = gmael(nf,5,1);
-       for (i=1; i<=nbl; i++)
+     long ex, prec = gprecision(M);
-         if (i!=l[j])
+     ex = gexpo(M) + gexpo(mulsi(8 * N, p));
-           coeff(a,i,k) = lsub(gmul((GEN)d[j],gcoeff(a,i,k)),
+     /* enough prec to use norm_by_embed */
-                               gmul(gcoeff(a,l[j],k),gcoeff(a,i,j)));
+     if (N * ex <= bit_accuracy(prec))
-     t=1; while (t<=nbl && (c[t] || gcmp0(gcoeff(a,t,k)))) t++;
+       return random_unif_loop_vec(nf, P, p, q);
-     if (t<=nbl) { d[k]=coeff(a,t,k); c[t]=k; l[k]=t; k++; }
    }
-   if (k>nbc) err(bugparier,"kerlens2");
-   y=cgetg(nbc+1,t_COL);
+   w = Q_remove_denom((GEN)nf[7], &D);
-   y[1]=(k>1)?coeff(a,l[1],k):un;
+   if (D)
-   for (q=gun,j=2; j<k; j++)
    {
-     q=gmul(q,(GEN)d[j-1]);
+     long v = pvaluation(D, p, &a);
-     y[j]=lmul(gcoeff(a,l[j],k),q);
+     D = gpowgs(p, v);
+     Dp = mulii(D, p);
+   } else {
+     w = dummycopy(w);
+     Dp = p;
    }
-   if (k>1) y[k]=lneg(gmul(q,(GEN)d[k-1]));
+   for (i=1; i<=N; i++)
-   for (j=k+1; j<=nbc; j++) y[j]=zero;
+     w[i] = (long)FpX_red((GEN)w[i], Dp); /* w[i] / D defined mod p */
-   av1=avma; return gerepile(av,av1,lift(y));
+   beta = gmul(w, P);
+   a = random_unif_loop_pol(nf,P,D,Dp,beta,T,p,q);
+   a = algtobasis_i(nf,a);
+   if (D) a = gdivexact(a,D);
+   a = centermod(a, p);
+   if (!is_uniformizer(D,gmul(w,a), T,q)) a[1] = laddii((GEN)a[1],p);
+   return a;
  }
+ /* Assuming P = (p,u) prime, return tau such that p Z + tau Z = p P^(-1)*/
  static GEN
- kerlens(GEN x, GEN pgen)
+ anti_uniformizer(GEN nf, GEN p, GEN u)
  {
-   long av = avma, i,j,k,t,nbc,nbl,p,q,*c,*l,*d,**a;
+   gpmem_t av = avma;
-   GEN y;
+   GEN mat = eltmul_get_table(nf, u);
+   return gerepileupto(av, FpM_deplin(mat,p));
+ }
-   if (cmpis(pgen, MAXHALFULONG>>1) > 0)
+ /*******************************************************************/
-     return kerlens2(x,pgen);
+ /*                                                                 */
-   /* ici p <= (MAXHALFULONG>>1) ==> long du C */
+ /*                   BUCHMANN-LENSTRA ALGORITHM                    */
-   p=itos(pgen); nbl=nbc=lg(x)-1;
+ /*                                                                 */
-   a=(long**)new_chunk(nbc+1);
+ /*******************************************************************/
-   for (j=1; j<=nbc; j++)
+ /* pr = (p,u) of ramification index e */
+ GEN
+ apply_kummer(GEN nf,GEN u,GEN e,GEN p)
+ {
+   GEN t, T = (GEN)nf[1], pr = cgetg(6,t_VEC);
+   long f = degpol(u), N = degpol(T);
+   pr[1] = (long)p;
+   pr[3] = (long)e;
+   pr[4] = lstoi(f);
+   if (f == N) /* inert */
    {
-     c=a[j]=new_chunk(nbl+1);
+     pr[2] = (long)gscalcol_i(p,N);
-     for (i=1; i<=nbl; i++) c[i]=smodis(gcoeff(x,i,j),p);
+     pr[5] = (long)gscalcol_i(gun,N);
    }
-   c=new_chunk(nbl+1); for (i=1; i<=nbl; i++) c[i]=0;
+   else
-   l=new_chunk(nbc+1);
+   { /* make sure v_pr(u) = 1 (automatic if e>1) */
-   d=new_chunk(nbc+1);
+     if (is_pm1(e) && !is_uniformizer(NULL,u, T, gpowgs(p,f+1)))
-   k = t = 1;
+       u[2] = laddii((GEN)u[2], p);
-   while (t<=nbl && k<=nbc)
+     t = algtobasis_i(nf, FpX_div(T,u,p));
-   {
+     pr[2] = (long)algtobasis_i(nf,u);
-     for (j=1; j<k; j++)
+     pr[5] = (long)centermod(t, p);
-       for (i=1; i<=nbl; i++)
-         if (i!=l[j])
-           a[k][i] = (d[j]*a[k][i] - a[j][i]*a[k][l[j]]) % p;
-     t=1; while (t<=nbl && (c[t] || !a[k][t])) t++;
-     if (t<=nbl) { d[k]=a[k][t]; c[t]=k; l[k++]=t; }
    }
-   if (k>nbc) err(bugparier,"kerlens");
+   return pr;
-   avma=av; y=cgetg(nbc+1,t_COL);
-   t=(k>1) ? a[k][l[1]]:1;
-   y[1]=(t>0)? lstoi(t):lstoi(t+p);
-   for (q=1,j=2; j<k; j++)
-   {
-     q = (q*d[j-1]) % p;
-     t = (a[k][l[j]]*q) % p;
-     y[j] = (t>0) ? lstoi(t) : lstoi(t+p);
-   }
-   if (k>1)
-   {
-     t = (q*d[k-1]) % p;
-     y[k] = (t>0) ? lstoi(p-t) : lstoi(-t);
-   }
-   for (j=k+1; j<=nbc; j++) y[j]=zero;
-   return y;
  }
- /* Calcule la constante de lenstra de l'ideal p.Z_K+a.Z_K ou a est un
+ /* return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
- vecteur sur la base d'entiers */
  static GEN
- lens(GEN nf, GEN p, GEN a)
+ pradical(GEN nf, GEN p, GEN *phi)
  {
-   long av=avma,tetpil,N=degpol(nf[1]),j;
+   long i,N = degpol(nf[1]);
-   GEN mat=cgetg(N+1,t_MAT);
+   GEN q,m,frob,rad;
-   for (j=1; j<=N; j++) mat[j]=(long)element_mulid(nf,a,j);
-   tetpil=avma; return gerepile(av,tetpil,kerlens(mat,p));
- }
- extern GEN det_mod_P_n(GEN a, GEN N, GEN P);
+   /* matrix of Frob: x->x^p over Z_K/p */
- GEN sylvestermatrix_i(GEN x, GEN y);
+   frob = cgetg(N+1,t_MAT);
+   for (i=1; i<=N; i++)
+     frob[i] = (long)element_powid_mod_p(nf,i,p, p);
- /* check if p^va doesnt divide norm x (or norm(x+p)) */
+   m = frob; q = p;
- #if 0
+   while (cmpis(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
- /* compute norm mod p^whatneeded using Sylvester's matrix */
+   rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
- /* looks slower than the new subresultant. Have to re-check this */
+   for (i=1; i<=N; i++)
+     coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
+   *phi = frob; return rad;
+ }
+ /* return powers of a: a^0, ... , a^d,  d = dim A */
  static GEN
- prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
+ get_powers(GEN mul, GEN p)
  {
-   GEN M,mod,x, c = denom(content(a));
+   long i, d = lg(mul[1]);
-   long v = pvaluation(c, p, &x); /* x is junk */
+   GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
-   mod = mulii(pf, gpowgs(p, degpol(pol)*v + 1));
+   P[0] = (long)gscalcol_i(gun, d-1);
+   z = (GEN)mul[1];
-   x = FpX_red(gmul(a,c), mod);
+   for (i=1; i<=d; i++)
-   M = sylvestermatrix_i(pol,x);
-   if (det_mod_P_n(M,mod,p) == gzero)
    {
-     x[2] = ladd((GEN)x[2], mulii(p,c));
+     P[i] = (long)z; /* a^i */
-     M = sylvestermatrix_i(pol,x);
+     if (i!=d) z = FpM_FpV_mul(mul, z, p);
-     if (det_mod_P_n(M,mod,p) == gzero) return NULL;
-     a[2] = ladd((GEN)a[2], p);
    }
-   return a;
+   return pow;
  }
- #else
- /* use subres to compute norm */
+ /* minimal polynomial of a in A (dim A = d).
+  * mul = multiplication table by a in A */
  static GEN
- prime_check_elt(GEN a, GEN pol, GEN p, GEN pf)
+ pol_min(GEN mul, GEN p)
  {
-   GEN norme=subres(pol,a);
+   gpmem_t av = avma;
-   if (resii(divii(norme,pf),p) != gzero) return a;
+   GEN z, pow = get_powers(mul, p);
-   /* Note: a+p can't succeed if e > 1, can we know this at this point ? */
+   z = FpM_deplin(pow, p);
-   a=gadd(a,p); norme=subres(pol,a);
+   return gerepileupto(av, gtopolyrev(z,0));
-   if (resii(divii(norme,pf),p) != gzero) return a;
-   return NULL;
  }
- #endif
- #if 0
  static GEN
- prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
+ get_pr(GEN nf, GEN L, long i, GEN p, int appr)
  {
-   long av, m = lg(beta)-1;
+   GEN pr, u, t, P = (GEN)L[i];
-   int i,j,K, *x = (int*)new_chunk(m+1);
+   long e, f;
-   GEN a;
+   gpmem_t av;
-   K = 1; av = avma;
+   if (typ(P) == t_VEC) return P; /* already done (Kummer) */
-   for(;;)
-   { /* x runs through strictly increasing sequences of length K,
-      * 1 <= x[i] <= m */
- nextK:
-     if (DEBUGLEVEL) fprintferr("K = %d\n", K);
-     for (i=1; i<=K; i++) x[i] = i;
-     for(;;)
-     {
-       if (DEBUGLEVEL > 1)
-       {
-         for (i=1; i<=K; i++) fprintferr("%d ",x[i]);
-         fprintferr("\n"); flusherr();
-       }
-       a = (GEN)beta[x[1]];
-       for (i=2; i<=K; i++) a = gadd(a, (GEN)beta[x[i]]);
-       if ((a = prime_check_elt(a,pol,p,pf))) return a;
-       avma = av;
-       /* start: i = K+1; */
+   av = avma;
-       do
+   if (appr) u = uniformizer_appr(nf, L, i, p);
-       {
+   else      u = uniformizer(nf, P, p);
-         if (--i == 0)
+   u = gerepilecopy(av, u);
-         {
+   t = anti_uniformizer(nf,p,u);
-           if (++K > m) return NULL; /* fail */
+   av = avma;
-           goto nextK;
+   e = 1 + int_elt_val(nf,t,p,t,NULL);
-         }
+   f = degpol(nf[1]) - (lg(P)-1);
-         x[i]++;
+   avma = av;
-       } while (x[i] > m - K + i);
+   pr = cgetg(6,t_VEC);
-       for (j=i; j<K; j++) x[j+1] = x[j]+1;
+   pr[1] = (long)p;
-     }
+   pr[2] = (long)u;
-   }
+   pr[3] = lstoi(e);
+   pr[4] = lstoi(f);
+   pr[5] = (long)t; return pr;
  }
- #endif
- static GEN
+ static int
- random_prime_two_elt_loop(GEN beta, GEN pol, GEN p, GEN pf)
+ use_appr(GEN L, GEN pp, long N)
  {
-   long av = avma, z,i, m = lg(beta)-1;
+   long i, f, l = lg(L);
-   long keep = getrand();
+   double prod = 1., NP;
-   int c = 0;
+   GEN P;
-   GEN a;
+   gpmem_t av = avma;
-   for(i=1; i<=m; i++)
+   for (i=1; i<l; i++)
-     if ((a = prime_check_elt((GEN)beta[i],pol,p,pf))) return a;
-   (void)setrand(1);
-   if (DEBUGLEVEL) fprintferr("prime_two_elt_loop, hard case: ");
-   for(;;avma=av)
    {
-     if (DEBUGLEVEL) fprintferr("%d ", ++c);
+     P = (GEN)L[i];
-     a = gzero;
+     f = typ(P)==t_VEC? itos((GEN)P[4]): N - (lg(P)-1);
-     for (i=1; i<=m; i++)
+     NP = gtodouble(gpowgs(pp, f));
-     {
+     prod *= (1 - 1./NP);
-       z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
-       if (z >= 9) z -= 7;
-       a = gadd(a,gmulsg(z,(GEN)beta[i]));
-     }
-     if ((a = prime_check_elt(a,pol,p,pf)))
-     {
-       if (DEBUGLEVEL) fprintferr("\n");
-       (void)setrand(keep); return a;
-     }
    }
+   avma = av;
+   return (prod * N * 10 < 1);
  }
- /* Input: an ideal mod p (!= Z_K)
+ /* prime ideal decomposition of p */
-  * Output: a 2-elt representation [p, x] */
  static GEN
- prime_two_elt(GEN nf, GEN p, GEN ideal)
+ _primedec(GEN nf, GEN p)
  {
-   GEN beta,a,pf, pol = (GEN)nf[1];
+   long i,k,c,iL,N;
-   long f, N=degpol(pol), m=lg(ideal)-1;
+   GEN ex,F,L,Lpr,Ip,H,phi,mat1,T,f,g,h,p1,UN;
-   ulong av;
+   int appr;
-   if (!m) return gscalcol_i(p,N);
+   if (DEBUGLEVEL>2) (void)timer2();
+   nf = checknf(nf); T = (GEN)nf[1];
+   F = factmod(T,p);
+   ex = (GEN)F[2];
+   F  = (GEN)F[1]; F = centerlift(F);
+   if (DEBUGLEVEL>5) msgtimer("factmod");
-   /* we want v_p(Norm(beta)) = p^f, f = N-m */
+   k = lg(F);
-   av = avma; f = N-m; pf = gpuigs(p,f);
+   if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
-   ideal = centerlift(ideal);
+   {
-   ideal = concatsp(gscalcol(p,N), ideal);
+     L = cgetg(k,t_VEC);
-   ideal = ideal_better_basis(nf, ideal, p);
+     for (i=1; i<k; i++)
-   beta = gmul((GEN)nf[7], ideal);
+       L[i] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
+     if (DEBUGLEVEL>5) msgtimer("simple primedec");
+     return L;
+   }
- #if 0
+   g = (GEN)F[1];
-   a = prime_two_elt_loop(beta,pol,p,pf);
+   for (i=2; i<k; i++) g = FpX_mul(g,(GEN)F[i], p);
-   if (!a) err(bugparier, "prime_two_elt (failed)");
+   h = FpX_div(T,g,p);
- #else
+   f = FpX_red(gdivexact(gsub(gmul(g,h), T), p), p);
-   a = random_prime_two_elt_loop(beta,pol,p,pf);
- #endif
-   a = centermod(algtobasis_intern(nf,a), p);
+   N = degpol(T);
-   if (resii(divii(subres(gmul((GEN)nf[7],a),pol),pf),p) == gzero)
+   L = cgetg(N+1,t_VEC); iL = 1;
-     a[1] = laddii((GEN)a[1],p);
+   for (i=1; i<k; i++)
-   return gerepilecopy(av,a);
+     if (is_pm1(ex[i]) || signe(FpX_res(f,(GEN)F[i],p)))
- }
+       L[iL++] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
+     else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
+       ex[i] = 0;
+   if (DEBUGLEVEL>2) msgtimer("%ld Kummer factors", iL-1);
- static GEN
+   /* phi matrix of x -> x^p - x in algebra Z_K/p */
- apply_kummer(GEN nf,GEN pol,GEN e,GEN p,long N,GEN *beta)
+   Ip = pradical(nf,p,&phi);
- {
+   if (DEBUGLEVEL>2) msgtimer("pradical");
-   GEN T,p1, res = cgetg(6,t_VEC);
-   long f = degpol(pol);
-   res[1]=(long)p;
+   /* split etale algebra Z_K / (p,Ip) */
-   res[3]=(long)e;
+   h = cgetg(N+1,t_VEC);
-   res[4]=lstoi(f);
+   if (iL > 1)
-   if (f == N) /* inert */
+   { /* split off Kummer factors */
-   {
+     GEN mulbeta, beta = NULL;
-     res[2]=(long)gscalcol_i(p,N);
+     for (i=1; i<k; i++)
-     res[5]=(long)gscalcol_i(gun,N);
+       if (!ex[i]) beta = beta? FpX_mul(beta, (GEN)F[i], p): (GEN)F[i];
+     beta = FpV_red(algtobasis_i(nf,beta), p);
+     mulbeta = FpM_red(eltmul_get_table(nf, beta), p);
+     p1 = concatsp(mulbeta, Ip);
+     /* Fp-base of ideal (Ip, beta) in ZK/p */
+     h[1] = (long)FpM_image(p1, p);
    }
    else
-   {
+     h[1] = (long)Ip;
-     T = (GEN) nf[1];
-     if (ggval(subres(pol,T),p) > f)
-       pol[2] = laddii((GEN)pol[2],p);
-     res[2] = (long) algtobasis_intern(nf,pol);
-     p1 = FpX_div(T,pol,p);
+   UN = gscalcol(gun, N);
-     res[5] = (long) centermod(algtobasis_intern(nf,p1), p);
+   for (c=1; c; c--)
+   { /* Let A:= (Z_K/p) / Ip; try to split A2 := A / Im H ~ Im M2
+        H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
+     GEN M, Mi, M2, Mi2, phi2;
+     long dim;
-     if (beta)
+     H = (GEN)h[c]; k = lg(H)-1;
-       *beta = *beta? FpX_div(*beta,pol,p): p1;
+     M   = FpM_suppl(concatsp(H,UN), p);
+     Mi  = FpM_inv(M, p);
+     M2  = vecextract_i(M, k+1,N); /* M = (H|M2) invertible */
+     Mi2 = rowextract_i(Mi,k+1,N);
+     /* FIXME: FpM_mul(,M2) could be done with vecextract_p */
+     phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
+     mat1 = FpM_ker(phi2, p);
+     dim = lg(mat1)-1; /* A2 product of 'dim' fields */
+     if (dim > 1)
+     { /* phi2 v = 0 <==> a = M2 v in Ker phi */
+       GEN I,R,r,a,mula,mul2, v = (GEN)mat1[2];
+       long n;
+       a = FpM_FpV_mul(M2,v, p);
+       mula = FpM_red(eltmul_get_table(nf, a), p);
+       mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
+       R = rootmod(pol_min(mul2,p), p); /* totally split mod p */
+       n = lg(R)-1;
+       for (i=1; i<=n; i++)
+       {
+         r = lift_intern((GEN)R[i]);
+         I = gaddmat_i(negi(r), mula);
+         h[c++] = (long)FpM_image(concatsp(H, I), p);
+       }
+       if (n == dim)
+         for (i=1; i<=n; i++) { H = (GEN)h[--c]; L[iL++] = (long)H; }
+     }
+     else /* A2 field ==> H maximal, f = N-k = dim(A2) */
+       L[iL++] = (long)H;
    }
-   return res;
+   setlg(L, iL);
+   Lpr = cgetg(iL, t_VEC);
+   if (DEBUGLEVEL>2) msgtimer("splitting %ld factors",iL-1);
+   appr = use_appr(L,p,N);
+   if (DEBUGLEVEL>2 && appr) fprintferr("using the approximation theorem\n");
+   for (i=1; i<iL; i++) Lpr[i] = (long)get_pr(nf, L, i, p, appr);
+   if (DEBUGLEVEL>2) msgtimer("finding uniformizers");
+   return Lpr;
  }
- /* prime ideal decomposition of p sorted by increasing residual degree */
  GEN
  primedec(GEN nf, GEN p)
  {
-   long av=avma,tetpil,i,j,k,kbar,np,c,indice,N,lp;
+   gpmem_t av = avma;
-   GEN ex,f,list,ip,elth,h;
+   return gerepileupto(av, gen_sort(_primedec(nf,p), 0, cmp_prime_over_p));
-   GEN modfrob,algebre,algebre1,b,mat1,T;
+ }
-   GEN alpha,beta,p1,p2,unmodp,zmodp,idmodp;
-   if (DEBUGLEVEL>=3) timer2();
+ /* return [Fp[x]: Fp] */
-   nf=checknf(nf); T=(GEN)nf[1]; N=degpol(T);
+ static long
-   f=factmod(T,p); ex=(GEN)f[2];
+ ffdegree(GEN x, GEN frob, GEN p)
-   f=centerlift((GEN)f[1]); np=lg(f);
+ {
-   if (DEBUGLEVEL>=6) msgtimer("factmod");
+   gpmem_t av = avma;
+   long d, f = lg(frob)-1;
+   GEN y = x;
-   if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
+   for (d=1; d < f; d++)
    {
-     list=cgetg(np,t_VEC);
+     y = FpM_FpV_mul(frob, y, p);
-     for (i=1; i<np; i++)
+     if (gegal(y, x)) break;
-       list[i]=(long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N, NULL);
-     if (DEBUGLEVEL>=6) msgtimer("simple primedec");
-     p1=stoi(4); tetpil=avma;
-     return gerepile(av,tetpil,vecsort(list,p1));
    }
+   avma = av; return d;
+ }
-   p1 = (GEN)f[1];
+ static GEN
-   for (i=2; i<np; i++)
+ lift_to_zk(GEN v, GEN c, long N)
-     p1 = FpX_red(gmul(p1, (GEN)f[i]), p);
+ {
-   p1 = FpX_red(gdiv(gadd(gmul(p1, FpX_div(T,p1,p)), gneg(T)), p), p);
+   GEN w = zerocol(N);
-   list = cgetg(N+1,t_VEC);
+   long i, l = lg(c);
-   indice=1; beta=NULL;
+   for (i=1; i<l; i++) w[c[i]] = v[i];
-   for (i=1; i<np; i++) /* e = 1 or f[i] does not divide p1 (mod p) */
+   return w;
-     if (is_pm1(ex[i]) || signe(FpX_res(p1,(GEN)f[i],p)))
+ }
-       list[indice++] = (long)apply_kummer(nf,(GEN)f[i],(GEN)ex[i],p,N,&beta);
-   if (DEBUGLEVEL>=3) msgtimer("unramified factors");
-   /* modfrob = modified Frobenius: x -> x^p - x mod p */
+ /* return integral x = 0 mod p/pr^e, (x,pr) = 1.
-   ip = pradical(nf,p,&modfrob);
+  * Don't reduce mod p here: caller may need result mod pr^k */
-   if (DEBUGLEVEL>=3) msgtimer("pradical");
+ 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), gpowgs(p,itos(e)-1));
+ }
-   if (beta)
+ /* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
+ static GEN
+ anti_uniformizer2(GEN nf, GEN pr)
+ {
+   GEN p = (GEN)pr[1], z;
+   z = gmod(special_anti_uniformizer(nf, pr), p);
+   z = eltmul_get_table(nf, z);
+   z = hnfmodid(z, p);
+   z = idealaddtoone_i(nf,pr,z);
+   return unnf_minus_x(z);
+ }
+ #define mpr_TAU 1
+ #define mpr_FFP 2
+ #define mpr_PR  3
+ #define mpr_T   4
+ #define mpr_NFP 5
+ #define SMALLMODPR 4
+ #define LARGEMODPR 6
+ static GEN
+ modpr_TAU(GEN modpr)
+ {
+   GEN tau = (GEN)modpr[mpr_TAU];
+   if (typ(tau) == t_INT && signe(tau) == 0) tau = NULL;
+   return tau;
+ }
+ /* prh = HNF matrix, which is identity but for the first line. Return a
+  * projector to ZK / prh ~ Z/prh[1,1] */
+ GEN
+ dim1proj(GEN prh)
+ {
+   long i, N = lg(prh)-1;
+   GEN ffproj = cgetg(N+1, t_VEC);
+   GEN x, q = gcoeff(prh,1,1);
+   ffproj[1] = un;
+   for (i=2; i<=N; i++)
    {
-     beta = algtobasis_intern(nf,beta);
+     x = gcoeff(prh,1,i);
-     lp=lg(ip)-1; p1=cgetg(2*lp+N+1,t_MAT);
+     if (signe(x)) x = subii(q,x);
-     for (i=1; i<=N; i++) p1[i]=(long)element_mulid(nf,beta,i);
+     ffproj[i] = (long)x;
-     for (   ; i<=N+lp; i++)
-     {
-       p2 = (GEN) ip[i-N];
-       p1[i+lp] = (long) p2;
-       p1[i] = ldiv(element_mul(nf,lift(p2),beta),p);
-     }
-     ip = FpM_image(p1, p);
-     if (lg(ip)>N) err(bugparier,"primedec (bad pradical)");
    }
-   unmodp=gmodulsg(1,p); zmodp=gmodulsg(0,p);
+   return ffproj;
-   idmodp = idmat_intern(N,unmodp,zmodp);
+ }
-   ip = gmul(ip, unmodp);
-   h=cgetg(N+1,t_VEC); h[1]=(long)ip;
+ static GEN
-   for (c=1; c; c--)
+ modprinit(GEN nf, GEN pr, int zk)
+ {
+   gpmem_t av = avma;
+   GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c, gf;
+   long N, i, k, f;
+   nf = checknf(nf); checkprimeid(pr);
+   gf = (GEN)pr[4];
+   f = itos(gf);
+   N = degpol(nf[1]);
+   prh = prime_to_ideal(nf, pr);
+   tau = zk? gzero: anti_uniformizer2(nf, pr);
+   p = (GEN)pr[1];
+   if (f == 1)
    {
-     elth=(GEN)h[c]; k=lg(elth)-1; kbar=N-k;
+     res = cgetg(SMALLMODPR, t_COL);
-     p1 = concatsp(elth,(GEN)idmodp[1]);
+     res[mpr_TAU] = (long)tau;
-     algebre = suppl_intern(p1,idmodp);
+     res[mpr_FFP] = (long)dim1proj(prh);
-     algebre1 = cgetg(kbar+1,t_MAT);
+     res[mpr_PR ] = (long)pr; return gerepilecopy(av, res);
-     for (i=1; i<=kbar; i++) algebre1[i]=algebre[i+k];
+   }
-     b = gmul(modfrob,algebre1);
-     for (i=1;i<=kbar;i++)
+   c = cgetg(f+1, t_VECSMALL);
-       b[i] = (long) project(algebre,(GEN) b[i],k,kbar);
+   ffproj = cgetg(N+1, t_MAT);
-     mat1 = FpM_ker(lift_intern(b), p);
+   for (k=i=1; i<=N; i++)
-     if (lg(mat1)>2)
+   {
+     x = gcoeff(prh, i,i);
+     if (!is_pm1(x)) { c[k] = i; ffproj[i] = (long)_ei(N, i); k++; }
+     else
+       ffproj[i] = lneg((GEN)prh[i]);
+   }
+   ffproj = rowextract_p(ffproj, c);
+   if (! divise((GEN)nf[4], p))
+   {
+     GEN basis, proj_modT;
+     if (N == f) T = (GEN)nf[1]; /* pr inert */
+     else
+       T = primpart(gmul((GEN)nf[7], (GEN)pr[2]));
+     basis = (GEN)nf[7];
+     proj_modT = cgetg(f+1, t_MAT);
+     for (i=1; i<=f; i++)
      {
-       GEN mat2 = cgetg(k+N+1,t_MAT);
+       GEN cx, w = Q_primitive_part((GEN)basis[c[i]], &cx);
-       for (i=1; i<=k; i++) mat2[i]=elth[i];
+       w = FpX_rem(w, T, p);
-       alpha=gmul(algebre1,(GEN)mat1[2]);
+       if (cx)
-       p1 = pol_min(alpha,nf,algebre,kbar,p);
-       p1 = (GEN)factmod(p1,p)[1];
-       for (i=1; i<lg(p1); i++)
        {
-         beta = eval_pol(nf,(GEN)p1[i],alpha,algebre,algebre1);
+         cx = gmod(cx, p);
-         beta = lift_intern(beta);
+         w = FpX_Fp_mul(w, cx, p);
-         for (j=1; j<=N; j++)
-           mat2[k+j] = (long)FpV(element_mulid(nf,beta,j), p);
-         h[c] = (long) image(mat2); c++;
        }
+       proj_modT[i] = (long)pol_to_vec(w, f); /* w_i mod (T,p) */
      }
-     else
+     ffproj = FpM_mul(proj_modT,ffproj,p);
+     res = cgetg(SMALLMODPR+1, t_COL);
+     res[mpr_TAU] = (long)tau;
+     res[mpr_FFP] = (long)ffproj;
+     res[mpr_PR ] = (long)pr;
+     res[mpr_T  ] = (long)T; return gerepilecopy(av, res);
+   }
+   if (isprime(gf))
+   {
+     mul = eltmulid_get_table(nf, c[2]);
+     mul = vecextract_p(mul, c);
+   }
+   else
+   {
+     GEN v, u, u2, frob;
+     long deg,deg1,deg2;
+     /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
+     frob = cgetg(f+1, t_MAT);
+     for (i=1; i<=f; i++)
      {
-       long av1; p1 = cgetg(6,t_VEC);
+       x = element_powid_mod_p(nf,c[i],p, p);
-       list[indice++] = (long)p1;
+       frob[i] = (long)FpM_FpV_mul(ffproj, x, p);
-       p1[1]=(long)p; p1[4]=lstoi(kbar);
-       p1[2]=(long)prime_two_elt(nf,p,elth);
-       p1[5]=(long)lens(nf,p,(GEN)p1[2]);
-       av1=avma;
-       i = int_elt_val(nf,(GEN)p1[5],p,(GEN)p1[5],NULL,N-1);
-       avma=av1;
-       p1[3]=lstoi(i+1);
      }
-     if (DEBUGLEVEL>=3) msgtimer("h[%ld]",c);
+     u = _ei(f,2); k = 2;
+     deg1 = ffdegree(u, frob, p);
+     while (deg1 < f)
+     {
+       k++; u2 = _ei(f, k);
+       deg2 = ffdegree(u2, frob, p);
+       deg = clcm(deg1,deg2);
+       if (deg == deg1) continue;
+       if (deg == deg2) { deg1 = deg2; u = u2; continue; }
+       u = gadd(u, u2);
+       while (ffdegree(u, frob, p) < deg) u = gadd(u, u2);
+       deg1 = deg;
+     }
+     v = lift_to_zk(u,c,N);
+     mul = cgetg(f+1,t_MAT);
+     mul[1] = (long)v; /* assume w_1 = 1 */
+     for (i=2; i<=f; i++) mul[i] = (long)element_mulid(nf,v,c[i]);
    }
-   setlg(list, indice); tetpil=avma;
-   return gerepile(av,tetpil,gen_sort(list,0,cmp_prime_over_p));
+   /* Z_K/pr = Fp(v), mul = mul by v */
+   mul = FpM_red(mul, p);
+   mul = FpM_mul(ffproj, mul, p);
+   pow = get_powers(mul, p);
+   T = gtopolyrev(FpM_deplin(pow, p), varn(nf[1]));
+   nfproj = cgetg(f+1, t_MAT);
+   for (i=1; i<=f; i++) nfproj[i] = (long)lift_to_zk((GEN)pow[i], c, N);
+   nfproj = gmul((GEN)nf[7], nfproj);
+   setlg(pow, f+1);
+   ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
+   res = cgetg(LARGEMODPR, t_COL);
+   res[mpr_TAU] = (long)tau;
+   res[mpr_FFP] = (long)ffproj;
+   res[mpr_PR ] = (long)pr;
+   res[mpr_T  ] = (long)T;
+   res[mpr_NFP] = (long)nfproj; return gerepilecopy(av, res);
  }
+ GEN
+ nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0); }
+ GEN
+ zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1); }
+ void
+ checkmodpr(GEN modpr)
+ {
+   if (typ(modpr) != t_COL || lg(modpr) < SMALLMODPR)
+     err(talker,"incorrect modpr format");
+   checkprimeid((GEN)modpr[mpr_PR]);
+ }
+ static GEN
+ to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
+ {
+   GEN modpr = (typ(*pr) == t_COL)? *pr: modprinit(nf, *pr, zk);
+   *T = lg(modpr)==SMALLMODPR? NULL: (GEN)modpr[mpr_T];
+   *pr = (GEN)modpr[mpr_PR];
+   *p = (GEN)(*pr)[1]; return modpr;
+ }
+ GEN
+ nf_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
+   GEN modpr = to_ff_init(nf,pr,T,p,0);
+   GEN tau = modpr_TAU(modpr);
+   if (!tau) modpr[mpr_TAU] = (long)anti_uniformizer2(nf, *pr);
+   return modpr;
+ }
+ GEN
+ zk_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
+   return to_ff_init(nf,pr,T,p,1);
+ }
+ /* assume x in 'basis' form (t_COL) */
+ GEN
+ zk_to_ff(GEN x, GEN modpr)
+ {
+   GEN pr = (GEN)modpr[mpr_PR];
+   GEN p = (GEN)pr[1];
+   GEN y = gmul((GEN)modpr[mpr_FFP], x);
+   if (lg(modpr) == SMALLMODPR) return modii(y,p);
+   y = FpV_red(y, p);
+   return col_to_ff(y, varn(modpr[mpr_T]));
+ }
  /* REDUCTION Modulo a prime ideal */
+ /* assume x in t_COL form, v_pr(x) >= 0 */
+ static GEN
+ kill_denom(GEN x, GEN nf, GEN p, GEN modpr)
+ {
+   GEN cx, den = denom(x);
+   long v;
+   if (gcmp1(den)) return x;
+   v = ggval(den,p);
+   if (v)
+   {
+     GEN tau = modpr_TAU(modpr);
+     if (!tau) err(talker,"modpr initialized for integers only!");
+     x = element_mul(nf,x, element_pow(nf,tau,stoi(v)));
+   }
+   x = Q_primitive_part(x, &cx);
+   x = FpV_red(x, p);
+   if (cx) x = FpV_red(gmul(gmod(cx, p), x), p);
+   return x;
+ }
  /* x integral, reduce mod prh in HNF */
  GEN
  nfreducemodpr_i(GEN x, GEN prh)
-Line 2323  nfreducemodpr_i(GEN x, GEN prh)
+Line 2615  nfreducemodpr_i(GEN x, GEN prh)
 Line 2323  nfreducemodpr_i(GEN x, GEN prh)
 Line 2615  nfreducemodpr_i(GEN x, GEN prh)
    x[1] = lresii((GEN)x[1], p); return x;
  }
- /* for internal use */
  GEN
- nfreducemodpr(GEN nf, GEN x, GEN prhall)
+ nfreducemodpr(GEN nf, GEN x, GEN modpr)
  {
-   long i,v;
+   gpmem_t av = avma;
-   GEN p,prh,den;
+   long i;
+   GEN pr, p;
+   checkmodpr(modpr);
+   pr = (GEN)modpr[mpr_PR];
+   p = (GEN)pr[1];
+   if (typ(x) != t_COL) x = algtobasis(nf,x);
    for (i=lg(x)-1; i>0; i--)
-     if (typ(x[i]) == t_INTMOD) { x=lift_intern(x); break; }
+     if (typ(x[i]) == t_INTMOD) { x = lift(x); break; }
-   prh=(GEN)prhall[1]; p=gcoeff(prh,1,1);
+   x = kill_denom(x, nf, p, modpr);
-   den=denom(x);
+   x = ff_to_nf(zk_to_ff(x,modpr), modpr);
-   if (!gcmp1(den))
+   if (typ(x) != t_COL) x = algtobasis(nf, x);
+   return gerepileupto(av, FpV(x, p));
+ }
+ GEN
+ nf_to_ff(GEN nf, GEN x, GEN modpr)
+ {
+   gpmem_t av = avma;
+   GEN pr = (GEN)modpr[mpr_PR];
+   GEN p = (GEN)pr[1];
+   long t = typ(x);
+   if (t == t_POLMOD) { x = (GEN)x[2]; t = typ(x); }
+   switch(t)
    {
-     v=ggval(den,p);
+     case t_INT: return modii(x, p);
-     if (v) x=element_mul(nf,x,element_pow(nf,(GEN)prhall[2],stoi(v)));
+     case t_FRAC: return gmod(x, p);
-     x = gmod(x,p);
+     case t_POL: x = algtobasis(nf, x); break;
+     case t_COL: break;
+     default: err(typeer,"nf_to_ff");
    }
-   return FpV(nfreducemodpr_i(x, prh), p);
+   x = kill_denom(x, nf, p, modpr);
+   return gerepilecopy(av, zk_to_ff(x, modpr));
  }
- /* public function */
  GEN
- nfreducemodpr2(GEN nf, GEN x, GEN prhall)
+ ff_to_nf(GEN x, GEN modpr)
  {
-   long av = avma; checkprhall(prhall);
+   if (lg(modpr) < LARGEMODPR) return x;
-   if (typ(x) != t_COL) x = algtobasis(nf,x);
+   return mulmat_pol((GEN)modpr[mpr_NFP], x);
-   return gerepileupto(av, nfreducemodpr(nf,x,prhall));
  }
+ GEN
+ modprM_lift(GEN x, GEN modpr)
+ {
+   long i,j,h,l = lg(x);
+   GEN y = cgetg(l, t_MAT);
+   if (l == 1) return y;
+   h = lg(x[1]);
+   for (j=1; j<l; j++)
+   {
+     GEN p1 = cgetg(h,t_COL); y[j] = (long)p1;
+     for (i=1; i<h; i++) p1[i]=(long)ff_to_nf(gcoeff(x,i,j), modpr);
+   }
+   return y;
+ }
+ GEN
+ modprX_lift(GEN x, GEN modpr)
+ {
+   long i, l;
+   GEN z;
+   if (typ(x)!=t_POL) return gcopy(x); /* scalar */
+   l = lgef(x);
+   z = cgetg(l, t_POL); z[1] = x[1];
+   for (i=2; i<l; i++) z[i] = (long)ff_to_nf((GEN)x[i], modpr);
+   return z;
+ }
+ /* reduce the coefficients of pol modulo modpr */
+ GEN
+ modprX(GEN x, GEN nf,GEN modpr)
+ {
+   long i, l;
+   GEN z;
+   if (typ(x)!=t_POL) return nf_to_ff(nf,x,modpr);
+   l = lgef(x);
+   z = cgetg(l,t_POL); z[1] = x[1];
+   for (i=2; i<l; i++) z[i] = (long)nf_to_ff(nf,(GEN)x[i],modpr);
+   return normalizepol(z);
+ }
+ GEN
+ modprV(GEN z, GEN nf,GEN modpr)
+ {
+   long i,l = lg(z);
+   GEN x;
+   x = cgetg(l,typ(z));
+   for (i=1; i<l; i++) x[i] = (long)nf_to_ff(nf,(GEN)z[i], modpr);
+   return x;
+ }
+ /* assume z a t_VEC/t_COL/t_MAT */
+ GEN
+ modprM(GEN z, GEN nf,GEN modpr)
+ {
+   long i,l = lg(z), m;
+   GEN x;
+   if (typ(z) != t_MAT) return modprV(z,nf,modpr);
+   x = cgetg(l,t_MAT); if (l==1) return x;
+   m = lg((GEN)z[1]);
+   for (i=1; i<l; i++) x[i] = (long)modprV((GEN)z[i],nf,modpr);
+   return x;
+ }
  /*                        relative ROUND 2
   *
   * input: nf = base field K
-Line 2363  nfreducemodpr2(GEN nf, GEN x, GEN prhall)
+Line 2736  nfreducemodpr2(GEN nf, GEN x, GEN prhall)
 Line 2363  nfreducemodpr2(GEN nf, GEN x, GEN prhall)
 Line 2736  nfreducemodpr2(GEN nf, GEN x, GEN prhall)
   */
  /* given MODULES x and y by their pseudo-bases in HNF, gives a
-  * pseudo-basis of the module generated by x and y. For internal use.
+  * pseudo-basis of the module generated by x and y. */
-  */
  static GEN
  rnfjoinmodules(GEN nf, GEN x, GEN y)
  {
-Line 2372  rnfjoinmodules(GEN nf, GEN x, GEN y)
+Line 2744  rnfjoinmodules(GEN nf, GEN x, GEN y)
 Line 2372  rnfjoinmodules(GEN nf, GEN x, GEN y)
 Line 2744  rnfjoinmodules(GEN nf, GEN x, GEN y)
    GEN p1,p2,z,Hx,Hy,Ix,Iy;
    if (x == NULL) return y;
-   Hx=(GEN)x[1]; lx=lg(Hx); Ix=(GEN)x[2];
+   Hx = (GEN)x[1]; lx = lg(Hx); Ix = (GEN)x[2];
-   Hy=(GEN)y[1]; ly=lg(Hy); Iy=(GEN)y[2];
+   Hy = (GEN)y[1]; ly = lg(Hy); Iy = (GEN)y[2];
    i = lx+ly-1;
    z = (GEN)gpmalloc(sizeof(long*)*(3+2*i));
    *z = evaltyp(t_VEC)|evallg(3);
    p1 =  z+3; z[1]=(long)p1; *p1 = evaltyp(t_MAT)|evallg(i);
    p2 = p1+i; z[2]=(long)p2; *p2 = evaltyp(t_VEC)|evallg(i);
-   for (i=1; i<lx; i++) { p1[i]=Hx[i]; p2[i]=Ix[i]; }
+   for (i=1; i<lx; i++) { p1[i] = Hx[i]; p2[i] = Ix[i]; }
-   for (   ; i<lx+ly-1; i++) { p1[i]=Hy[i-lx+1]; p2[i]=Iy[i-lx+1]; }
+   p1 += lx-1;
+   p2 += lx-1;
+   for (i=1; i<ly; i++) { p1[i] = Hy[i]; p2[i] = Iy[i]; }
    x = nfhermite(nf,z); free(z); return x;
  }
- /* a usage interne, pas de gestion de pile : x et y sont des vecteurs dont
+ typedef struct {
-  * les coefficients sont les composantes sur nf[7]; avec reduction mod pr sauf
+   GEN nf, multab, modpr,T,p;
-  * si prhall=NULL
+   long h;
-  */
+ } rnfeltmod_muldata;
  static GEN
- rnfelement_mulidmod(GEN nf, GEN multab, GEN unnf, GEN x, long h, GEN prhall)
+ _mul(void *data, GEN x, GEN y/* base; ignored */)
  {
-   long j,k,N;
+   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
-   GEN p1,c,v,s,znf;
+   GEN z = x? element_mulid(D->multab,x,D->h)
+            : element_mulidid(D->multab,D->h,D->h);
-   if (h==1) return gcopy(x);
+   (void)y;
-   N = lg(x)-1; multab += (h-1)*N;
+   return FpVQX_red(z,D->T,D->p);
-   x = lift(x); v = cgetg(N+1,t_COL);
-   znf = gscalcol_i(gzero,lg(unnf)-1);
-   for (k=1; k<=N; k++)
-   {
-     s = gzero;
-     for (j=1; j<=N; j++)
-       if (!gcmp0(p1 = (GEN)x[j]) && !gcmp0(c = gcoeff(multab,k,j)))
-       {
-         if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-         s = gadd(s,p1);
-       }
-     if (s == gzero) s = znf;
-     else
-       if (prhall) s = nfreducemodpr(nf,s,prhall);
-     v[k] = (long)s;
-   }
-   return v;
  }
- /* a usage interne, pas de gestion de pile : x est un vecteur dont
-  * les coefficients sont les composantes sur nf[7]
-  */
  static GEN
- rnfelement_sqrmod(GEN nf, GEN multab, GEN unnf, GEN x, GEN prhall)
+ _sqr(void *data, GEN x)
  {
-   long i,j,k,n;
+   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
-   GEN p1,c,z,s;
+   GEN z = x? sqr_by_tab(D->multab,x)
+            : element_mulidid(D->multab,D->h,D->h);
-   n=lg(x)-1; x=lift(x); z=cgetg(n+1,t_COL);
+   return FpVQX_red(z,D->T,D->p);
-   for (k=1; k<=n; k++)
-   {
-     if (k == 1)
-       s = element_sqr(nf,(GEN)x[1]);
-     else
-       s = gmul2n(element_mul(nf,(GEN)x[1],(GEN)x[k]), 1);
-     for (i=2; i<=n; i++)
-     {
-       c = gcoeff(multab,k,(i-1)*n+i);
-       if (!gcmp0(c))
-       {
-         p1=element_sqr(nf,(GEN)x[i]);
-         if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-         s = gadd(s,p1);
-       }
-       for (j=i+1; j<=n; j++)
-       {
-         c = gcoeff(multab,k,(i-1)*n+j);
-         if (!gcmp0(c))
-         {
-           p1=gmul2n(element_mul(nf,(GEN)x[i],(GEN)x[j]),1);
-           if (!gegal(c,unnf)) p1 = element_mul(nf,p1,c);
-           s = gadd(s,p1);
-         }
-       }
-     }
-     if (prhall) s = nfreducemodpr(nf,s,prhall);
-     z[k]=(long)s;
-   }
-   return z;
  }
- /* Compute x^n mod pr in the extension, assume n >= 0 [cf puissii()]*/
+ /* Compute x^n mod pr in the extension, assume n >= 0 */
  static GEN
- rnfelementid_powmod(GEN nf, GEN multab, GEN matId, long h, GEN n, GEN prhall)
+ rnfelementid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long i,k,m;
+   GEN y;
-   GEN y,p1, unrnf=(GEN)matId[1], unnf=(GEN)unrnf[1];
+   rnfeltmod_muldata D;
-   if (!signe(n)) return unrnf;
+   if (!signe(n)) return gun;
-   y = (GEN)matId[h]; p1 = n+2; m = *p1;
-   k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k;
+   D.multab = multab;
-   for (i=lgefint(n)-2;;)
+   D.h = h;
-   {
+   D.T = T;
-     for (; k; m<<=1,k--)
+   D.p = p;
-     {
+   y = leftright_pow(NULL, n, (void*)&D, &_sqr, &_mul);
-       y = rnfelement_sqrmod(nf,multab,unnf,y,prhall);
-       if (m < 0) y = rnfelement_mulidmod(nf,multab,unnf,y,h,prhall);
-     }
-     if (--i == 0) break;
-     m = *++p1; k = BITS_IN_LONG;
-   }
    return gerepilecopy(av, y);
  }
+ #if 0
  GEN
  mymod(GEN x, GEN p)
  {
-Line 2501  mymod(GEN x, GEN p)
+Line 2820  mymod(GEN x, GEN p)
 Line 2501  mymod(GEN x, GEN p)
 Line 2820  mymod(GEN x, GEN p)
    for (i=1; i<lx; i++) y[i] = (long)mymod((GEN)x[i],p);
    return y;
  }
+ #endif
+ #define FqX_div(x,y,T,p) FpXQX_divres((x),(y),(T),(p),NULL)
+ #define FqX_mul FpXQX_mul
+ /* relative Dedekind criterion over nf, applied to the order defined by a
+  * root of irreducible polynomial P, modulo the prime ideal pr. Returns
+  * [flag,basis,val], where basis is a pseudo-basis of the enlarged order,
+  * flag is 1 iff this order is pr-maximal, and val is the valuation at pr of
+  * the order discriminant */
  static GEN
- rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, GEN matId)
+ rnfdedekind_i(GEN nf,GEN P,GEN pr, GEN disc)
  {
-   long av=avma,tetpil,av1,lim,i,j,k,n,v1,v2,vpol,m,cmpt,sep;
+   long vt, r, d, n, m, i, j;
-   GEN p,q,q1,prhall,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
+   gpmem_t av = avma;
+   GEN Prd,p1,p2,p,tau,g,matid;
+   GEN modpr,res,h,k,base,nfT,T,gzk,hzk;
+   if (!disc) disc = discsr(P);
+   P = lift(P);
+   nf = checknf(nf); nfT = (GEN)nf[1];
+   res = cgetg(4,t_VEC);
+   modpr = nf_to_ff_init(nf,&pr, &T, &p);
+   tau = gmul((GEN)nf[7], (GEN)pr[5]);
+   n = degpol(nfT);
+   m = degpol(P);
+   Prd = modprX(P, nf, modpr);
+   p1 = (GEN)FqX_factor(Prd,T,p)[1];
+   r = lg(p1); if (r < 2) err(constpoler,"rnfdedekind");
+   g = (GEN)p1[1];
+   for (i=2; i<r; i++) g = FqX_mul(g, (GEN)p1[i], T, p);
+   h = FqX_div(Prd,g, T, p);
+   gzk = modprX_lift(g, modpr);
+   hzk = modprX_lift(h, modpr);
+   k = gsub(P, RXQX_mul(gzk,hzk, nfT));
+   k = gdiv(RXQX_mul(tau,k,nfT), p);
+   k = modprX(k, nf, modpr);
+   p2 = FqX_gcd(g,h,  T,p);
+   k  = FqX_gcd(p2,k, T,p); d = degpol(k);  /* <= m */
+   vt = idealval(nf,disc,pr) - 2*d;
+   res[3] = lstoi(vt);
+   res[1] = (!d || vt<=1)? un: zero;
+   base = cgetg(3,t_VEC);
+   p1 = cgetg(m+d+1,t_MAT); base[1] = (long)p1;
+   p2 = cgetg(m+d+1,t_VEC); base[2] = (long)p2;
+  /* if d > 0, base[2] temporarily multiplied by p, for the final nfhermitemod
+   * [ which requires integral ideals ] */
+   matid = gscalmat(d? p: gun, n);
+   for (j=1; j<=m; j++)
+   {
+     p1[j] = (long)_ei(m, j);
+     p2[j] = (long)matid;
+   }
+   if (d)
+   {
+     GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */
+     GEN X = polx[varn(P)], pal;
+     pal = FqX_div(Prd,k, T,p);
+     pal = modprX_lift(pal, modpr);
+     for (   ; j<=m+d; j++)
+     {
+       p1[j] = (long)pol_to_vec(pal,m);
+       p2[j] = (long)prinvp;
+       pal = RXQX_rem(RXQX_mul(pal,X,nfT),P,nfT);
+     }
+     /* the modulus is integral */
+     base = nfhermitemod(nf,base, gmul(gpowgs(p, m-d),
+                                       idealpows(nf, prinvp, d)));
+     base[2] = ldiv((GEN)base[2], p); /* cancel the factor p */
+   }
+   res[2] = (long)base; return gerepilecopy(av, res);
+ }
+ GEN
+ rnfdedekind(GEN nf,GEN P0,GEN pr)
+ {
+   return rnfdedekind_i(nf,P0,pr,NULL);
+ }
+ static GEN
+ rnfordmax(GEN nf, GEN pol, GEN pr, GEN disc)
+ {
+   gpmem_t av=avma,av1,lim;
+   long i,j,k,n,v1,v2,vpol,m,cmpt,sep;
+   GEN p,T,q,q1,modpr,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
    GEN pip,baseIp,baseOp,alpha,matprod,alphainv,matC,matG,vecpro,matH;
-   GEN neworder,H,Hid,alphalistinv,alphalist,prhinv;
+   GEN neworder,H,Hid,alphalistinv,alphalist,prhinv,nfT,id,rnfId;
    if (DEBUGLEVEL>1) fprintferr(" treating %Z\n",pr);
-   prhall=nfmodprinit(nf,pr);
+   modpr = nf_to_ff_init(nf,&pr,&T,&p);
-   q=cgetg(3,t_VEC); q[1]=(long)pr;q[2]=(long)prhall;
+   p1 = rnfdedekind_i(nf,pol,modpr,disc);
-   p1=rnfdedekind(nf,pol,q);
    if (gcmp1((GEN)p1[1])) return gerepilecopy(av,(GEN)p1[2]);
-   sep=itos((GEN)p1[3]);
+   sep = itos((GEN)p1[3]);
-   A=gmael(p1,2,1);
+   A = gmael(p1,2,1);
-   I=gmael(p1,2,2);
+   I = gmael(p1,2,2);
-   n=degpol(pol); vpol=varn(pol);
+   pip = basistoalg(nf, (GEN)pr[2]);
-   p=(GEN)pr[1]; q=powgi(p,(GEN)pr[4]); pip=(GEN)pr[2];
+   nfT = (GEN)nf[1];
-   q1=q; while (cmpis(q1,n)<0) q1=mulii(q1,q);
+   n = degpol(pol); vpol = varn(pol);
+   q = T? gpowgs(p,degpol(T)): p;
+   q1 = q; while (cmpis(q1,n) < 0) q1 = mulii(q1,q);
+   rnfId = idmat(degpol(pol));
+   id    = idmat(degpol(nfT));
-   multab=cgetg(n*n+1,t_MAT);
+   multab = cgetg(n*n+1,t_MAT);
-   for (j=1; j<=n*n; j++) multab[j]=lgetg(n+1,t_COL);
+   for (j=1; j<=n*n; j++) multab[j] = lgetg(n+1,t_COL);
-   prhinv = idealinv(nf,(GEN)prhall[1]);
+   prhinv = idealinv(nf, pr);
-   alphalistinv=cgetg(n+1,t_VEC);
+   alphalistinv = cgetg(n+1,t_VEC);
-   alphalist=cgetg(n+1,t_VEC);
+   alphalist    = cgetg(n+1,t_VEC);
-   A1=cgetg(n+1,t_MAT);
+   A1 = cgetg(n+1,t_MAT);
-   av1=avma; lim=stack_lim(av1,1);
+   av1 = avma; lim = stack_lim(av1,1);
    for(cmpt=1; ; cmpt++)
    {
-     if (DEBUGLEVEL>1)
+     if (DEBUGLEVEL>1) fprintferr("    %ld%s pass\n",cmpt,eng_ord(cmpt));
+     for (j=1; j<=n; j++)
      {
-       fprintferr("    %ld%s pass\n",cmpt,eng_ord(cmpt));
+       if (gegal((GEN)I[j],id))
-       flusherr();
+       {
-     }
+         alphalist[j] = alphalistinv[j] = un;
-     for (i=1; i<=n; i++)
+         A1[j] = A[j]; continue;
-     {
+       }
-       if (gegal((GEN)I[i],id)) alphalist[i]=alphalistinv[i]=(long)unnf;
+       p1 = ideal_two_elt(nf,(GEN)I[j]);
+       if (gcmp0((GEN)p1[1]))
+         alpha = (GEN)p1[2];
        else
        {
-         p1=ideal_two_elt(nf,(GEN)I[i]);
+         v1 = element_val(nf,(GEN)p1[1],pr);
-         v1=gcmp0((GEN)p1[1])? EXP220
+         v2 = element_val(nf,(GEN)p1[2],pr);
-                             : element_val(nf,(GEN)p1[1],pr);
+         alpha = (v1>v2)? (GEN)p1[2]: (GEN)p1[1];
-         v2=element_val(nf,(GEN)p1[2],pr);
-         if (v1>v2) p2=(GEN)p1[2]; else p2=(GEN)p1[1];
-         alphalist[i]=(long)p2;
-         alphalistinv[i]=(long)element_inv(nf,p2);
        }
-     }
+       alphalist[j]    = (long)alpha;
-     for (j=1; j<=n; j++)
+       alphalistinv[j] = (long)element_inv(nf,alpha);
-     {
-       p1=cgetg(n+1,t_COL); A1[j]=(long)p1;
+       p1 = cgetg(n+1,t_COL); A1[j] = (long)p1;
        for (i=1; i<=n; i++)
-         p1[i] = (long)element_mul(nf,gcoeff(A,i,j),(GEN)alphalist[j]);
+         p1[i] = (long)element_mul(nf,gcoeff(A,i,j),alpha);
      }
-     Aa=basistoalg(nf,A1); Aainv=lift_intern(ginv(Aa));
+     Aa = basistoalg(nf,A1);
-     Aaa=mat_to_vecpol(Aa,vpol);
+     Aainv = lift_intern(ginv(Aa));
-     for (i=1; i<=n; i++) for (j=i; j<=n; j++)
+     Aaa = lift_intern(mat_to_vecpol(Aa,vpol));
-     {
+     for (i=1; i<=n; i++)
-       long tp;
+       for (j=i; j<=n; j++)
-       p1 = lift_intern(gres(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol));
-       tp = typ(p1);
-       if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
-         p2 = gmul(p1, (GEN)Aainv[1]);
-       else
-         p2 = gmul(Aainv, pol_to_vec(p1, n));
-       p3 = algtobasis(nf,p2);
-       for (k=1; k<=n; k++)
        {
-         coeff(multab,k,(i-1)*n+j) = p3[k];
+         long tp;
-         coeff(multab,k,(j-1)*n+i) = p3[k];
+         p1 = RXQX_rem(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol, nfT);
+         tp = typ(p1);
+         if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
+           p2 = gmul(p1, (GEN)Aainv[1]);
+         else
+           p2 = mulmat_pol(Aainv, p1);
+         for (k=1; k<=n; k++)
+         {
+           GEN c = gres((GEN)p2[k], nfT);
+           coeff(multab,k,(i-1)*n+j) = (long)c;
+           coeff(multab,k,(j-1)*n+i) = (long)c;
+         }
        }
-     }
+     multabmod = modprM(multab,nf,modpr);
-     R=cgetg(n+1,t_MAT); multabmod = mymod(multab,p);
+     R = cgetg(n+1,t_MAT);
-     R[1] = matId[1];
+     R[1] = rnfId[1];
      for (j=2; j<=n; j++)
-       R[j] = (long) rnfelementid_powmod(nf,multabmod,matId, j,q1,prhall);
+       R[j] = (long) rnfelementid_powmod(multabmod,j,q1,T,p);
-     baseIp = nfkermodpr(nf,R,prhall);
+     baseIp = FqM_ker(R,T,p);
-     baseOp = lift_intern(nfsuppl(nf,baseIp,n,prhall));
+     baseOp = lg(baseIp)==1? rnfId: FqM_suppl(baseIp,T,p);
-     alpha=cgetg(n+1,t_MAT);
+     alpha = modprM_lift(baseOp, modpr);
-     for (j=1; j<lg(baseIp); j++) alpha[j]=baseOp[j];
+     for (j=1; j<lg(baseIp); j++)
+     {
+       p1 = (GEN)alpha[j];
+       for (i=1; i<=n; i++) p1[i] = (long)to_polmod((GEN)p1[i], nfT);
+     }
      for (   ; j<=n; j++)
      {
-       p1=cgetg(n+1,t_COL); alpha[j]=(long)p1;
+       p1 = (GEN)alpha[j];
-       for (i=1; i<=n; i++)
+       for (i=1; i<=n; i++) p1[i] = lmul(pip, (GEN)p1[i]);
-         p1[i]=(long)element_mul(nf,pip,gcoeff(baseOp,i,j));
      }
-     matprod=cgetg(n+1,t_MAT);
+     alphainv = lift_intern(ginv(alpha));
+     alpha = lift_intern(alpha);
+     matprod = cgetg(n+1,t_MAT);
      for (j=1; j<=n; j++)
      {
-       p1=cgetg(n+1,t_COL); matprod[j]=(long)p1;
+       p1 = cgetg(n+1,t_COL); matprod[j] = (long)p1;
        for (i=1; i<=n; i++)
-       {
+         p1[i] = (long)gmod(element_mulid(multab, (GEN)alpha[i],j), nfT);
-         p2 = rnfelement_mulidmod(nf,multab,unnf, (GEN)alpha[i],j, NULL);
-         for (k=1; k<=n; k++)
-           p2[k] = lmul((GEN)nf[7], (GEN)p2[k]);
-         p1[i] = (long)p2;
-       }
      }
-     alphainv = lift_intern(ginv(basistoalg(nf,alpha)));
      matC = cgetg(n+1,t_MAT);
      for (j=1; j<=n; j++)
      {
-       p1=cgetg(n*n+1,t_COL); matC[j]=(long)p1;
+       p1 = cgetg(n*n+1,t_COL); matC[j] = (long)p1;
        for (i=1; i<=n; i++)
        {
          p2 = gmul(alphainv, gcoeff(matprod,i,j));
          for (k=1; k<=n; k++)
-           p1[(i-1)*n+k]=(long)nfreducemodpr(nf,algtobasis(nf,(GEN)p2[k]),prhall);
+         {
+           GEN c = gres((GEN)p2[k], nfT);
+           p1[(i-1)*n+k] = (long)nf_to_ff(nf,c,modpr);
+         }
        }
      }
-     matG=nfkermodpr(nf,matC,prhall); m=lg(matG)-1;
+     matG = FqM_ker(matC,T,p); m = lg(matG)-1;
-     vecpro=cgetg(3,t_VEC);
+     matG = modprM_lift(matG, modpr);
-     p1=cgetg(n+m+1,t_MAT); vecpro[1]=(long)p1;
+     vecpro = cgetg(3,t_VEC);
-     p2=cgetg(n+m+1,t_VEC); vecpro[2]=(long)p2;
+     vecpro[1] = (long)concatsp(matG,rnfId);
-     for (j=1; j<=m; j++)
-     {
+     p2 = cgetg(n+m+1,t_VEC);
-       p1[j] = llift((GEN)matG[j]);
+     vecpro[2] = (long)p2;
-       p2[j] = (long)prhinv;
+     for (j=1; j<=m; j++) p2[j] = (long)prhinv;
-     }
-     p1 += m;
      p2 += m;
      for (j=1; j<=n; j++)
-     {
-       p1[j] = matId[j];
        p2[j] = (long)idealmul(nf,(GEN)I[j],(GEN)alphalistinv[j]);
-     }
-     matH=nfhermite(nf,vecpro);
-     p1=algtobasis(nf,gmul(Aa,basistoalg(nf,(GEN)matH[1])));
-     p2=(GEN)matH[2];
-     tetpil=avma; neworder=cgetg(3,t_VEC);
+     matH = nfhermite(nf,vecpro);
-     H=cgetg(n+1,t_MAT); Hid=cgetg(n+1,t_VEC);
+     p1 = algtobasis(nf, gmul(Aa, basistoalg(nf,(GEN)matH[1])));
+     p2 = (GEN)matH[2];
+     neworder = cgetg(3,t_VEC);
+     H   = cgetg(n+1,t_MAT);
+     Hid = cgetg(n+1,t_VEC);
      for (j=1; j<=n; j++)
      {
+       alphainv = (GEN)alphalistinv[j];
+       if (alphainv == gun) { H[j] = p1[j]; Hid[j] = p2[j]; continue; }
        p3=cgetg(n+1,t_COL); H[j]=(long)p3;
        for (i=1; i<=n; i++)
-         p3[i]=(long)element_mul(nf,gcoeff(p1,i,j),(GEN)alphalistinv[j]);
+         p3[i] = (long)element_mul(nf,gcoeff(p1,i,j),alphainv);
-       Hid[j]=(long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
+       Hid[j] = (long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
      }
-     if (DEBUGLEVEL>3)
+     if (DEBUGLEVEL>3) { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
-       { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
      if (sep == 2 || gegal(I,Hid))
      {
-       neworder[1]=(long)H; neworder[2]=(long)Hid;
+       neworder[1] = (long)H;
-       return gerepile(av,tetpil,neworder);
+       neworder[2] = (long)Hid;
+       return gerepilecopy(av, neworder);
      }
-     A=H; I=Hid;
+     A = H; I = Hid;
      if (low_stack(lim, stack_lim(av1,1)) || (cmpt & 3) == 0)
      {
        GEN *gptr[2]; gptr[0]=&A; gptr[1]=&I;
-Line 2665  rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, G
+Line 3077  rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, G
 Line 2665  rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, G
 Line 3077  rnfordmax(GEN nf, GEN pol, GEN pr, GEN unnf, GEN id, G
  static void
  check_pol(GEN x, long v)
  {
-   long i,tx, lx = lg(x);
+   long i,tx, lx = lgef(x);
    if (varn(x) != v)
      err(talker,"incorrect variable in rnf function");
    for (i=2; i<lx; i++)
-Line 2680  GEN
+Line 3092  GEN
 Line 2680  GEN
 Line 3092  GEN
  fix_relative_pol(GEN nf, GEN x, int chk_lead)
  {
    GEN xnf = (typ(nf) == t_POL)? nf: (GEN)nf[1];
-   long i, vnf = varn(xnf), lx = lg(x);
+   long i, vnf = varn(xnf), lx = lgef(x);
    if (typ(x) != t_POL || varn(x) >= vnf)
      err(talker,"incorrect polynomial in rnf function");
    x = dummycopy(x);
-Line 2703  fix_relative_pol(GEN nf, GEN x, int chk_lead)
+Line 3115  fix_relative_pol(GEN nf, GEN x, int chk_lead)
 Line 2703  fix_relative_pol(GEN nf, GEN x, int chk_lead)
 Line 3115  fix_relative_pol(GEN nf, GEN x, int chk_lead)
  static GEN
  rnfround2all(GEN nf, GEN pol, long all)
  {
-   long av=avma,tetpil,i,j,n,N,nbidp,vpol,*ep;
+   gpmem_t av = avma;
-   GEN p1,p2,p3,p4,polnf,list,unnf,id,matId,I,W,pseudo,y,discpol,d,D,sym;
+   long i,j,n,N,nbidp,vpol,*ep;
+   GEN A,p1,p2,nfT,list,id,I,W,pseudo,y,d,D,sym,unnf,disc;
-   nf=checknf(nf); polnf=(GEN)nf[1]; vpol=varn(pol);
+   nf=checknf(nf); nfT=(GEN)nf[1]; vpol=varn(pol);
    pol = fix_relative_pol(nf,pol,1);
-   N=degpol(polnf); n=degpol(pol); discpol=discsr(pol);
+   N = degpol(nfT);
-   list=idealfactor(nf,discpol); ep=(long*)list[2]; list=(GEN)list[1];
+   n = degpol(pol);
-   nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i]=itos((GEN)ep[i]);
+   disc = discsr(pol); pol = lift(pol);
+   list = idealfactor(nf, disc);
+   ep  = (long*)list[2];
+   list= (GEN)  list[1];
+   nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i] = itos((GEN)ep[i]);
    if (DEBUGLEVEL>1)
    {
      fprintferr("Ideals to consider:\n");
-Line 2718  rnfround2all(GEN nf, GEN pol, long all)
+Line 3135  rnfround2all(GEN nf, GEN pol, long all)
 Line 2718  rnfround2all(GEN nf, GEN pol, long all)
 Line 3135  rnfround2all(GEN nf, GEN pol, long all)
        if (ep[i]>1) fprintferr("%Z^%ld\n",list[i],ep[i]);
      flusherr();
    }
-   id=idmat(N); unnf=gscalcol_i(gun,N);
+   id = idmat(N); unnf = (GEN)id[1];
-   matId=idmat_intern(n,unnf, gscalcol_i(gzero,N));
    pseudo = NULL;
    for (i=1; i<=nbidp; i++)
      if (ep[i] > 1)
      {
-       y=rnfordmax(nf,pol,(GEN)list[i],unnf,id,matId);
+       y = rnfordmax(nf, pol, (GEN)list[i], disc);
        pseudo = rnfjoinmodules(nf,pseudo,y);
      }
-   if (!pseudo)
+   if (pseudo)
    {
-     I=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i]=(long)id;
+     A = (GEN)pseudo[1];
-     pseudo=cgetg(3,t_VEC); pseudo[1]=(long)matId; pseudo[2]=(long)I;
+     I = (GEN)pseudo[2];
    }
-   W=gmodulcp(mat_to_vecpol(basistoalg(nf,(GEN)pseudo[1]),vpol),pol);
+   else
-   p2=cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j]=lgetg(n+1,t_COL);
+   {
-   sym=polsym(pol,n-1);
+     I = cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i] = (long)id;
+     A = idmat_intern(n, unnf, zerocol(N));
+   }
+   W = mat_to_vecpol(lift_intern(basistoalg(nf,A)), vpol);
+   sym = polsym_gen(pol, NULL, n-1, nfT, NULL);
+   p2 = cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j] = lgetg(n+1,t_COL);
    for (j=1; j<=n; j++)
      for (i=j; i<=n; i++)
      {
-       p1 = lift_intern(gmul((GEN)W[i],(GEN)W[j]));
+       p1 = RXQX_mul((GEN)W[i],(GEN)W[j], nfT);
-       coeff(p2,j,i)=coeff(p2,i,j)=(long)quicktrace(p1,sym);
+       p1 = RXQX_rem(p1, pol, nfT);
+       p1 = simplify_i(quicktrace(p1,sym));
+       coeff(p2,j,i)=coeff(p2,i,j) = (long)p1;
      }
-   d = algtobasis_intern(nf,det(p2));
+   d = algtobasis_i(nf, det(p2));
-   I=(GEN)pseudo[2];
+   i=1; while (i<=n && gegal((GEN)I[i], id)) i++;
-   i=1; while (i<=n && gegal((GEN)I[i],id)) i++;
+   if (i > n) D = id;
-   if (i>n) D=id;
    else
    {
      D = (GEN)I[i];
      for (i++; i<=n; i++)
-       if (!gegal((GEN)I[i],id)) D = idealmul(nf,D,(GEN)I[i]);
+       if (!gegal((GEN)I[i], id)) D = idealmul(nf,D,(GEN)I[i]);
      D = idealpow(nf,D,gdeux);
    }
-   p4=gun; p3=auxdecomp(content(d),0);
+   p2 = core2partial(content(d));
-   for (i=1; i<lg(p3[1]); i++)
+   p2 = sqri((GEN)p2[2]);
-     p4 = gmul(p4, gpuigs(gcoeff(p3,i,1), itos(gcoeff(p3,i,2))>>1));
-   p4 = gsqr(p4); tetpil=avma;
    i = all? 2: 0;
-   p1=cgetg(3 + i,t_VEC);
+   p1=cgetg(3 + i, t_VEC);
-   if (i) { p1[1]=lcopy((GEN)pseudo[1]); p1[2]=lcopy(I); }
+   if (i) { p1[1] = lcopy(A); p1[2] = lcopy(I); }
    p1[1+i] = (long)idealmul(nf,D,d);
-   p1[2+i] = ldiv(d,p4);
+   p1[2+i] = (long)gdiv(d, p2);
-   return gerepile(av,tetpil,p1);
+   return gerepileupto(av,p1);
  }
  GEN
-Line 2777  rnfdiscf(GEN nf, GEN pol)
+Line 3198  rnfdiscf(GEN nf, GEN pol)
 Line 2777  rnfdiscf(GEN nf, GEN pol)
 Line 3198  rnfdiscf(GEN nf, GEN pol)
    return rnfround2all(nf,pol,0);
  }
- /* given bnf as output by buchinit and a pseudo-basis of an order
-  * in HNF [A,I] (or [A,I,D,d] it does not matter), tries to simplify the
-  * HNF as much as possible. The resulting matrix will be upper triangular
-  * but the diagonal coefficients will not be equal to 1. The ideals
-  * are guaranteed to be integral and primitive.
-  */
  GEN
+ gen_if_principal(GEN bnf, GEN x)
+ {
+   gpmem_t av = avma;
+   GEN z = isprincipalall(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
+   if (typ(z) == t_INT) { avma = av; return NULL; }
+   return z;
+ }
+ /* given bnf and a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it
+  * does not matter), tries to simplify the HNF as much as possible. The
+  * resulting matrix will be upper triangular but the diagonal coefficients
+  * will not be equal to 1. The ideals are guaranteed to be integral and
+  * primitive. */
+ GEN
  rnfsimplifybasis(GEN bnf, GEN order)
  {
-   long av=avma,tetpil,j,N,n;
+   gpmem_t av = avma;
+   long j, n, l;
    GEN p1,id,Az,Iz,nf,A,I;
-   bnf = checkbnf(bnf);
+   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
    if (typ(order)!=t_VEC || lg(order)<3)
      err(talker,"not a pseudo-basis in nfsimplifybasis");
-   A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1; nf=(GEN)bnf[7];
+   A = (GEN)order[1];
-   N=degpol(nf[1]); id=idmat(N); Iz=cgetg(n+1,t_VEC); Az=cgetg(n+1,t_MAT);
+   I = (GEN)order[2]; n = lg(A)-1;
+   id = idmat(degpol(nf[1]));
+   Iz = cgetg(n+1,t_VEC);
+   Az = cgetg(n+1,t_MAT);
    for (j=1; j<=n; j++)
    {
-     if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; }
+     if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; continue; }
-     else
+     Iz[j] = (long)Q_primitive_part((GEN)I[j], &p1);
+     Az[j] = p1? lmul((GEN)A[j],p1): A[j];
+     if (p1 && gegal((GEN)Iz[j], id)) continue;
+     p1 = gen_if_principal(bnf, (GEN)Iz[j]);
+     if (p1)
      {
-       p1=content((GEN)I[j]);
+       Iz[j] = (long)id;
-       if (!gcmp1(p1))
+       Az[j] = (long)element_mulvec(nf,p1,(GEN)Az[j]);
-       {
-         Iz[j]=(long)gdiv((GEN)I[j],p1); Az[j]=lmul((GEN)A[j],p1);
-       }
-       else Az[j]=A[j];
-       if (!gegal((GEN)Iz[j],id))
-       {
-         p1=isprincipalgen(bnf,(GEN)Iz[j]);
-         if (gcmp0((GEN)p1[1]))
-         {
-           p1=(GEN)p1[2]; Iz[j]=(long)id;
-           Az[j]=(long)element_mulvec(nf,p1,(GEN)Az[j]);
-         }
-       }
      }
    }
-   tetpil=avma; p1=cgetg(lg(order),t_VEC); p1[1]=lcopy(Az); p1[2]=lcopy(Iz);
+   l = lg(order); p1 = cgetg(l, t_VEC);
-   for (j=3; j<lg(order); j++) p1[j]=lcopy((GEN)order[j]);
+   p1[1] = (long)Az;
-   return gerepile(av,tetpil,p1);
+   p1[2] = (long)Iz; for (j=3; j<l; j++) p1[j] = order[j];
+   return gerepilecopy(av, p1);
  }
  GEN
  rnfdet2(GEN nf, GEN A, GEN I)
  {
-   long av,tetpil,i;
+   gpmem_t av,tetpil;
+   long i;
    GEN p1;
    nf=checknf(nf); av = tetpil = avma;
-Line 2847  rnfdet0(GEN nf, GEN x, GEN y)
+Line 3275  rnfdet0(GEN nf, GEN x, GEN y)
 Line 2847  rnfdet0(GEN nf, GEN x, GEN y)
 Line 3275  rnfdet0(GEN nf, GEN x, GEN y)
    return y? rnfdet2(nf,x,y): rnfdet(nf,x);
  }
- /* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it does
+ /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
-  * not matter), gives an nxn matrix (not in HNF) of a pseudo-basis and
+    t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
-  * an ideal vector [id,id,...,id,I] such that order=nf[7]^(n-1)xI.
+ static GEN
-  * Since it uses the approximation theorem, can be long.
+ nfidealdet1(GEN nf, GEN a, GEN b)
-  */
+ {
+   gpmem_t av = avma;
+   GEN x,p1,res,u,v,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);
+   u = (GEN)p1[1];
+   v = (GEN)p1[2];
+   res = cgetg(5,t_VEC);
+   res[1] = (long)gmul(x,da);
+   res[2] = (long)gdiv(v,db);
+   res[3] = lnegi(db);
+   res[4] = (long)element_div(nf, u, (GEN)res[1]);
+   return gerepilecopy(av,res);
+ }
+ static GEN
+ get_order(GEN nf, GEN O, char *s)
+ {
+   if (typ(O) == t_POL)
+     return rnfpseudobasis(nf, O);
+   if (typ(O)!=t_VEC || lg(O) < 3)
+     err(talker,"not a pseudo-matrix in %s", s);
+   return O;
+ }
+ /* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d]), gives an
+  * n x n matrix (not in HNF) of a pseudo-basis and an ideal vector
+  * [id,id,...,id,I] such that order = nf[7]^(n-1) x I.
+  * Uses the approximation theorem ==> slow. */
  GEN
  rnfsteinitz(GEN nf, GEN order)
  {
-   long av=avma,tetpil,i,n;
+   gpmem_t av = avma;
+   long i,n,l;
    GEN Id,A,I,p1,a,b;
    nf = checknf(nf);
    Id = idmat(degpol(nf[1]));
-   if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
+   order = get_order(nf, order, "rnfsteinitz");
-   if (typ(order)!=t_VEC || lg(order)<3)
+   A = matalgtobasis(nf, (GEN)order[1]);
-     err(talker,"not a pseudo-matrix in rnfsteinitz");
+   I = dummycopy((GEN)order[2]); n=lg(A)-1;
-   A=dummycopy((GEN)order[1]);
-   I=dummycopy((GEN)order[2]); n=lg(A)-1;
    if (typ(A) != t_MAT || typ(I) != t_VEC || lg(I) != n+1)
      err(typeer,"rnfsteinitz");
    for (i=1; i<n; i++)
    {
-     a = (GEN)I[i];
+     GEN c1,c2;
-     if (!gegal(a,Id))
+     a = (GEN)I[i]; if (gegal(a,Id)) continue;
+     c1 = (GEN)A[i];
+     c2 = (GEN)A[i+1];
+     b = (GEN)I[i+1];
+     if (gegal(b,Id))
      {
-       GEN c1 = (GEN)A[i];
+       A[i]  = (long)c2;
-       GEN c2 = (GEN)A[i+1];
+       A[i+1]= lneg(c1);
-       b = (GEN)I[i+1];
+       I[i]  = (long)b;
-       if (gegal(b,Id))
+       I[i+1]= (long)a;
-       {
-         A[i]  = (long)c2;
-         A[i+1]= lneg(c1);
-         I[i]  = (long)b;
-         I[i+1]= (long)a;
-       }
-       else
-       {
-         p1 = nfidealdet1(nf,a,b);
-         A[i]  = ladd(element_mulvec(nf,(GEN)p1[1], c1),
-                      element_mulvec(nf,(GEN)p1[2], c2));
-         A[i+1]= ladd(element_mulvec(nf,(GEN)p1[3], c1),
-                      element_mulvec(nf,(GEN)p1[4], c2));
-         I[i]  =(long)Id;
-         I[i+1]=(long)idealmul(nf,a,b); p1 = content((GEN)I[i+1]);
-         if (!gcmp1(p1))
-         {
-           I[i+1] = ldiv((GEN)I[i+1],p1);
-           A[i+1] = lmul(p1,(GEN)A[i+1]);
-         }
-       }
      }
+     else
+     {
+       p1 = nfidealdet1(nf,a,b);
+       A[i]  = ladd(element_mulvec(nf, (GEN)p1[1], c1),
+                    element_mulvec(nf, (GEN)p1[2], c2));
+       A[i+1]= ladd(element_mulvec(nf, (GEN)p1[3], c1),
+                    element_mulvec(nf, (GEN)p1[4], c2));
+       I[i]  = (long)Id;
+       I[i+1]= (long)Q_primitive_part(idealmul(nf,a,b), &p1);
+       if (p1) A[i+1] = (long)element_mulvec(nf, p1,(GEN)A[i+1]);
+     }
    }
-   tetpil=avma; p1=cgetg(lg(order),t_VEC);
+   l = lg(order);
-   p1[1]=lcopy(A); p1[2]=lcopy(I);
+   p1 = cgetg(l,t_VEC);
-   for (i=3; i<lg(order); i++) p1[i]=lcopy((GEN)order[i]);
+   p1[1]=(long)A;
-   return gerepile(av,tetpil,p1);
+   p1[2]=(long)I; for (i=3; i<l; i++) p1[i]=order[i];
+   return gerepilecopy(av, p1);
  }
- /* Given bnf as output by buchinit and either an order as output by
+ /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
-  * rnfpseudobasis or a polynomial, and outputs a basis if it is free,
+  * and outputs a basis if it is free, an n+1-generating set if it is not */
-  * an n+1-generating set if it is not
-  */
  GEN
  rnfbasis(GEN bnf, GEN order)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long j,N,n;
+   long j, n;
-   GEN nf,A,I,classe,p1,p2,id;
+   GEN nf, A, I, cl, col, a, id;
-   bnf = checkbnf(bnf);
+   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
-   nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
+   id = idmat(degpol(nf[1]));
-   if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
+   order = get_order(nf, order, "rnfbasis");
-   if (typ(order)!=t_VEC || lg(order)<3)
+   I = (GEN)order[2]; n = lg(I)-1;
-     err(talker,"not a pseudo-matrix in rnfbasis");
-   A=(GEN)order[1]; I=(GEN)order[2]; n=lg(A)-1;
    j=1; while (j<n && gegal((GEN)I[j],id)) j++;
-   if (j<n) order=rnfsteinitz(nf,order);
+   if (j<n)
-   A=(GEN)order[1]; I=(GEN)order[2]; classe=(GEN)I[n];
-   p1=isprincipalgen(bnf,classe);
-   if (gcmp0((GEN)p1[1]))
    {
-     p2=cgetg(n+1,t_MAT);
+     order = rnfsteinitz(nf,order);
-     p2[n]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
+     I = (GEN)order[2];
    }
-   else
+   A = (GEN)order[1];
+   col= (GEN)A[n]; A = vecextract_i(A, 1, n-1);
+   cl = (GEN)I[n];
+   a = gen_if_principal(bnf, cl);
+   if (!a)
    {
-     p1=ideal_two_elt(nf,classe);
+     GEN p1 = ideal_two_elt(nf, cl);
-     p2=cgetg(n+2,t_MAT);
+     A = concatsp(A, gmul((GEN)p1[1], col));
-     p2[n]=lmul((GEN)p1[1],(GEN)A[n]);
+     a = (GEN)p1[2];
-     p2[n+1]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[n]);
    }
-   for (j=1; j<n; j++) p2[j]=A[j];
+   A = concatsp(A, element_mulvec(nf, a, col));
-   return gerepilecopy(av,p2);
+   return gerepilecopy(av, A);
  }
- /* Given bnf as output by buchinit and either an order as output by
+ /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
-  * rnfpseudobasis or a polynomial, and outputs a basis (not pseudo)
+  * and outputs a basis (not pseudo) in Hermite Normal Form if it exists, zero
-  * in Hermite Normal Form if it exists, zero if not
+  * if not
   */
  GEN
  rnfhermitebasis(GEN bnf, GEN order)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long j,N,n;
+   long j, n;
-   GEN nf,A,I,p1,id;
+   GEN nf, A, I, a, id;
-   bnf = checkbnf(bnf); nf=(GEN)bnf[7];
+   bnf = checkbnf(bnf); nf = (GEN)bnf[7];
-   N=degpol(nf[1]); id=idmat(N);
+   id = idmat(degpol(nf[1]));
-   if (typ(order)==t_POL)
+   order = get_order(nf, order, "rnfbasis");
-   {
+   A = (GEN)order[1]; A = dummycopy(A);
-     order=rnfpseudobasis(nf,order);
+   I = (GEN)order[2]; n = lg(A)-1;
-     A=(GEN)order[1];
-   }
-   else
-   {
-     if (typ(order)!=t_VEC || lg(order)<3)
-       err(talker,"not a pseudo-matrix in rnfbasis");
-     A=gcopy((GEN)order[1]);
-   }
-   I=(GEN)order[2]; n=lg(A)-1;
    for (j=1; j<=n; j++)
    {
-     if (!gegal((GEN)I[j],id))
+     if (gegal((GEN)I[j],id)) continue;
-     {
-       p1=isprincipalgen(bnf,(GEN)I[j]);
+     a = gen_if_principal(bnf, (GEN)I[j]);
-       if (gcmp0((GEN)p1[1]))
+     if (!a) { avma = av; return gzero; }
-         A[j]=(long)element_mulvec(nf,(GEN)p1[2],(GEN)A[j]);
+     A[j] = (long)element_mulvec(nf, a, (GEN)A[j]);
-       else { avma=av; return gzero; }
-     }
    }
    return gerepilecopy(av,A);
  }
- long
+ static long
- rnfisfree(GEN bnf, GEN order)
+ _rnfisfree(GEN bnf, GEN order)
  {
-   long av=avma,n,N,j;
+   long n, j;
-   GEN nf,p1,id,I;
+   GEN nf, p1, id, I;
    bnf = checkbnf(bnf);
    if (gcmp1(gmael3(bnf,8,1,1))) return 1;
-   nf=(GEN)bnf[7]; N=degpol(nf[1]); id=idmat(N);
+   nf = (GEN)bnf[7]; id = idmat(degpol(nf[1]));
-   if (typ(order)==t_POL) order=rnfpseudobasis(nf,order);
+   order = get_order(nf, order, "rnfisfree");
-   if (typ(order)!=t_VEC || lg(order)<3)
+   I = (GEN)order[2]; n = lg(I)-1;
-     err(talker,"not a pseudo-matrix in rnfisfree");
-   I=(GEN)order[2]; n=lg(I)-1;
    j=1; while (j<=n && gegal((GEN)I[j],id)) j++;
-   if (j>n) { avma=av; return 1; }
+   if (j>n) return 1;
-   p1=(GEN)I[j];
+   p1 = (GEN)I[j];
    for (j++; j<=n; j++)
-     if (!gegal((GEN)I[j],id)) p1=idealmul(nf,p1,(GEN)I[j]);
+     if (!gegal((GEN)I[j],id)) p1 = idealmul(nf,p1,(GEN)I[j]);
-   j = gcmp0(isprincipal(bnf,p1));
+   return gcmp0( isprincipal(bnf,p1) );
-   avma=av; return j;
  }
+ long
+ rnfisfree(GEN bnf, GEN order)
+ {
+   gpmem_t av = avma;
+   long n = _rnfisfree(bnf, order);
+   avma = av; return n;
+ }
  /**********************************************************************/
  /**                                                                  **/
  /**                   COMPOSITUM OF TWO NUMBER FIELDS                **/
  /**                                                                  **/
  /**********************************************************************/
- extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
+ /* modular version */
- extern GEN squff2(GEN x, long klim, long hint);
- extern GEN to_polmod(GEN x, GEN mod);
- /* modular version. TODO: check that compositum2 is not slower */
  GEN
  polcompositum0(GEN A, GEN B, long flall)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long v,k;
+   long v, k;
    GEN C, LPRS;
    if (typ(A)!=t_POL || typ(B)!=t_POL) err(typeer,"polcompositum0");
-Line 3034  polcompositum0(GEN A, GEN B, long flall)
+Line 3475  polcompositum0(GEN A, GEN B, long flall)
 Line 3034  polcompositum0(GEN A, GEN B, long flall)
 Line 3475  polcompositum0(GEN A, GEN B, long flall)
    if (!ZX_is_squarefree(B)) err(talker,"compositum: %Z not separable", B);
    k = 1; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
-   C = squff2(C,0,0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
+   C = DDF2(C, 0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
+   C = sort_vecpol(C);
    if (flall)
    {
      long i,l = lg(C);
      GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
      for (i=1; i<l; i++)
      { /* invmod possibly very costly */
-       a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], (GEN)C[i]));
+       a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], (GEN)C[i]));
        a = gneg_i(gmod(a, (GEN)C[i]));
        b = gadd(polx[v], gmulsg(k,a));
        w = cgetg(5,t_VEC); /* [C, a, b, n ] */
        w[1] = C[i];
        w[2] = (long)to_polmod(a, (GEN)w[1]);
        w[3] = (long)to_polmod(b, (GEN)w[1]);
        w[4] = lstoi(-k); C[i] = (long)w;
-Line 3066  compositum2(GEN pol1,GEN pol2)
+Line 3508  compositum2(GEN pol1,GEN pol2)
 Line 3066  compositum2(GEN pol1,GEN pol2)
 Line 3508  compositum2(GEN pol1,GEN pol2)
    return polcompositum0(pol1,pol2,1);
  }
- extern int isrational(GEN x);
- extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
  int
  nfissquarefree(GEN nf, GEN x)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
    GEN g, y = derivpol(x);
    if (isrational(x))
      g = modulargcd(x, y);
-Line 3082  nfissquarefree(GEN nf, GEN x)
+Line 3521  nfissquarefree(GEN nf, GEN x)
 Line 3082  nfissquarefree(GEN nf, GEN x)
 Line 3521  nfissquarefree(GEN nf, GEN x)
  }
  GEN
- rnfequation0(GEN nf, GEN B, long flall)
+ _rnfequation(GEN A, GEN B, long *pk, GEN *pLPRS)
  {
-   ulong av = avma;
+   long k, lA, lB;
-   long v,vpol,k,lA,lB;
+   GEN nf, C;
-   GEN cC,A,C,LPRS;
-   if (typ(nf)==t_POL) A=nf; else { nf=checknf(nf); A=(GEN)nf[1]; }
+   A = get_nfpol(A, &nf);
-   B = fix_relative_pol(nf,B,1);
+   B = fix_relative_pol(A,B,1);
-   v   = varn(A); lA = lgef(A);
+   lA = lgef(A);
-   vpol= varn(B); lB = lgef(B);
+   lB = lgef(B);
    if (lA<=3 || lB<=3) err(constpoler,"rnfequation");
    check_pol_int(A,"rnfequation");
-   B = lift_intern(B); B = gdiv(B, content(B));
+   B = primpart(lift_intern(B));
    for (k=2; k<lB; k++)
      if (lgef(B[k]) >= lA) B[k] = lres((GEN)B[k],A);
    if (!nfissquarefree(A,B))
      err(talker,"not k separable relative equation in rnfequation");
-   k = 0; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
+   *pk = 0; C = ZY_ZXY_resultant_all(A, B, pk, pLPRS);
    if (gsigne(leadingcoeff(C)) < 0) C = gneg_i(C);
-   C = primitive_part(C, &cC);
+   *pk = -*pk; return primpart(C);
+ }
+ GEN
+ rnfequation0(GEN A, GEN B, long flall)
+ {
+   gpmem_t av = avma;
+   long k;
+   GEN LPRS, nf, C;
+   A = get_nfpol(A, &nf);
+   C = _rnfequation(A, B, &k, flall? &LPRS: NULL);
    if (flall)
    {
-     GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
+     GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b + k a */
      /* invmod possibly very costly */
-     a = gmul((GEN)LPRS[1], ZX_invmod((GEN)LPRS[2], C));
+     a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], C));
      a = gneg_i(gmod(a, C));
-     b = gadd(polx[v], gmulsg(k,a));
+     b = gadd(polx[varn(A)], gmulsg(k,a));
      w = cgetg(4,t_VEC); /* [C, a, n ] */
      w[1] = (long)C;
      w[2] = (long)to_polmod(a, (GEN)w[1]);
-     w[3] = lstoi(-k); C = w;
+     w[3] = lstoi(k); C = w;
    }
    return gerepilecopy(av, C);
  }
-Line 3186  static GEN
+Line 3635  static GEN
 Line 3186  static GEN
 Line 3635  static GEN
  rnfdiv(GEN x, GEN y)
  {
    long i, ru = lg(x);
-   GEN z;
+   GEN z = cgetg(ru,t_COL);
-   z=cgetg(ru,t_COL);
    for (i=1; i<ru; i++) z[i]=(long)gdiv((GEN)x[i],(GEN)y[i]);
    return z;
  }
-Line 3197  static GEN
+Line 3645  static GEN
 Line 3197  static GEN
 Line 3645  static GEN
  rnfmul(GEN x, GEN y)
  {
    long i, ru = lg(x);
-   GEN z;
+   GEN z = cgetg(ru,t_COL);
-   z=cgetg(ru,t_COL);
    for (i=1; i<ru; i++) z[i]=(long)gmul((GEN)x[i],(GEN)y[i]);
    return z;
  }
-Line 3208  static GEN
+Line 3655  static GEN
 Line 3208  static GEN
 Line 3655  static GEN
  rnfvecmul(GEN x, GEN v)
  {
    long i, lx = lg(v);
-   GEN y;
+   GEN y = cgetg(lx,typ(v));
-   y=cgetg(lx,typ(v));
    for (i=1; i<lx; i++) y[i]=(long)rnfmul(x,(GEN)v[i]);
    return y;
  }
  static GEN
- allonge(GEN v, long N)
+ allonge(GEN v, long l)
  {
-   long r,r2,i;
+   long r, r2, i;
-   GEN y;
+   GEN y = cgetg(l,t_COL);
-   r=lg(v)-1; r2=N-r;
+   r = lg(v); r2 = l-r;
-   y=cgetg(N+1,t_COL);
+   for (i=1; i < r; i++) y[i] = v[i];
-   for (i=1; i<=r; i++) y[i]=v[i];
+   for (   ; i < l; i++) y[i] = lconj((GEN)v[i-r2]);
-   for ( ; i<=N; i++) y[i]=(long)gconj((GEN)v[i-r2]);
    return y;
  }
  static GEN
  findmin(GEN nf, GEN ideal, GEN muf,long prec)
  {
-   long av=avma,N,tetpil,i;
+   gpmem_t av = avma;
-   GEN m,y;
+   long i, l;
+   GEN m,y, G = gmael(nf,5,2);
-   m = qf_base_change(gmael(nf,5,3), ideal, 0); /* nf[5][3] = T2 */
+   m = gmul(G, ideal);
-   m = lllgramintern(m,4,1,prec);
+   m = lllintern(m,4,1,prec);
    if (!m)
    {
-     m = lllint(ideal);
+     ideal = lllint_ip(ideal,4);
-     m = qf_base_change(gmael(nf,5,3), gmul(ideal,m), 0);
+     m = gmul(G, ideal);
-     m = lllgramintern(m,4,1,prec);
+     m = lllintern(m,4,1,prec);
      if (!m) err(precer,"rnflllgram");
    }
-   ideal=gmul(ideal,m);
+   ideal = gmul(ideal,m);
-   N=lg(ideal)-1; y=cgetg(N+1,t_MAT);
+   l = lg(ideal); y = cgetg(l,t_MAT);
-   for (i=1; i<=N; i++)
+   for (i=1; i<l; i++)
-     y[i] = (long) allonge(nftocomplex(nf,(GEN)ideal[i]),N);
+     y[i] = (long)allonge(nftocomplex(nf,(GEN)ideal[i]),l);
-   m=ground(greal(gauss(y,allonge(muf,N))));
+   m = ground(greal(gauss(y, allonge(muf,l))));
-   tetpil=avma; return gerepile(av,tetpil,gmul(ideal,m));
+   return gerepileupto(av, gmul(ideal,m));
  }
  #define swap(x,y) { long _t=x; x=y; y=_t; }
-Line 3260  findmin(GEN nf, GEN ideal, GEN muf,long prec)
+Line 3706  findmin(GEN nf, GEN ideal, GEN muf,long prec)
 Line 3260  findmin(GEN nf, GEN ideal, GEN muf,long prec)
 Line 3706  findmin(GEN nf, GEN ideal, GEN muf,long prec)
  GEN
  rnflllgram(GEN nf, GEN pol, GEN order,long prec)
  {
-   long av=avma,tetpil,i,j,k,l,kk,kmax,r1,ru,lx,vnf;
+   gpmem_t av=avma,tetpil;
+   long i,j,k,l,kk,kmax,r1,ru,lx,vnf;
    GEN p1,p2,M,I,U,ronf,poll,unro,roorder,powreorder,mth,s,MC,MPOL,MCS;
    GEN B,mu,Bf,temp,ideal,x,xc,xpol,muf,mufc,muno,y,z,Ikk_inv;
-Line 3431  rnflllgram(GEN nf, GEN pol, GEN order,long prec)
+Line 3878  rnflllgram(GEN nf, GEN pol, GEN order,long prec)
 Line 3431  rnflllgram(GEN nf, GEN pol, GEN order,long prec)
 Line 3878  rnflllgram(GEN nf, GEN pol, GEN order,long prec)
  GEN
  rnfpolred(GEN nf, GEN pol, long prec)
  {
-   ulong av = avma;
+   gpmem_t av = avma;
-   long i,j,k,n,N, vpol = varn(pol);
+   long i, j, n, v = varn(pol);
-   GEN id,id2,newid,newor,p1,p2,al,newpol,w,z;
+   GEN id, al, w, I, O, bnf;
-   GEN bnf,zk,newideals,ideals,order,neworder;
    if (typ(pol)!=t_POL) err(typeer,"rnfpolred");
-   if (typ(nf)!=t_VEC) err(idealer1);
+   bnf = nf; nf = checknf(bnf);
-   switch(lg(nf))
+   bnf = (nf == bnf)? NULL: checkbnf(bnf);
-   {
+   if (degpol(pol) <= 1) return _vec(polx[v]);
-     case 10: bnf = NULL; break;
-     case 11: bnf = nf; nf = checknf((GEN)nf[7]); break;
+   id = rnfpseudobasis(nf,pol);
-     default: err(idealer1);
-       return NULL; /* not reached */
-   }
-   if (degpol(pol) <= 1)
-   {
-     w=cgetg(2,t_VEC);
-     w[1]=lpolx[vpol]; return w;
-   }
-   id=rnfpseudobasis(nf,pol); N=degpol(nf[1]);
    if (bnf && gcmp1(gmael3(bnf,8,1,1))) /* if bnf is principal */
    {
-     ideals=(GEN)id[2]; n=lg(ideals)-1; order=(GEN)id[1];
+     GEN newI, newO, zk = idmat(degpol(nf[1]));
-     newideals=cgetg(n+1,t_VEC); neworder=cgetg(n+1,t_MAT);
+     O = (GEN)id[1];
-     zk=idmat(N);
+     I = (GEN)id[2]; n = lg(I)-1;
+     newI = cgetg(n+1,t_VEC);
+     newO = cgetg(n+1,t_MAT);
      for (j=1; j<=n; j++)
      {
-       newideals[j]=(long)zk; p1=cgetg(n+1,t_COL); neworder[j]=(long)p1;
+       newI[j] = (long)zk; al = gen_if_principal(bnf,(GEN)I[j]);
-       p2=(GEN)order[j];
+       newO[j] = (long)element_mulvec(nf, al, (GEN)O[j]);
-       al=(GEN)isprincipalgen(bnf,(GEN)ideals[j])[2];
-       for (k=1; k<=n; k++)
-         p1[k]=(long)element_mul(nf,(GEN)p2[k],al);
      }
-     id=cgetg(3,t_VEC); id[1]=(long)neworder; id[2]=(long)newideals;
+     id = cgetg(3,t_VEC);
+     id[1] = (long)newO;
+     id[2] = (long)newI;
    }
-   id2=rnflllgram(nf,pol,id,prec);
-   z=(GEN)id2[1]; newid=(GEN)z[2]; newor=(GEN)z[1];
+   id = (GEN)rnflllgram(nf,pol,id,prec)[1];
-   n=lg(newor)-1; w=cgetg(n+1,t_VEC);
+   O = (GEN)id[1];
+   I = (GEN)id[2]; n = lg(I)-1;
+   w = cgetg(n+1,t_VEC);
+   pol = lift(pol);
    for (j=1; j<=n; j++)
    {
-     p1=(GEN)newid[j]; al=gmul(gcoeff(p1,1,1),(GEN)newor[j]);
+     GEN p1, newpol;
-     p1=basistoalg(nf,(GEN)al[n]);
+     p1 = (GEN)I[j]; al = gmul(gcoeff(p1,1,1),(GEN)O[j]);
+     p1 = basistoalg(nf,(GEN)al[n]);
      for (i=n-1; i; i--)
-       p1=gadd(basistoalg(nf,(GEN)al[i]),gmul(polx[vpol],p1));
+       p1 = gadd(basistoalg(nf,(GEN)al[i]), gmul(polx[v],p1));
-     newpol=gtopoly(gmodulcp(gtovec(caract2(lift(pol),lift(p1),vpol)),
+     p1 = gmodulcp(gtovec(caract2(pol,lift(p1),v)), (GEN)nf[1]);
-                             (GEN) nf[1]), vpol);
+     newpol = gtopoly(p1, v);
      p1 = ggcd(newpol, derivpol(newpol));
-     if (degpol(p1)>0)
+     if (degpol(p1) > 0)
      {
-       newpol=gdiv(newpol,p1);
+       newpol = gdiv(newpol, p1);
-       newpol=gdiv(newpol,leading_term(newpol));
+       newpol = gdiv(newpol, leading_term(newpol));
      }
-     w[j]=(long)newpol;
+     w[j] = (long)newpol;
-     if (DEBUGLEVEL>=4) outerr(newpol);
    }
    return gerepilecopy(av,w);
  }
- extern GEN vecpol_to_mat(GEN v, long n);
  /* given a relative polynomial pol over nf, compute a pseudo-basis for the
   * extension, then an absolute basis */
- GEN
+ static GEN
- makebasis(GEN nf,GEN pol)
+ makebasis(GEN nf, GEN pol, GEN rnfeq)
  {
-   GEN elts,ids,polabs,plg,B,bs,p1,p2,a,den,vbs,vbspro,vpro,rnf;
+   GEN elts,ids,polabs,plg,plg0,B,bs,p1,den,vbs,d,vpro;
-   long av=avma,n,N,m,i,j, v = varn(pol);
+   gpmem_t av = avma;
+   long n,N,m,i,j,k, v = varn(pol);
-   p1 = rnfequation2(nf,pol);
+   polabs= (GEN)rnfeq[1];
-   polabs= (GEN)p1[1];
+   plg   = (GEN)rnfeq[2]; plg = lift_intern(plg);
-   plg   = (GEN)p1[2];
+   p1 = rnfpseudobasis(nf,pol);
-   a     = (GEN)p1[3];
-   rnf = cgetg(12,t_VEC);
-   for (i=2;i<=9;i++) rnf[i]=zero;
-   rnf[1] =(long)pol;
-   rnf[10]=(long)nf; p2=cgetg(4,t_VEC);
-   rnf[11]=(long)p2; p2[1]=p2[2]=zero; p2[3]=(long)a;
-   if (signe(a))
-     pol = gsubst(pol,v,gsub(polx[v],
-                             gmul(a,gmodulcp(polx[varn(nf[1])],(GEN)nf[1]))));
-   p1=rnfpseudobasis(nf,pol);
    elts= (GEN)p1[1];
    ids = (GEN)p1[2];
-   if (DEBUGLEVEL>1) { fprintferr("relative basis computed\n"); flusherr(); }
+   if (DEBUGLEVEL>1) fprintferr("relative basis computed\n");
-   N=degpol(pol); n=degpol((GEN)nf[1]); m=n*N;
+   N = degpol(pol);
-   den = denom(content(lift(plg)));
+   n = degpol(nf[1]); m = n*N;
-   vbs = cgetg(n+1,t_VEC);
-   vbs[1] = un;
+   plg0= Q_remove_denom(plg, &den); /* plg = plg0/den */
-   vbs[2] = (long)plg; vbspro = gmul(den,plg);
+   /* nf = K = Q(a), vbs[i+1] = a^i as an elt of L = Q[X] / polabs */
-   for(i=3;i<=n;i++)
+   vbs = RXQ_powers(plg0, polabs, n-1);
-     vbs[i] = ldiv(gmul((GEN)vbs[i-1],vbspro),den);
+   if (den)
+   { /* restore denominators */
+     vbs[2] = (long)plg; d = den;
+     for (i=3; i<=n; i++) { d = mulii(d,den); vbs[i] = ldiv((GEN)vbs[i], d); }
+   }
+   /* bs = integer basis of K, as elements of L */
    bs = gmul(vbs, vecpol_to_mat((GEN)nf[7],n));
-   vpro=cgetg(N+1,t_VEC);
+   vpro = cgetg(N+1,t_VEC);
-   for (i=1;i<=N;i++)
+   for (i=1; i<=N; i++) vpro[i] = lpowgs(polx[v], i-1);
+   vpro = gmul(vpro, elts);
+   B = cgetg(m+1, t_MAT);
+   for(i=k=1; i<=N; i++)
    {
-     p1=cgetg(3,t_POLMOD);
+     GEN w = element_mulvec(nf, (GEN)vpro[i], (GEN)ids[i]);
-     p1[1]=(long)polabs;
+     for(j=1; j<=n; j++)
-     p1[2]=lpuigs(polx[v],i-1); vpro[i]=(long)p1;
+     {
+       p1 = gres(gmul(bs, (GEN)w[j]), polabs);
+       B[k++] = (long)pol_to_vec(p1, m);
+     }
    }
-   vpro=gmul(vpro,elts); B = cgetg(m+1, t_MAT);
+   B = Q_remove_denom(B, &den);
-   for(i=1;i<=N;i++)
+   if (den) { B = hnfmodid(B, den); B = gdiv(B, den); } else B = idmat(m);
-     for(j=1;j<=n;j++)
+   p1 = cgetg(3,t_VEC);
+   p1[1] = (long)polabs;
+   p1[2] = (long)B; return gerepilecopy(av, p1);
+ }
+ /* relative polredabs. Returns relative polynomial by default (flag = 0)
+  * flag & nf_ORIG: + element (base change)
+  * flag & nf_ADDZK: + integer basis
+  * flag & nf_ABSOLUTE: absolute polynomial */
+ GEN
+ rnfpolredabs(GEN nf, GEN relpol, long flag)
+ {
+   GEN red, bas, z, elt, POL, pol, T, a;
+   long v, fl;
+   gpmem_t av = avma;
+   if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
+   nf = checknf(nf); v = varn(relpol);
+   if (DEBUGLEVEL>1) (void)timer2();
+   relpol = unifpol(nf,relpol,1);
+   T = (GEN)nf[1];
+   if ((flag & nf_ADDZK) && (flag != (nf_ADDZK|nf_ABSOLUTE)))
+     err(impl,"this combination of flags in rnfpolredabs");
+   fl = (flag & nf_ADDZK)? nf_ORIG: nf_RAW;
+   if (flag & nf_PARTIALFACT)
+   {
+     long sa;
+     bas = NULL; /* -Wall */
+     POL = _rnfequation(nf, relpol, &sa, NULL);
+     red = polredabs0(POL, fl | nf_PARTIALFACT);
+     a = stoi(sa);
+   }
+   else
+   {
+     GEN rel, eq = rnfequation2(nf,relpol);
+     POL = (GEN)eq[1];
+     a   = (GEN)eq[3];
+     rel = poleval(relpol, gsub(polx[v],
+                                gmul(a, gmodulcp(polx[varn(T)],T))));
+     bas = makebasis(nf, rel, eq);
+     if (DEBUGLEVEL>1)
      {
-       p1 = gmul(bs, element_mul(nf,(GEN)vpro[i],gmael(ids,i,j)));
+       msgtimer("absolute basis");
-       B[(i-1)*n+j] = (long)pol_to_vec(lift_intern(p1), m);
+       fprintferr("original absolute generator: %Z\n", POL);
      }
-   p1 = denom(B); B = gmul(B,p1);
+     red = polredabs0(bas, fl);
-   B = hnfmodid(B, p1); B = gdiv(B,p1);
+   }
-   p1=cgetg(4,t_VEC);
+   pol = (GEN)red[1];
-   p1[1]=(long)polabs;
+   if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",pol);
-   p1[2]=(long)B;
+   if (flag & nf_ABSOLUTE)
-   p1[3]=(long)rnf; return gerepilecopy(av, p1);
+   {
+     if (flag & nf_ADDZK)
+     {
+       GEN t = (GEN)red[2], B = (GEN)bas[2];
+       GEN v = RXQ_powers(lift_intern(t), pol, degpol(pol)-1);
+       z = cgetg(3, t_VEC);
+       z[1] = (long)pol;
+       z[2] = lmul(v, B);
+       return gerepilecopy(av, z);
+     }
+     return gerepilecopy(av, pol);
+   }
+   elt = eltabstorel((GEN)red[2], T, relpol, a);
+   pol = rnfcharpoly(nf,relpol,elt,v);
+   if (!(flag & nf_ORIG)) return gerepileupto(av, pol);
+   z = cgetg(3,t_VEC);
+   z[1] = (long)pol;
+   z[2] = (long)to_polmod(modreverse_i((GEN)elt[2],(GEN)elt[1]), pol);
+   return gerepilecopy(av, z);
  }

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