[BACK]Return to nffactor.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / modules

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

version 1.1, 2001/10/02 11:17:11 version 1.2, 2002/09/11 07:27:06
Line 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"  #include "pari.h"
 #include "parinf.h"  #include "parinf.h"
   
   extern void init_dalloc();
   extern double *dalloc(size_t n);
   extern GEN gmul_mati_smallvec(GEN x, GEN y);
   extern GEN GS_norms(GEN B, long prec);
   extern GEN RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr);
   extern GEN RXQX_red(GEN P, GEN T);
   extern GEN apply_kummer(GEN nf,GEN pol,GEN e,GEN p);
   extern GEN centermod_i(GEN x, GEN p, GEN ps2);
   extern GEN dim1proj(GEN prh);
 extern GEN hensel_lift_fact(GEN pol, GEN fact, GEN T, GEN p, GEN pev, long e);  extern GEN hensel_lift_fact(GEN pol, GEN fact, GEN T, GEN p, GEN pev, long e);
 extern GEN nf_get_T2(GEN base, GEN polr);  extern GEN initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis);
 extern GEN nfreducemodpr_i(GEN x, GEN prh);  extern GEN max_modulus(GEN p, double tau);
   extern GEN mulmat_pol(GEN A, GEN x);
   extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
   extern GEN polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N);
 extern GEN sort_factor(GEN y, int (*cmp)(GEN,GEN));  extern GEN sort_factor(GEN y, int (*cmp)(GEN,GEN));
 extern GEN pidealprimeinv(GEN nf, GEN x);  extern GEN special_pivot(GEN x);
   extern GEN trivfact(void);
   extern GEN vandermondeinverse(GEN L, GEN T, GEN den, GEN prep);
   extern GEN vconcat(GEN A, GEN B);
   extern int cmbf_precs(GEN q, GEN A, GEN B, long *a, long *b, GEN *qa, GEN *qb);
   extern int isrational(GEN x);
   extern GEN LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, pari_timer *T, long *ti_LLL);
   extern void remake_GM(GEN nf, nffp_t *F, long prec);
   #define RXQX_div(x,y,T) RXQX_divrem((x),(y),(T),NULL)
   #define RXQX_rem(x,y,T) RXQX_divrem((x),(y),(T),ONLY_REM)
   
 static GEN nffactormod2(GEN nf,GEN pol,GEN pr);  
 static GEN nfmod_split2(GEN nf, GEN prhall, GEN polb, GEN v, GEN q);  
 static GEN nf_pol_mul(GEN nf,GEN pol1,GEN pol2);  
 static GEN nf_pol_divres(GEN nf,GEN pol1,GEN pol2, GEN *pr);  
 static GEN nf_pol_subres(GEN nf,GEN pol1,GEN pol2);  
 static GEN nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol);  
 static GEN nfmod_pol_divres(GEN nf,GEN prhall,GEN pol1,GEN pol2, GEN *pr);  
 static GEN nfmod_pol_gcd(GEN nf,GEN prhall,GEN pol1,GEN pol2);  
 static GEN nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt);  
 static GEN nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol);  
 static GEN nfsqff(GEN nf,GEN pol,long fl);  static GEN nfsqff(GEN nf,GEN pol,long fl);
 static int nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint);  
   
 typedef struct Nfcmbf /* for use in nfsqff */  /* for nf_bestlift: reconstruction of algebraic integers known mod \wp^k */
   /* FIXME: assume \wp has degree 1 for now */
   typedef struct {
     long k;    /* input known mod \wp^k */
     GEN pk;    /* p^k = Norm(\wp^k)*/
     GEN den;   /* denom(prk^-1) = p^k */
     GEN prk;   /* |.|^2 LLL-reduced basis (b_i) of \wp^k  (NOT T2-reduced) */
     GEN iprk;  /* den * prk^-1 */
     GEN GSmin; /* min |b_i^*|^2 */
     GEN ZpProj;/* projector to Zp / \wp^k = Z/p^k  (\wp unramified, degree 1) */
     GEN tozk;
     GEN topow;
   } nflift_t;
   
   typedef struct /* for use in nfsqff */
 {  {
   GEN pol, h, hinv, fact, res, lt, den;    nflift_t *L;
   long nfact, nfactmod;    GEN nf, pol, fact, res, bound, pr, pk, Br, ZC, dn, polbase, BS_2;
 } Nfcmbf;    long hint;
 static Nfcmbf nfcmbf;  } nfcmbf_t;
   
   GEN
   RXQX_mul(GEN x, GEN y, GEN T)
   {
     return RXQX_red(gmul(x,y), T);
   }
   
 static GEN  static GEN
 unifpol0(GEN nf,GEN pol,long flag)  unifpol0(GEN nf,GEN pol,long flag)
 {  {
   static long n = 0;    static long n = 0;
   static GEN vun = NULL;    static GEN vun = NULL;
   GEN f = (GEN) nf[1];    GEN f = (GEN) nf[1];
   long i = degpol(f), av;    long i = degpol(f);
     gpmem_t av;
   
   if (i != n)    if (i != n)
   {    {
Line 67  unifpol0(GEN nf,GEN pol,long flag)
Line 98  unifpol0(GEN nf,GEN pol,long flag)
   switch(typ(pol))    switch(typ(pol))
   {    {
     case t_INT: case t_FRAC: case t_RFRAC:      case t_INT: case t_FRAC: case t_RFRAC:
         if (flag) return gcopy(pol);
       pol = gmul(pol,vun); break;        pol = gmul(pol,vun); break;
   
     case t_POL:      case t_POL:
         if (flag && !degpol(pol)) return gcopy(constant_term(pol));
       pol = gmodulcp(pol,f); /* fall through */        pol = gmodulcp(pol,f); /* fall through */
     case t_POLMOD:      case t_POLMOD:
       pol = algtobasis(nf,pol);        pol = algtobasis(nf,pol);
Line 100  unifpol(GEN nf,GEN pol,long flag)
Line 133  unifpol(GEN nf,GEN pol,long flag)
   return unifpol0(nf,(GEN) pol, flag);    return unifpol0(nf,(GEN) pol, flag);
 }  }
   
 #if 0  /* factorization of x modulo pr. Assume x already reduced */
 /* return a monic polynomial of degree d with random coefficients in Z_nf */  GEN
 static GEN  FqX_factor(GEN x, GEN T, GEN p)
 random_pol(GEN nf,long d)  
 {  {
   long i,j, n = degpol(nf[1]);    GEN rep;
   GEN pl,p;    if (!T)
   
   pl=cgetg(d+3,t_POL);  
   for (i=2; i<d+2; i++)  
   {    {
     p=cgetg(n+1,t_COL); pl[i]=(long)p;      rep = factmod0(x, p);
     for (j=1; j<=n; j++)      rep[2] = (long)small_to_vec((GEN)rep[2]);
       p[j] = lstoi(mymyrand()%101 - 50);  
   }    }
   p=cgetg(n+1,t_COL); pl[i]=(long)p;    else
   p[1]=un; for (i=2; i<=n; i++) p[i]=zero;    {
       rep = factmod9(x, p, T);
   pl[1] = evalsigne(1) | evallgef(d+3) | evalvarn(0);      rep = lift_intern(lift_intern(rep));
   return pl;    }
     return rep;
 }  }
 #endif  
   
 /* multiplication of x by y */  GEN
 static GEN  nffactormod(GEN nf, GEN x, GEN pr)
 nf_pol_mul(GEN nf,GEN x,GEN y)  
 {  {
   long tetpil,av=avma;    long j, l, vx = varn(x), vn;
   GEN res = gmul(unifpol(nf,x,1), unifpol(nf,y,1));    gpmem_t av = avma;
     GEN z, rep, xrd, modpr, T, p;
   
   tetpil = avma;    nf = checknf(nf);
   return gerepile(av,tetpil,unifpol(nf,res,0));    vn = varn((GEN)nf[1]);
     if (typ(x)!=t_POL) err(typeer,"nffactormod");
     if (vx >= vn)
       err(talker,"polynomial variable must have highest priority in nffactormod");
   
     modpr = nf_to_ff_init(nf, &pr, &T, &p);
     xrd = modprX(x, nf, modpr);
     rep = FqX_factor(xrd,T,p);
     z = (GEN)rep[1]; l = lg(z);
     for (j = 1; j < l; j++) z[j] = (long)modprX_lift((GEN)z[j], modpr);
     return gerepilecopy(av, rep);
 }  }
   
 /* compute x^2 in nf */  /* If p doesn't divide bad and has a divisor of degree 1, return the latter. */
 static GEN  static GEN
 nf_pol_sqr(GEN nf,GEN x)  choose_prime(GEN nf, GEN bad, GEN *p, byteptr *PT)
 {  {
   long tetpil,av=avma;    GEN  q = icopy(gun), r, x = (GEN)nf[1];
   GEN res = gsqr(unifpol(nf,x,1));    ulong pp = *p? itou(*p): 0;
     byteptr pt = *PT;
   tetpil = avma;    gpmem_t av = avma;
   return gerepile(av,tetpil,unifpol(nf,res,0));    for (;;)
     {
       NEXT_PRIME_VIADIFF_CHECK(pp, pt);
       if (! smodis(bad,pp)) continue;
       affui(pp, q);
       r = rootmod(x, q); if (lg(r) > 1) break;
       avma = av;
     }
     r = gsub(polx[varn(x)], lift_intern((GEN)r[1]));
     *p = q;
     *PT = pt; return apply_kummer(nf,r,gun,q);
 }  }
   
 /* reduce the coefficients of pol modulo prhall */  
 static GEN  static GEN
 nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol)  QXQ_normalize(GEN P, GEN T)
 {  {
   long av=avma,tetpil,i;    GEN t = leading_term(P);
   GEN p1;    if (!gcmp1(t))
     {
   if (typ(pol)!=t_POL) return nfreducemodpr(nf,pol,prhall);      if (is_rational_t(typ(t)))
   pol=unifpol(nf,pol,0);        P = gdiv(P, t);
       else
   tetpil=avma; i=lgef(pol);        P = RXQX_mul(QX_invmod(t,T), P, T);
   p1=cgetg(i,t_POL); p1[1]=pol[1];    }
   for (i--; i>=2; i--)    return P;
     p1[i] = (long) nfreducemodpr(nf,(GEN)pol[i],prhall);  
   return gerepile(av,tetpil, normalizepol(p1));  
 }  }
   
 /* x^2 modulo prhall */  /* return the roots of pol in nf */
 static GEN  GEN
 nfmod_pol_sqr(GEN nf,GEN prhall,GEN x)  nfroots(GEN nf,GEN pol)
 {  {
   long av=avma,tetpil;    gpmem_t av = avma;
   GEN px;    int d = degpol(pol);
     GEN A,g, T;
   
   px = nfmod_pol_reduce(nf,prhall,x);    nf = checknf(nf); T = (GEN)nf[1];
   px = unifpol(nf,lift(px),1);    if (typ(pol) != t_POL) err(notpoler,"nfroots");
   px = unifpol(nf,nf_pol_sqr(nf,px),0);    if (varn(pol) >= varn(T))
   tetpil=avma;      err(talker,"polynomial variable must have highest priority in nfroots");
   return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));    if (d == 0) return cgetg(1,t_VEC);
     if (d == 1)
     {
       A = gneg_i(gdiv((GEN)pol[2],(GEN)pol[3]));
       return gerepilecopy(av, _vec( basistoalg(nf,A) ));
     }
     A = fix_relative_pol(nf,pol,0);
     A = primpart( lift_intern(A) );
     if (DEBUGLEVEL>3) fprintferr("test if polynomial is square-free\n");
     g = nfgcd(A, derivpol(A), T, NULL);
   
     if (degpol(g))
     { /* not squarefree */
       g = QXQ_normalize(g, T);
       A = RXQX_div(A,g,T);
     }
     A = QXQ_normalize(A, T);
     A = primpart(A);
     A = nfsqff(nf,A,1);
     return gerepileupto(av, gen_sort(A, 0, cmp_pol));
 }  }
   
 /* multiplication of x by y modulo prhall */  int
 static GEN  nfissplit(GEN nf, GEN x)
 nfmod_pol_mul(GEN nf,GEN prhall,GEN x,GEN y)  
 {  {
   long av=avma,tetpil;    gpmem_t av = avma;
   GEN px,py;    long l = lg(nfsqff(checknf(nf), x, 2));
     avma = av; return l != 1;
   px = nfmod_pol_reduce(nf,prhall,x); px = unifpol(nf,lift(px),1);  
   py = nfmod_pol_reduce(nf,prhall,y); py = unifpol(nf,lift(py),1);  
   px = unifpol(nf,nf_pol_mul(nf,px,py),0);  
   tetpil=avma;  
   return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));  
 }  }
   
 /* Euclidan division of x by y */  /* nf = K[y] / (P). Galois over K ? */
 static GEN  int
 nf_pol_divres(GEN nf,GEN x,GEN y,GEN *pr)  nfisgalois(GEN nf, GEN P)
 {  {
   long av = avma,tetpil;    return degpol(P) <= 2 || nfissplit(nf, P);
   GEN nq = poldivres(unifpol(nf,x,1),unifpol(nf,y,1),pr);  
   GEN *gptr[2];  
   
   tetpil=avma; nq=unifpol(nf,nq,0);  
   if (pr) *pr = unifpol(nf,*pr,0);  
   gptr[0]=&nq; gptr[1]=pr;  
   gerepilemanysp(av,tetpil,gptr,pr ? 2:1);  
   return nq;  
 }  }
   
 /* Euclidan division of x by y modulo prhall */  /* return a minimal lift of elt modulo id */
 static GEN  static GEN
 nfmod_pol_divres(GEN nf,GEN prhall,GEN x,GEN y, GEN *pr)  nf_bestlift(GEN elt, GEN bound, nflift_t *T)
 {  {
   long av=avma,dx,dy,dz,i,j,k,l,n,tetpil;    GEN u;
   GEN z,p1,p3,px,py;    long i,l = lg(T->prk), t = typ(elt);
     if (t != t_INT)
   py = nfmod_pol_reduce(nf,prhall,y);  
   if (gcmp0(py))  
     err(talker, "division by zero in nfmod_pol_divres");  
   
   tetpil=avma;  
   px=nfmod_pol_reduce(nf,prhall,x);  
   dx=degpol(px); dy=degpol(py); dz=dx-dy;  
   if (dx<dy)  
   {    {
     GEN vzero;      if (t == t_POL) elt = mulmat_pol(T->tozk, elt);
       u = gmul(T->iprk,elt);
     if (pr) *pr = gerepile(av,tetpil,px);      for (i=1; i<l; i++) u[i] = (long)diviiround((GEN)u[i], T->pk);
     else avma = av;  
   
     n=degpol(nf[1]);  
     vzero = cgetg(n+1,t_COL);  
     n=degpol(nf[1]);  
     for (i=1; i<=n; i++) vzero[i]=zero;  
   
     z=cgetg(3,t_POL); z[2]=(long)vzero;  
     z[1]=evallgef(2) | evalvarn(varn(px));  
     return z;  
   }    }
   p1 = NULL; /* gcc -Wall */    else
   
   z=cgetg(dz+3,t_POL); z[1]=evalsigne(1) | evallgef(3+dz);  
   setvarn(z,varn(px));  
   z[dz+2] = (long) element_divmodpr(nf,(GEN)px[dx+2],(GEN)py[dy+2],prhall);  
   for (i=dx-1; i>=dy; --i)  
   {    {
     l=avma; p1=(GEN)px[i+2];      u = gmul(elt, (GEN)T->iprk[1]);
     for (j=i-dy+1; j<=i && j<=dz; j++)      for (i=1; i<l; i++) u[i] = (long)diviiround((GEN)u[i], T->pk);
       p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));      elt = gscalcol(elt, l-1);
     p1 = nfreducemodpr(nf,p1,prhall);  
     tetpil=avma; p3=element_divmodpr(nf,p1,(GEN)py[dy+2],prhall);  
     z[i-dy+2]=lpile(l,tetpil,p3);  
     z[i-dy+2]=(long)nfreducemodpr(nf,(GEN)z[i-dy+2],prhall);  
   }    }
   l=avma;    u = gsub(elt, gmul(T->prk, u));
   for (i=dy-1; i>=0; --i)    if (bound && gcmp(QuickNormL2(u,DEFAULTPREC), bound) > 0) u = NULL;
   {    return u;
     l=avma; p1=((GEN)px[i+2]);  
     for (j=0; j<=i && j<=dz; j++)  
       p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));  
     p1 = gerepileupto(l, nfreducemodpr(nf,p1,prhall));  
     if (!gcmp0(p1)) break;  
   }  
   
   if (!pr) { avma = l; return z; }  
   
   if (i<0)  
   {  
     avma=l;  
     p3 = cgetg(3,t_POL); p3[2]=zero;  
     p3[1] = evallgef(2) | evalvarn(varn(px));  
     *pr=p3; return z;  
   }  
   
   p3=cgetg(i+3,t_POL);  
   p3[1]=evalsigne(1) | evallgef(i+3) | evalvarn(varn(px));  
   p3[i+2]=(long)nfreducemodpr(nf,p1,prhall);  
   for (k=i-1; k>=0; --k)  
   {  
     l=avma; p1=((GEN)px[k+2]);  
     for (j=0; j<=k && j<=dz; j++)  
       p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[k-j+2]));  
     p3[k+2]=lpileupto(l,nfreducemodpr(nf,p1,prhall));  
   }  
   *pr=p3; return z;  
 }  }
   
 /* GCD of x and y */  
 static GEN  static GEN
 nf_pol_subres(GEN nf,GEN x,GEN y)  nf_bestlift_to_pol(GEN elt, GEN bound, nflift_t *T)
 {  {
   long av=avma,tetpil;    GEN u = nf_bestlift(elt,bound,T);
   GEN s = srgcd(unifpol(nf,x,1), unifpol(nf,y,1));    if (u) u = gmul(T->topow, u);
     return u;
   tetpil=avma; return gerepile(av,tetpil,unifpol(nf,s,1));  
 }  }
   
 /* GCD of x and y modulo prhall */  /* return the lift of pol with coefficients of T2-norm <= C (if possible) */
 static GEN  static GEN
 nfmod_pol_gcd(GEN nf,GEN prhall,GEN x,GEN y)  nf_pol_lift(GEN pol, GEN bound, nfcmbf_t *T)
 {  {
   long av=avma;    long i, d = lgef(pol);
   GEN p1,p2;    GEN x = cgetg(d,t_POL);
     nflift_t *L = T->L;
   
   if (lgef(x)<lgef(y)) { p1=y; y=x; x=p1; }    x[1] = pol[1];
   p1=nfmod_pol_reduce(nf,prhall,x);    for (i=2; i<d; i++)
   p2=nfmod_pol_reduce(nf,prhall,y);  
   while (!isexactzero(p2))  
   {    {
     GEN p3;      x[i] = (long)nf_bestlift_to_pol((GEN)pol[i], bound, L);
       if (!x[i]) return NULL;
     nfmod_pol_divres(nf,prhall,p1,p2,&p3);  
     p1=p2; p2=p3;  
   }    }
   return gerepileupto(av,p1);    return x;
 }  }
   
 /* return pol^e modulo prhall and the polynomial pmod */  /* return the factorization of x in nf */
 static GEN  GEN
 nfmod_pol_pow(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN e)  nffactor(GEN nf,GEN pol)
 {  {
   long i, av = avma, n = degpol(nf[1]);    GEN A,g,y,p1,rep,T;
   GEN p1,p2,vun;    long l, j, d = degpol(pol);
     gpmem_t av;
     if (DEBUGLEVEL>3) (void)timer2();
   
   vun=cgetg(n+1,t_COL); vun[1]=un; for (i=2; i<=n; i++) vun[i]=zero;    nf = checknf(nf); T = (GEN)nf[1];
   p1=gcopy(polun[varn(pol)]); p1[2]=(long)vun;    if (typ(pol) != t_POL) err(notpoler,"nffactor");
   if (gcmp0(e)) return p1;    if (varn(pol) >= varn(T))
       err(talker,"polynomial variable must have highest priority in nffactor");
   
   p2=nfmod_pol_reduce(nf,prhall,pol);    if (d == 0) return trivfact();
   for(;;)    rep = cgetg(3, t_MAT); av = avma;
     if (d == 1)
   {    {
     if (!vali(e))      rep[1] = (long)_col( gcopy(pol) );
     {      rep[2] = (long)_col( gun );
       p1=nfmod_pol_mul(nf,prhall,p1,p2);      return rep;
       nfmod_pol_divres(nf,prhall,p1,pmod,&p1);  
     }  
     if (gcmp1(e)) break;  
   
     e=shifti(e,-1);  
     p2=nfmod_pol_sqr(nf,prhall,p2);  
     nfmod_pol_divres(nf,prhall,p2,pmod,&p2);  
   }    }
   return gerepileupto(av,p1);  
 }  
   
 static long    A = fix_relative_pol(nf,pol,0);
 isdivbyprime(GEN nf, GEN x, GEN pr)    if (degpol(nf[1]) == 1)
 {      return gerepileupto(av, factpol(simplify(pol), 0));
   GEN elt, p = (GEN)pr[1], tau = (GEN)pr[5];  
   
   elt = element_mul(nf, x, tau);    A = primpart( lift_intern(A) );
   if (divise(content(elt), p)) return 1;    if (DEBUGLEVEL>3) fprintferr("test if polynomial is square-free\n");
     g = nfgcd(A, derivpol(A), T, NULL);
   
   return 0;    A = QXQ_normalize(A, T);
 }    A = primpart(A);
   
 /* return the factor of nf.pol modulo p corresponding to pr */    if (degpol(g))
 static GEN    { /* not squarefree */
 localpol(GEN nf, GEN pr)      gpmem_t av1;
 {      GEN ex;
   long i, l;      g = QXQ_normalize(g, T);
   GEN fct, pol = (GEN)nf[1], p = (GEN)pr[1];      A = RXQX_div(A,g, T);
   
   fct = lift((GEN)factmod(pol, p)[1]);      y = nfsqff(nf,A,0); av1 = avma;
   l = lg(fct) - 1;      l = lg(y);
   for (i = 1; i <= l; i++)      ex=(GEN)gpmalloc(l * sizeof(long));
     if (isdivbyprime(nf, (GEN)fct[i], pr)) return (GEN)fct[i];      for (j=l-1; j>=1; j--)
       {
   err(talker, "cannot find a suitable factor in localpol");        GEN fact = lift((GEN)y[j]), quo = g, q;
   return NULL; /* not reached */        long e = 0;
         for(e = 1;; e++)
         {
           q = RXQX_divrem(quo,fact,T, ONLY_DIVIDES);
           if (!q) break;
           quo = q;
         }
         ex[j] = e;
       }
       avma = av1; y = gerepileupto(av, y);
       p1 = cgetg(l, t_COL); for (j=l-1; j>=1; j--) p1[j] = lstoi(ex[j]);
       free(ex);
     }
     else
     {
       y = gerepileupto(av, nfsqff(nf,A,0));
       l = lg(y);
       p1 = cgetg(l, t_COL); for (j=l-1; j>=1; j--) p1[j] = un;
     }
     if (DEBUGLEVEL>3)
       fprintferr("number of factor(s) found: %ld\n", lg(y)-1);
     rep[1] = (long)y;
     rep[2] = (long)p1; return sort_factor(rep, cmp_pol);
 }  }
   
 /* factorization of x modulo pr */  /* return a bound for T_2(P), P | polbase in C[X]
    * NB: Mignotte bound: A | S ==>
    *  |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S)
    *
    * Apply to sigma(S) for all embeddings sigma, then take the L_2 norm over
    * sigma, then take the sup over i.
    **/
 static GEN  static GEN
 nffactormod0(GEN nf, GEN x, GEN pr)  nf_Mignotte_bound(GEN nf, GEN polbase)
 {  {
   long av = avma, j, l, vx = varn(x), vn;    GEN G = gmael(nf,5,2), lS = leading_term(polbase); /* t_INT */
   GEN rep, pol, xrd, prh, p1;    GEN p1, C, N2, matGS, binlS, bin;
     long prec, i, j, d = degpol(polbase), n = degpol(nf[1]), r1 = nf_get_r1(nf);
   
   nf=checknf(nf);    if (typ(lS) == t_COL) lS = (GEN)lS[1];
   vn = varn((GEN)nf[1]);    binlS = bin = vecbinome(d-1);
   if (typ(x)!=t_POL) err(typeer,"nffactormod");    if (!gcmp1(lS)) binlS = gmul(lS, bin);
   if (vx >= vn)  
     err(talker,"polynomial variable must have highest priority in nffactormod");  
   
   prh = nfmodprinit(nf, pr);    N2 = cgetg(n+1, t_VEC);
     prec = gprecision(G);
     for (;;)
     {
       nffp_t F;
   
   if (divise((GEN)nf[4], (GEN)pr[1]))      matGS = cgetg(d+2, t_MAT);
     return gerepileupto(av, nffactormod2(nf,x,pr));      for (j=0; j<=d; j++) matGS[j+1] = lmul(G, (GEN)polbase[j+2]);
       matGS = gtrans_i(matGS);
       for (j=1; j <= r1; j++) /* N2[j] = || sigma_j(S) ||_2 */
       {
         N2[j] = lsqrt( QuickNormL2((GEN)matGS[j], DEFAULTPREC), DEFAULTPREC );
         if (lg(N2[j]) < DEFAULTPREC) goto PRECPB;
       }
       for (   ; j <= n; j+=2)
       {
         GEN q1 = QuickNormL2((GEN)matGS[j  ], DEFAULTPREC);
         GEN q2 = QuickNormL2((GEN)matGS[j+1], DEFAULTPREC);
         p1 = gmul2n(mpadd(q1, q2), -1);
         N2[j] = N2[j+1] = lsqrt( p1, DEFAULTPREC );
         if (lg(N2[j]) < DEFAULTPREC) goto PRECPB;
       }
       if (j > n) break; /* done */
   PRECPB:
       prec = (prec<<1)-2;
       remake_GM(nf, &F, prec); G = F.G;
       if (DEBUGLEVEL>1) err(warnprec, "nf_factor_bound", prec);
     }
   
   xrd = nfmod_pol_reduce(nf, prh, x);    /* Take sup over 0 <= i <= d of
   if (gcmp1((GEN)pr[4]))     * sum_sigma | binom(d-1, i-1) ||sigma(S)||_2 + binom(d-1,i) lc(S) |^2 */
   
     /* i = 0: n lc(S)^2 */
     C = mulsi(n, sqri(lS));
     /* i = d: sum_sigma ||sigma(S)||_2^2 */
     p1 = gnorml2(N2); if (gcmp(C, p1) < 0) C = p1;
     for (i = 1; i < d; i++)
   {    {
     pol = gun; /* dummy */      GEN s = gzero;
     for (j = 2; j < lg(xrd); j++)      for (j = 1; j <= n; j++)
       xrd[j] = mael(xrd, j, 1);  
     rep = factmod(xrd, (GEN)pr[1]);  
     rep = lift(rep);  
   }  
   else  
   {  
     pol = localpol(nf, pr);  
     xrd = lift(unifpol(nf, xrd, 1));  
     p1  = gun;  
     for (j = 2; j < lg(xrd); j++)  
     {      {
       xrd[j] = lmod((GEN)xrd[j], pol);        p1 = mpadd( mpmul((GEN)bin[i], (GEN)N2[j]), (GEN)binlS[i+1] );
       p1 = mpppcm(p1, denom(content((GEN)xrd[j])));        s = mpadd(s, gsqr(p1));
     }      }
     rep = factmod9(gmul(xrd, p1), (GEN)pr[1], pol);      if (gcmp(C, s) < 0) C = s;
     rep = lift(lift(rep));  
   }    }
     return C;
   }
   
   l = lg((GEN)rep[1]);  /* return a bound for T_2(P), P | polbase
   for (j = 1; j < l; j++)   * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2,
     mael(rep, 1, j) = (long)unifpol(nf, gmael(rep, 1, j), 1);   * where [P]_2 is Bombieri's 2-norm
    * Sum over conjugates
   */
   static GEN
   nf_Beauzamy_bound(GEN nf, GEN polbase)
   {
     GEN lt,C,run,s, G = gmael(nf,5,2), POL, bin;
     long i,prec,precnf, d = degpol(polbase), n = degpol(nf[1]);
   
   return gerepilecopy(av, rep);    precnf = gprecision(G);
     prec   = MEDDEFAULTPREC;
     bin = vecbinome(d);
     POL = polbase + 2;
     /* compute [POL]_2 */
     for (;;)
     {
       run= realun(prec);
       s = realzero(prec);
       for (i=0; i<=d; i++)
       {
         GEN p1 = gnorml2(gmul(G, gmul(run, (GEN)POL[i]))); /* T2(POL[i]) */
         if (!signe(p1)) continue;
         if (lg(p1) == 3) break;
         /* s += T2(POL[i]) / binomial(d,i) */
         s = addrr(s, gdiv(p1, (GEN)bin[i+1]));
       }
       if (i > d) break;
   
       prec = (prec<<1)-2;
       if (prec > precnf)
       {
         nffp_t F; remake_GM(nf, &F, prec); G = F.G;
         if (DEBUGLEVEL>1) err(warnprec, "nf_factor_bound", prec);
       }
     }
     lt = (GEN)leading_term(polbase)[1];
     s = gmul(s, mulis(sqri(lt), n));
     C = gpow(stoi(3), dbltor(1.5 + d), DEFAULTPREC); /* 3^{3/2 + d} */
     return gdiv(gmul(C, s), gmulsg(d, mppi(DEFAULTPREC)));
 }  }
   
 GEN  static GEN
 nffactormod(GEN nf, GEN x, GEN pr)  nf_factor_bound(GEN nf, GEN polbase)
 {  {
   return nffactormod0(nf, x, pr);    gpmem_t av = avma;
     GEN a = nf_Mignotte_bound(nf, polbase);
     GEN b = nf_Beauzamy_bound(nf, polbase);
     if (DEBUGLEVEL>2)
     {
       fprintferr("Mignotte bound: %Z\n",a);
       fprintferr("Beauzamy bound: %Z\n",b);
     }
     return gerepileupto(av, gmin(a, b));
 }  }
   
 extern GEN trivfact(void);  static long
   ZXY_get_prec(GEN P)
 /* factorization of x modulo pr */  
 GEN  
 nffactormod2(GEN nf,GEN pol,GEN pr)  
 {  {
   long av = avma, tetpil,lb,nbfact,psim,N,n,i,j,k,d,e,vf,r,kk;    long i, j, z, prec = 0;
   GEN y,ex,*t,f1,f2,f3,df1,g1,polb,pold,polu,vker;    for (i=2; i<lgef(P); i++)
   GEN Q,f,x,u,v,v2,v3,vz,q,vun,vzero,prhall;  
   
   nf=checknf(nf);  
   if (typ(pol)!=t_POL) err(typeer,"nffactormod");  
   if (varn(pol) >= varn(nf[1]))  
     err(talker,"polynomial variable must have highest priority in nffactormod");  
   
   prhall=nfmodprinit(nf,pr); n=degpol(nf[1]);  
   vun = gscalcol_i(gun, n);  
   vzero = gscalcol_i(gzero, n);  
   q = vker = NULL; /* gcc -Wall */  
   
   f=unifpol(nf,pol,0); f=nfmod_pol_reduce(nf,prhall,f);  
   if (!signe(f)) err(zeropoler,"nffactormod");  
   d=degpol(f); vf=varn(f);  
   if (d == 0) return trivfact();  
   t  = (GEN*)new_chunk(d+1); ex = cgetg(d+1, t_VECSMALL);  
   x=gcopy(polx[vf]); x[3]=(long)vun; x[2]=(long)vzero;  
   if (d<=1) { nbfact=2; t[1] = f; ex[1]=1; }  
   else  
   {    {
     q = (GEN)pr[1]; psim = VERYBIGINT;      GEN p = (GEN)P[i];
     if (!is_bigint(q)) psim = itos(q);      if (typ(p) == t_INT)
    /* psim has an effect only when p is small. If too big, set it to a huge  
     * number (i.e ignore it) to avoid an error in itos on next line.  
     */  
     q=gpuigs(q, itos((GEN)pr[4]));  
     f1=f; e=1; nbfact=1;  
     while (lgef(f1)>3)  
     {      {
       df1=deriv(f1,vf); f2=nfmod_pol_gcd(nf,prhall,f1,df1);        z = lgefint(p);
       g1=nfmod_pol_divres(nf,prhall,f1,f2,NULL); k=0;        if (z > prec) prec = z;
       while (lgef(g1)>3)      }
       else
       {
         for (j=2; j<lgef(p); j++)
       {        {
         k++;          z = lgefint(p[j]);
         if (k%psim == 0)          if (z > prec) prec = z;
         {  
           k++; f2=nfmod_pol_divres(nf,prhall,f2,g1,NULL);  
         }  
         f3=nfmod_pol_gcd(nf,prhall,f2,g1);  
         u = nfmod_pol_divres(nf,prhall,g1,f3,NULL);  
         f2= nfmod_pol_divres(nf,prhall,f2,f3,NULL);  
         g1=f3;  
         if (lgef(u)>3)  
         {  
           N=degpol(u); Q=cgetg(N+1,t_MAT);  
           v3=cgetg(N+1,t_COL); Q[1]=(long)v3;  
           v3[1]=(long)vun; for (i=2; i<=N; i++) v3[i]=(long)vzero;  
   
           v2 = v = nfmod_pol_pow(nf,prhall,u,x,q);  
           for (j=2; j<=N; j++)  
           {  
             v3=cgetg(N+1,t_COL); Q[j]=(long)v3;  
             for (i=1; i<=lgef(v2)-2; i++) v3[i]=v2[i+1];  
             for (; i<=N; i++) v3[i]=(long)vzero;  
             if (j<N)  
             {  
               v2=nfmod_pol_mul(nf,prhall,v2,v);  
               nfmod_pol_divres(nf,prhall,v2,u,&v2);  
             }  
           }  
           for (i=1; i<=N; i++)  
             coeff(Q,i,i)=lsub((GEN)coeff(Q,i,i),vun);  
           v2=nfkermodpr(nf,Q,prhall); r=lg(v2)-1; t[nbfact]=gcopy(u); kk=1;  
           if (r>1)  
           {  
             vker=cgetg(r+1,t_COL);  
             for (i=1; i<=r; i++)  
             {  
               v3=cgetg(N+2,t_POL);  
               v3[1]=evalsigne(1)+evallgef(2+N); setvarn(v3,vf);  
               vker[i]=(long)v3; for (j=1; j<=N; j++) v3[j+1]=coeff(v2,j,i);  
               normalizepol(v3);  
             }  
           }  
           while (kk<r)  
           {  
             v=gcopy(polun[vf]); v[2]=(long)vzero;  
             for (i=1; i<=r; i++)  
             {  
               vz=cgetg(n+1,t_COL);  
               for (j=1; j<=n; j++)  
                 vz[j] = lmodsi(mymyrand()>>8, q);  
               vz=nfreducemodpr(nf,vz,prhall);  
               v=gadd(v,nfmod_pol_mul(nf,prhall,vz,(GEN)vker[i]));  
             }  
             for (i=1; i<=kk && kk<r; i++)  
             {  
               polb=t[nbfact+i-1]; lb=lgef(polb);  
               if (lb>4)  
               {  
                 if(psim==2)  
                   polu=nfmod_split2(nf,prhall,polb,v,q);  
                 else  
                 {  
                   polu=nfmod_pol_pow(nf,prhall,polb,v,shifti(q,-1));  
                   polu=gsub(polu,vun);  
                 }  
                 pold=nfmod_pol_gcd(nf,prhall,polb,polu);  
                 if (lgef(pold)>3 && lgef(pold)<lb)  
                 {  
                   t[nbfact+i-1]=pold; kk++;  
                   t[nbfact+kk-1]=nfmod_pol_divres(nf,prhall,polb,pold,NULL);  
                 }  
               }  
             }  
           }  
           for (i=nbfact; i<nbfact+r; i++) ex[i]=e*k;  
           nbfact+=r;  
         }  
       }        }
       e*=psim; j=(degpol(f2))/psim+3; f1=cgetg(j,t_POL);  
       f1[1] = evalsigne(1) | evallgef(j) | evalvarn(vf);  
       for (i=2; i<j; i++)  
         f1[i]=(long)element_powmodpr(nf,(GEN)f2[psim*(i-2)+2],  
                                      gdiv(q,(GEN)pr[1]),prhall);  
     }      }
   }    }
   if (nbfact < 2)    return prec + 1;
     err(talker, "%Z not a prime (nffactormod)", q);  }
   for (j=1; j<nbfact; j++)  
   long
   ZM_get_prec(GEN x)
   {
     long i, j, l, k = 2, lx = lg(x);
   
     for (j=1; j<lx; j++)
   {    {
     v = element_divmodpr(nf,vun,gmael(t,j,lgef(t[j])-1),prhall);      GEN c = (GEN)x[j];
     t[j] = unifpol(nf,nfmod_pol_mul(nf,prhall,v,(GEN)t[j]),1);      for (i=1; i<lx; i++) { l = lgefint(c[i]); if (l > k) k = l; }
   }    }
     return k;
   tetpil=avma; y=cgetg(3,t_MAT);  
   u=cgetg(nbfact,t_COL); y[1]=(long)u;  
   v=cgetg(nbfact,t_COL); y[2]=(long)v;  
   for (j=1,k=0; j<nbfact; j++)  
     if (ex[j])  
     {  
       k++;  
       u[k]=lcopy((GEN)t[j]);  
       v[k]=lstoi(ex[j]);  
     }  
   return gerepile(av,tetpil,y);  
 }  }
   long
 /* return pol + pol^2 + ... + pol^(q/2) modulo prhall and  ZX_get_prec(GEN x)
    the polynomial pmod */  
 static GEN  
 nfmod_split2(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN exp)  
 {  {
   long av = avma;    long j, l, k = 2, lx = lgef(x);
   GEN p1,p2,q;  
   
   if (cmpis(exp,2)<=0) return pol;    for (j=2; j<lx; j++)
   p2=p1=pol; q=shifti(exp,-1);  
   while (!gcmp1(q))  
   {    {
     p2=nfmod_pol_sqr(nf,prhall,p2);      l = lgefint(x[j]); if (l > k) k = l;
     nfmod_pol_divres(nf,prhall,p2,pmod,&p2);  
     q=shifti(q,-1); p1=gadd(p1,p2);  
   }    }
   return gerepileupto(av,p1);    return k;
 }  }
   
 /* If p doesn't divide either a or b and has a divisor of degree 1, return it.  
  * Return NULL otherwise.  
  */  
 static GEN  static GEN
 p_ok(GEN nf, GEN p, GEN a)  complex_bound(GEN P)
 {  {
   long av,m,i;    return gmul(max_modulus(P, 0.01), dbltor(1.0101)); /* exp(0.01) = 1.0100 */
   GEN dec;  
   
   if (divise(a,p)) return NULL;  
   av = avma; dec = primedec(nf,p); m=lg(dec);  
   for (i=1; i<m; i++)  
   {  
     GEN pr = (GEN)dec[i];  
     if (is_pm1(pr[4]))  
       return pr;  
   }  
   avma = av; return NULL;  
 }  }
   
 /* for each new prime ct--, if ct = 0, return NULL */  /* return Bs: if r a root of sigma_i(P), |r| < Bs[i] */
 static GEN  static GEN
 choose_prime(GEN nf, GEN dk, GEN lim, long ct)  nf_root_bounds(GEN P, GEN T)
 {  {
   GEN p, pr;    long lR, i, j, l, prec;
     GEN Ps, R, V, nf;
   
   p = nextprime(lim);    if (isrational(P)) return complex_bound(P);
   for (;;)  
     T = get_nfpol(T, &nf);
   
     prec = ZXY_get_prec(P);
     l = lgef(P);
     if (nf && nfgetprec(nf) >= prec)
       R = (GEN)nf[6];
     else
       R = roots(T, prec);
     lR = lg(R);
     V = cgetg(lR, t_VEC);
     Ps = cgetg(l, t_POL); /* sigma (P) */
     Ps[1] = P[1];
     for (j=1; j<lg(R); j++)
   {    {
     if ((pr = p_ok(nf,p,dk))) break;      GEN r = (GEN)R[j];
     ct--;      for (i=2; i<l; i++) Ps[i] = (long)poleval((GEN)P[i], r);
     if (!ct) return NULL;      V[j] = (long)complex_bound(Ps);
     p = nextprime(addis(p,2));  
   }    }
     return V;
   return pr;  
 }  }
   
 /* test if the discriminant of polbase modulo some few primes  /* return B such that if x in O_K, K = Z[X]/(T), then the L2-norm of the
    is non-zero. Return 1 if it is so (=> polbase is square-free)   * coordinates of the numerator of x [on the power, resp. integral, basis if T
    and 0 otherwise (=> polbase may or may not be square-free) */   * is a polynomial, resp. an nf] is  <= B T_2(x)
 static int   * *ptden = multiplicative bound for denom(x) */
 is_sqf(GEN nf, GEN polbase)  static GEN
   L2_bound(GEN T, GEN tozk, GEN *ptden)
 {  {
   GEN lt, pr, prh, p2, p;    GEN nf, M, L, prep, den, u;
   long i, d = lgef(polbase), ct = 5;    long prec;
   
   lt = (GEN)leading_term(polbase)[1];    T = get_nfpol(T, &nf);
   p  = stoi(101);    u = NULL; /* gcc -Wall */
   
   while (ct > 0)    prec = ZX_get_prec(T) + ZM_get_prec(tozk);
   {    den = nf? gun: NULL;
     /* small primes tend to divide discriminants more often    den = initgaloisborne(T, den, prec, &L, &prep, NULL);
        than large ones so we look at primes >= 101 */    M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
     pr = choose_prime(nf,lt,p,30);    if (nf) M = gmul(tozk, M);
     if (!pr) break;    if (gcmp1(den)) den = NULL;
     *ptden = den;
     return QuickNormL2(M, DEFAULTPREC);
   }
   
     p=(GEN)pr[1];  /* || L ||_p^p in dimension n (L may be a scalar) */
     prh=prime_to_ideal(nf,pr);  static GEN
   normlp(GEN L, long p, long n)
   {
     long i,l, t = typ(L);
     GEN z;
   
     p2=gcopy(polbase);    if (!is_vec_t(t)) return gmulsg(n, gpowgs(L, p));
     lt=mpinvmod(lt,p);  
   
     for (i=2; i<d; i++)    l = lg(L); z = gzero;
       p2[i] = nfreducemodpr_i(gmul(lt,(GEN)p2[i]), prh)[1];    /* assert(n == l-1); */
     p2 = normalizepol(p2);    for (i=1; i<l; i++)
       z = gadd(z, gpowgs((GEN)L[i], p));
     /* discriminant is non-zero => polynomial is square-free */    return z;
     if (!gcmp0(p2) && !divise(discsr(p2),p))  { return 1; }  
   
     ct--;  
     p=addis(p,1);  
   }  
   
   return 0;  
 }  }
   
 /* rescale p in K[X] (coeffs in algtobasis form) --> primitive in O_K[X] */  /* S = S0 + qS1, P = P0 + qP1 (Euclidean div. by q integer). For a true
 GEN   * factor (vS, vP), we have:
 nf_pol_to_int(GEN p, GEN *den)   *    | S vS + P vP |^2 < Btra
    * This implies | S1 vS + P1 vP |^2 < Blow, assuming q > sqrt(Btra).
    * d = dimension of low part (= [nf:Q])
    * n0 = bound for |vS|^2
    * */
   static double
   get_Blow(long n0, long d, GEN q)
 {  {
   long i, l = lgef(p);  #if 0
   GEN d = gun;    double sqrtn0 = sqrt((double)n0);
   for (i=2; i<l; i++)    double t = sqrtn0 + sqrt((double)d) * (sqrtn0 + 1);
     d = mpppcm(d,denom((GEN)p[i]));  #else
   if (! gcmp1(d)) p = gmul(p, d);    double sqrtd = sqrt((double)d);
   *den = d; return p;    double t, aux;
   
     if (gexpo(q)>30)
       aux = 1.0001;
     else
     {
       aux = gtodouble(q);
       aux = sqrt(1 + 4/(aux*aux));
     }
     t = n0*sqrtd + sqrtd/2 * aux * (sqrtd * (n0+1)); /* assume pr degree 1 */
   #endif
     t = 1. + 0.5 * t;
     return t * t;
 }  }
   
 /* return the roots of pol in nf */  typedef struct {
 GEN    GEN d;
 nfroots(GEN nf,GEN pol)    GEN dPinvS;   /* d P^(-1) S   [ integral ] */
     double **PinvSdbl; /* P^(-1) S as double */
     GEN S1, P1;   /* S = S0 + S1 q, idem P */
   } trace_data;
   
   /* S1 * u - P1 * round(P^-1 S u). K non-zero coords in u given by ind */
   static GEN
   get_trace(GEN ind, trace_data *T)
 {  {
   long av=avma,tetpil,d=lgef(pol),fl;    long i, j, l, K = lg(ind)-1;
   GEN p1,p2,polbase,polmod,den;    GEN z, s, v;
   
   p2=NULL; /* gcc -Wall */    s = (GEN)T->S1[ ind[1] ];
   nf=checknf(nf);    if (K == 1) return s;
   if (typ(pol)!=t_POL) err(talker,"not a polynomial in nfroots");  
   if (varn(pol) >= varn(nf[1]))  
     err(talker,"polynomial variable must have highest priority in nfroots");  
   
   polbase=unifpol(nf,pol,0);    /* compute s = S1 u */
     for (j=2; j<=K; j++) s = gadd(s, (GEN)T->S1[ ind[j] ]);
   
   if (d==3)    /* compute v := - round(P^1 S u) */
     l = lg(s);
     v = cgetg(l, t_VECSMALL);
     for (i=1; i<l; i++)
   {    {
     tetpil=avma; p1=cgetg(1,t_VEC);      double r, t = 0.;
     return gerepile(av,tetpil,p1);      /* quick approximate computation */
       for (j=1; j<=K; j++) t += T->PinvSdbl[ ind[j] ][i];
       r = floor(t + 0.5);
       if (fabs(t + 0.5 - r) < 0.01)
       { /* dubious, compute exactly */
         z = gzero;
         for (j=1; j<=K; j++) z = addii(z, ((GEN**)T->dPinvS)[ ind[j] ][i]);
         v[i] = - itos( diviiround(z, T->d) );
       }
       else
         v[i] = - (long)r;
   }    }
     return gadd(s, gmul_mati_smallvec(T->P1, v));
   }
   
   if (d==4)  static trace_data *
   {  init_trace(trace_data *T, GEN S, nflift_t *L, GEN q)
     tetpil=avma; p1=cgetg(2,t_VEC);  {
     p1[1] = (long)basistoalg(nf,gneg_i(    long e = gexpo((GEN)S), i,j, l,h;
       element_div(nf,(GEN)polbase[2],(GEN)polbase[3])));    GEN qgood, S1, invd;
     return gerepile(av,tetpil,p1);  
   }  
   
   p1=element_inv(nf,leading_term(polbase));    if (e < 0) return NULL; /* S = 0 */
   polbase=nf_pol_mul(nf,p1,polbase);  
   
   polbase = nf_pol_to_int(polbase, &den);    qgood = shifti(gun, e - 32); /* single precision check */
   polmod=unifpol(nf,polbase,1);    if (cmpii(qgood, q) > 0) q = qgood;
   
   if (DEBUGLEVEL>=4)    S1 = gdivround(S, q);
     fprintferr("test if the polynomial is square-free\n");    if (gcmp0(S1)) return NULL;
   
   fl = is_sqf(nf, polbase);    invd = ginv(itor(L->den, DEFAULTPREC));
   
   /* the polynomial may not be square-free ... */    T->dPinvS = gmul(L->iprk, S);
   if (!fl)    l = lg(S);
     h = lg(T->dPinvS[1]);
     T->PinvSdbl = (double**)cgetg(l, t_MAT);
     init_dalloc();
     for (j = 1; j < l; j++)
   {    {
     p1=derivpol(polmod);      double *t = dalloc(h * sizeof(double));
     p2=nf_pol_subres(nf,polmod,p1);      GEN c = (GEN)T->dPinvS[j];
     if (degpol(p2) == 0) fl = 1;      T->PinvSdbl[j] = t;
       for (i=1; i < h; i++) t[i] = rtodbl(gmul(invd, (GEN)c[i]));
   }    }
   
   if (!fl)    T->d  = L->den;
   {    T->P1 = gdivround(L->prk, q);
     p1=element_inv(nf,leading_term(p2));    T->S1 = S1; return T;
     p2=nf_pol_mul(nf,p1,p2);  
     polmod=nf_pol_divres(nf,polmod,p2,NULL);  
   
     p1=element_inv(nf,leading_term(polmod));  
     polmod=nf_pol_mul(nf,p1,polmod);  
   
     polmod = nf_pol_to_int(polmod, &den);  
     polmod=unifpol(nf,polmod,1);  
   }  
   
   p1 = nfsqff(nf,polmod,1);  
   tetpil=avma; return gerepile(av, tetpil, gen_sort(p1, 0, cmp_pol));  
 }  }
   
 /* return a minimal lift of elt modulo id */  static void
 static GEN  update_trace(trace_data *T, long k, long i)
 nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt)  
 {  {
   GEN u = gmul(idinv,elt);    if (!T) return;
   long i, l = lg(u);    T->S1[k]       = T->S1[i];
   for (i=1; i<l; i++) u[i] = (long)gdivround((GEN)u[i], den);    T->dPinvS[k]   = T->dPinvS[i];
   return gsub(elt, gmul(id, u));    T->PinvSdbl[k] = T->PinvSdbl[i];
 }  }
   
 /* return the lift of pol with coefficients of T2-norm <= C (if possible) */  /* Naive recombination of modular factors: combine up to maxK modular
    * factors, degree <= klim and divisible by hint
    *
    * target = polynomial we want to factor
    * famod = array of modular factors.  Product should be congruent to
    * target/lc(target) modulo p^a
    * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
 static GEN  static GEN
 nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol)  nfcmbf(nfcmbf_t *T, GEN p, long a, long maxK, long klim)
 {  {
   long i, d = lgef(pol);    long Sbound;
   GEN x = cgetg(d,t_POL);    GEN pol = T->pol, nf = T->nf, famod = T->fact;
   x[1] = pol[1];    GEN bound = T->bound;
   for (i=2; i<d; i++)    GEN nfpol = (GEN)nf[1];
     x[i] = (long) nf_bestlift(id,idinv,den,(GEN)pol[i]);    long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol);
   return x;    GEN lc, lcpol;
 }    GEN pk = gpowgs(p,a), pas2 = shifti(pk,-1);
   
 #if 0    GEN trace1   = cgetg(lfamod+1, t_MAT);
 /* return pol(elt) */    GEN trace2   = cgetg(lfamod+1, t_MAT);
 static GEN    GEN ind      = cgetg(lfamod+1, t_VECSMALL);
 nf_pol_eval(GEN nf,GEN pol,GEN elt)    GEN degpol   = cgetg(lfamod+1, t_VECSMALL);
 {    GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
   long av=avma,tetpil,i;    GEN listmod  = cgetg(lfamod+1, t_COL);
   GEN p1;    GEN fa       = cgetg(lfamod+1, t_COL);
     GEN res = cgetg(3, t_VEC);
     GEN q = ceil_safe(mpsqrt(T->BS_2));
     const double Blow = get_Blow(lfamod, dnf, q);
     trace_data _T1, _T2, *T1, *T2;
   
   i=lgef(pol)-1; if (i==2) return gcopy((GEN)pol[2]);    if (maxK < 0) maxK = lfamod-1;
   
   p1=element_mul(nf,(GEN)pol[i],elt);    lc = absi(leading_term(pol));
   for (i--; i>=3; i--)    if (gcmp1(lc)) lc = NULL;
     p1=element_mul(nf,elt,gadd((GEN)pol[i],p1));    lcpol = lc? gmul(lc,pol): pol;
   tetpil=avma; return gerepile(av,tetpil,gadd(p1,(GEN)pol[2]));  
 }  
 #endif  
   
 /* return the factorization of x in nf */  
 GEN  
 nffactor(GEN nf,GEN pol)  
 {  
   GEN y,p1,p2,den,p3,p4,quot,rep=cgetg(3,t_MAT);  
   long av = avma,tetpil,i,j,d,fl;  
   
   if (DEBUGLEVEL >= 4) timer2();  
   
   p3=NULL; /* gcc -Wall */  
   nf=checknf(nf);  
   if (typ(pol)!=t_POL) err(typeer,"nffactor");  
   if (varn(pol) >= varn(nf[1]))  
     err(talker,"polynomial variable must have highest priority in nffactor");  
   
   d=lgef(pol);  
   if (d==3)  
   {    {
     rep[1]=lgetg(1,t_COL);      GEN t1,t2, lc2 = lc? sqri(lc): NULL;
     rep[2]=lgetg(1,t_COL);  
     return rep;  
   }  
   if (d==4)  
   {  
     p1=cgetg(2,t_COL); rep[1]=(long)p1; p1[1]=lcopy(pol);  
     p1=cgetg(2,t_COL); rep[2]=(long)p1; p1[1]=un;  
     return rep;  
   }  
   
   p1=element_inv(nf,leading_term(pol));      for (i=1; i <= lfamod; i++)
   pol=nf_pol_mul(nf,p1,pol);      {
         GEN P = (GEN)famod[i];
         long d = degpol(P);
   
   p1=unifpol(nf,pol,0);        degpol[i] = d; P += 2;
   p1 = nf_pol_to_int(p1, &den);        t1 = (GEN)P[d-1];/* = - S_1 */
         t2 = sqri(t1);
         if (d > 1) t2 = subii(t2, shifti((GEN)P[d-2],1));
         t2 = modii(t2, pk); /* = S_2 Newton sum */
         if (lc)
         {
           t1 = modii(mulii(lc, t1), pk);
           t2 = modii(mulii(lc2,t2), pk);
         }
         trace1[i] = (long)nf_bestlift(t1, NULL, T->L);
         trace2[i] = (long)nf_bestlift(t2, NULL, T->L);
       }
   
   if (DEBUGLEVEL>=4)      T1 = init_trace(&_T1, trace1, T->L, q);
     fprintferr("test if the polynomial is square-free\n");      T2 = init_trace(&_T2, trace2, T->L, q);
   
   fl = is_sqf(nf, p1);  
   
   /* polynomial may not be square-free ... */  
   if (!fl)  
   {  
     p2=derivpol(p1);  
     p3=nf_pol_subres(nf,p1,p2);  
     if (degpol(p3) == 0) fl = 1;  
   }    }
     degsofar[0] = 0; /* sentinel */
   
   if (!fl)    /* ind runs through strictly increasing sequences of length K,
   {     * 1 <= ind[i] <= lfamod */
     p4=element_inv(nf,leading_term(p3));  nextK:
     p3=nf_pol_mul(nf,p4,p3);    if (K > maxK || 2*K > lfamod) goto END;
     if (DEBUGLEVEL > 3)
       fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K));
     setlg(ind, K+1); ind[1] = 1;
     Sbound = ((K+1)>>1);
     i = 1; curdeg = degpol[ind[1]];
     for(;;)
     { /* try all combinations of K factors */
       for (j = i; j < K; j++)
       {
         degsofar[j] = curdeg;
         ind[j+1] = ind[j]+1; curdeg += degpol[ind[j+1]];
       }
       if (curdeg <= klim && curdeg % T->hint == 0) /* trial divide */
       {
         GEN t, y, q, list;
         gpmem_t av;
   
     p2=nf_pol_divres(nf,p1,p3,NULL);        av = avma;
     p4=element_inv(nf,leading_term(p2));        /* d - 1 test */
     p2=nf_pol_mul(nf,p4,p2);        if (T1)
         {
           t = get_trace(ind, T1);
           if (rtodbl(QuickNormL2(t,DEFAULTPREC)) > Blow)
           {
             if (DEBUGLEVEL>6) fprintferr(".");
             avma = av; goto NEXT;
           }
         }
         /* d - 2 test */
         if (T2)
         {
           t = get_trace(ind, T2);
           if (rtodbl(QuickNormL2(t,DEFAULTPREC)) > Blow)
           {
             if (DEBUGLEVEL>3) fprintferr("|");
             avma = av; goto NEXT;
           }
         }
         avma = av;
         y = lc; /* full computation */
         for (i=1; i<=K; i++)
         {
           GEN q = (GEN)famod[ind[i]];
           if (y) q = gmul(y, q);
           y = centermod_i(q, pk, pas2);
         }
         y = nf_pol_lift(y, bound, T);
         if (!y)
         {
           if (DEBUGLEVEL>3) fprintferr("@");
           avma = av; goto NEXT;
         }
         /* try out the new combination: y is the candidate factor */
         q = RXQX_divrem(lcpol,y, nfpol, ONLY_DIVIDES);
         if (!q)
         {
           if (DEBUGLEVEL>3) fprintferr("*");
           avma = av; goto NEXT;
         }
   
     p2 = nf_pol_to_int(p2, &den);        /* found a factor */
         list = cgetg(K+1, t_VEC);
         listmod[cnt] = (long)list;
         for (i=1; i<=K; i++) list[i] = famod[ind[i]];
   
     p2=unifpol(nf,p2,1);        y = primpart(y);
     tetpil = avma; y = nfsqff(nf,p2,0);        fa[cnt++] = (long)QXQ_normalize(y, nfpol);
     i = nfcmbf.nfact;        /* fix up pol */
         pol = q;
         if (lc) pol = primpart(pol);
         for (i=j=k=1; i <= lfamod; i++)
         { /* remove used factors */
           if (j <= K && i == ind[j]) j++;
           else
           {
             famod[k] = famod[i];
             update_trace(T1, k, i);
             update_trace(T2, k, i);
             degpol[k] = degpol[i]; k++;
           }
         }
         lfamod -= K;
         if (lfamod < 2*K) goto END;
         i = 1; curdeg = degpol[ind[1]];
         if (lc) lc = absi(leading_term(pol));
         lcpol = lc? gmul(lc,pol): pol;
         if (DEBUGLEVEL > 2)
         {
           fprintferr("\n"); msgtimer("to find factor %Z",y);
           fprintferr("remaining modular factor(s): %ld\n", lfamod);
         }
         continue;
       }
   
     quot=nf_pol_divres(nf,p1,p2,NULL);  NEXT:
     p3=(GEN)gpmalloc((i+1) * sizeof(long));      for (i = K+1;;)
     for (j=i; j>=1; j--)  
     {      {
       GEN fact=(GEN)y[j], quo = quot, rem;        if (--i == 0) { K++; goto nextK; }
       long e = 0;        if (++ind[i] <= lfamod - K + i)
   
       do  
       {        {
         quo = nf_pol_divres(nf,quo,fact,&rem);          curdeg = degsofar[i-1] + degpol[ind[i]];
         e++;          if (curdeg <= klim) break;
       }        }
       while (gcmp0(rem));  
       p3[j]=lstoi(e);  
     }      }
     avma = (long)y; y = gerepile(av, tetpil, y);  
     p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=lcopy((GEN)p3[i]);  
     free(p3);  
   }    }
   else  END:
   {    if (degpol(pol) > 0)
     tetpil=avma; y = gerepile(av,tetpil,nfsqff(nf,p1,0));    { /* leftover factor */
     i = nfcmbf.nfact;      if (signe(leading_term(pol)) < 0) pol = gneg_i(pol);
     p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=un;  
   }  
   if (DEBUGLEVEL>=4)  
     fprintferr("number of factor(s) found: %ld\n", nfcmbf.nfact);  
   rep[1]=(long)y;  
   rep[2]=(long)p2; return sort_factor(rep, cmp_pol);  
 }  
   
 /* test if the matrix M is suitable */      if (lc && lfamod < 2*K) pol = QXQ_normalize(primpart(pol), nfpol);
 static long      setlg(famod, lfamod+1);
 test_mat(GEN M, GEN p, GEN C2, long k)      listmod[cnt] = (long)dummycopy(famod);
 {      fa[cnt++] = (long)pol;
   long av = avma, i, N = lg(M);  
   GEN min, prod, L2, R;  
   
   min = prod = gcoeff(M,1,1);  
   for (i = 2; i < N; i++)  
   {  
     L2 = gcoeff(M,i,i); prod = mpmul(prod,L2);  
     if (mpcmp(L2,min) < 0) min = L2;  
   }    }
   R = mpmul(min, gpowgs(p, k<<1));    if (DEBUGLEVEL>6) fprintferr("\n");
   i = mpcmp(mpmul(C2,prod), R); avma = av;    setlg(listmod, cnt); setlg(fa, cnt);
   return (i < 0);    res[1] = (long)fa;
     res[2] = (long)listmod; return res;
 }  }
   
 /* return the matrix corresponding to pr^e with R(pr^e) > C */  
 static GEN  static GEN
 T2_matrix_pow(GEN nf, GEN pre, GEN p, GEN C, long *kmax, long prec)  nf_check_factors(nfcmbf_t *T, GEN P, GEN M_L, GEN famod, GEN pk)
 {  {
   long av=avma,av1,lim, k = *kmax, N = degpol((GEN)nf[1]);    GEN nf = T->nf, bound = T->bound;
   int tot_real = !signe(gmael(nf,2,2));    GEN nfT = (GEN)nf[1];
   GEN p1,p2,p3,u,C2,T2, x = (GEN)nf[1];    long i, j, r, n0;
     GEN pol = P, lcpol, lc, list, piv, y, pas2;
   
   C2 = gdiv(gmul2n(C,2), absi((GEN)nf[3]));    piv = special_pivot(M_L);
   p1 = gmul(pre,lllintpartial(pre)); av1 = avma;    if (!piv) return NULL;
   T2 = tot_real? gmael(nf,5,4)    if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv);
                : nf_get_T2((GEN) nf[7], roots(x,prec));  
   p3 = qf_base_change(T2,p1,1);  
   
   if (N <= 6 && test_mat(p3,p,C2,k))    pas2 = shifti(pk,-1);
     r  = lg(piv)-1;
     n0 = lg(piv[1])-1;
     list = cgetg(r+1, t_COL);
     lc = absi(leading_term(pol));
     if (is_pm1(lc)) lc = NULL;
     lcpol = lc? gmul(lc, pol): pol;
     for (i=1; i<r; i++)
   {    {
     avma = av1; return gerepileupto(av,p1);      GEN c = (GEN)piv[i];
   }      if (DEBUGLEVEL) fprintferr("nf_LLL_cmbf: checking factor %ld\n",i);
   
   av1=avma; lim = stack_lim(av1,2);      y = lc;
   for (;;)      for (j=1; j<=n0; j++)
   {        if (signe(c[j]))
     if (DEBUGLEVEL>2) fprintferr("exponent: %ld\n",k);        {
           GEN q = (GEN)famod[j];
           if (y) q = gmul(y, q);
           y = centermod_i(q, pk, pas2);
         }
       y = nf_pol_lift(y, bound, T);
       if (!y) return NULL;
   
     for(;;)      /* y is the candidate factor */
     {      pol = RXQX_divrem(lcpol,y,nfT, ONLY_DIVIDES);
       u = tot_real? lllgramall(p3,2,lll_IM) : lllgramintern(p3,2,1,prec);      if (!pol) return NULL;
       if (u) break;  
   
       prec=(prec<<1)-2;      y = primpart(y);
       if (DEBUGLEVEL > 1) err(warnprec,"nffactor[1]",prec);      if (lc)
       T2 = nf_get_T2((GEN) nf[7], roots(x,prec));  
       p3 = qf_base_change(T2,p1,1);  
     }  
     if (DEBUGLEVEL>2) msgtimer("lllgram + base change");  
     p3 = qf_base_change(p3,u,1);  
     if (test_mat(p3,p,C2,k))  
     {      {
       *kmax = k;        pol = primpart(pol);
       return gerepileupto(av,gmul(p1,u));        lc = absi(leading_term(pol));
     }      }
       lcpol = lc? gmul(lc, pol): pol;
     /* we also need to increase the precision */      list[i] = (long)QXQ_normalize(y, nfT);
     p2=shifti(gceil(mulsr(k, glog(p,DEFAULTPREC))),-1);  
     prec += (long)(itos(p2)*pariK1);  
     if (DEBUGLEVEL > 1) err(warnprec,"nffactor[2]",prec);  
     k = k<<1; p1 = idealmullll(nf,p1,p1);  
     if (low_stack(lim, stack_lim(av1,2)))  
     {  
       if (DEBUGMEM>1) err(warnmem,"T2_matrix_pow");  
       p1 = gerepileupto(av1,p1);  
     }  
     if (!tot_real) T2 = nf_get_T2((GEN) nf[7], roots(x,prec));  
     p3 = qf_base_change(T2,p1,1);  
   }    }
     y = primpart(pol);
     list[i] = (long)QXQ_normalize(y, nfT); return list;
 }  }
   
 /* return the factorization of the square-free polynomial x.  
    The coeff of x are in Z_nf and its leading term is a rational  
    integer. If fl = 1,return only the roots of x in nf */  
 static GEN  static GEN
 nfsqff(GEN nf,GEN pol, long fl)  nf_to_Zp(GEN x, GEN pk, GEN pas2, GEN ffproj)
 {  {
   long d=lgef(pol),i,k,m,n,av=avma,tetpil,newprec,prec,nbf=BIGINT,anbf,ct=5;    return centermodii(gmul(ffproj, x), pk, pas2);
   GEN p1,pr,p2,rep,k2,C,h,dk,dki,p,prh,p3,T2,polbase,fact,pk,ap,apr;  }
   GEN polmod,polred,hinv,lt,minp,den,maxp=shifti(gun,32),run,aprh;  
   
   if (DEBUGLEVEL>=4) msgtimer("square-free");  /* assume lc(pol) != 0 mod prh. Reduce P to Zp[X] mod p^a */
   static GEN
   ZpX(GEN P, GEN pk, GEN ffproj)
   {
     long i, l;
     GEN z, pas2 = shifti(pk,-1);
   
   dk=absi((GEN)nf[3]);    if (typ(P)!=t_POL) return nf_to_Zp(P,pk,pas2,ffproj);
   dki=mulii(dk,(GEN)nf[4]);    l = lgef(P);
   n=degpol(nf[1]);    z = cgetg(l,t_POL); z[1] = P[1];
     for (i=2; i<l; i++) z[i] = (long)nf_to_Zp((GEN)P[i],pk,pas2,ffproj);
     return normalizepol(z);
   }
   
   polbase = unifpol(nf,pol,0);  /* We want to be able to reconstruct x, |x|^2 < C, from x mod pr^k
   polmod  = unifpol(nf,pol,1);   * We want B := min B_i >= C. Let alpha the LLL parameter,
   dki=mulii(dki,gnorm((GEN)polmod[d-1]));   * y = 1/(alpha-1/4) > 4/3, the theoretical guaranteed bound follows from
    *    (Npr)^(2k/n) = det(L)^(2/n) <= 1/n sum B_i
    *                                <= 1/n B sum y^i
    *                                <= 1/n B y^(n-1) / (y-1)
    *
    *  Hence log(B/n) >= 2k/n log(Npr) - (n-1)log(y) + log(y-1)
    *                 >= log(C/n), provided
    *   k >= [ log(C/n) + (n-1)log(y) - log(y-1) ] * n/(2log(Npr))
    */
   static double
   bestlift_bound(GEN C, long n, double alpha, GEN Npr)
   {
     const double y = 1 / (alpha - 0.25); /* = 2 if alpha = 3/4 */
     double t;
     if (typ(C) != t_REAL) C = gmul(C, realun(DEFAULTPREC));
     setlg(C, DEFAULTPREC);
     t = rtodbl(mplog(divrs(C,n))) + (n-1) * log(y) - log(y-1);
     return ceil((t * n) / (2.* log(gtodouble(Npr))));
   }
   
   prec = DEFAULTPREC;  static void
   for (;;)  bestlift_init(long a, GEN nf, GEN pr, GEN C, nflift_t *T)
   {  {
     if (prec <= gprecision(nf))    const int D = 4;
       T2 = gprec_w(gmael(nf,5,3), prec);    const double alpha = ((double)D-1) / D; /* LLL parameter */
     else    gpmem_t av = avma;
       T2 = nf_get_T2((GEN)nf[7], roots((GEN)nf[1], prec));    GEN prk, PRK, B, GSmin, pk;
   
     run=realun(prec);    if (!a) a = (long)bestlift_bound(C, degpol(nf[1]), alpha, idealnorm(nf,pr));
     p1=realzero(prec);  
     for (i=2; i<d; i++)  
     {  
       p2 = gmul(run, (GEN)polbase[i]);  
       p2 = qfeval(T2, p2);  
       p1 = addrr(p1, gdiv(p2, binome(stoi(d), i-2)));  
       if (signe(p1) < 0) break;  
     }  
   
     if (signe(p1) > 0) break;    for (;; avma = av, a<<=1)
     {
       if (DEBUGLEVEL>2) fprintferr("exponent: %ld\n",a);
       prk = idealpows(nf, pr, a);
       pk = gcoeff(prk,1,1);
       PRK = hnfmodid(prk, pk);
   
     prec = (prec<<1)-2;      PRK = lllint_i(PRK, D, 0, NULL, NULL, &B);
     if (DEBUGLEVEL > 1) err(warnprec, "nfsqff", prec);      GSmin = vecmin(GS_norms(B, DEFAULTPREC));
       if (gcmp(GSmin, C) >= 0) break;
   }    }
     if (DEBUGLEVEL>2) fprintferr("for this exponent, GSmin = %Z\n",GSmin);
     T->k = a;
     T->pk = T->den = pk;
     T->prk = PRK;
     T->iprk = ZM_inv(PRK, pk);
     T->GSmin= GSmin;
     T->ZpProj = dim1proj(prk); /* nf --> Zp */
   }
   
   lt = (GEN)leading_term(polbase)[1];  /* assume pr has degree 1 */
   p1 = gmul(p1, mulis(sqri(lt), n));  static GEN
   C = gpow(stoi(3), gadd(gmulsg(3, ghalf), stoi(d)), prec);  nf_LLL_cmbf(nfcmbf_t *T, GEN p, long k, long rec)
   C = gdiv(gmul(C, p1), gmulsg(d, mppi(prec)));  {
     nflift_t *L = T->L;
   if (DEBUGLEVEL>=4)    GEN pk = L->pk, PRK = L->prk, PRKinv = L->iprk, GSmin = L->GSmin;
     fprintferr("the bound on the T2-norm of the coeff. is: %Z\n", C);  
   
   /* the theoretical bound for the exponent should be:    GEN famod = T->fact, nf = T->nf, ZC = T->ZC, Br = T->Br;
      k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.347))); */    GEN Pbase = T->polbase, P = T->pol, dn = T->dn;
   k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.15)));    GEN nfT = (GEN)nf[1];
   k2=gmul2n(gmulgs(k2,n),-1);    GEN Btra;
     long dnf = degpol(nfT), dP = degpol(P);
   
   minp=gmin(gexp(gmul2n(k2,-6),BIGDEFAULTPREC), maxp);    double BitPerFactor = 0.5; /* nb bits / modular factor */
   minp=gceil(minp);    long i, C, tmax, n0;
     GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
   if (DEBUGLEVEL>=4)    double Blow;
   {    gpmem_t av, av2, lim;
     fprintferr("lower bound for the prime numbers: %Z\n", minp);    long ti_LLL = 0, ti_CF = 0;
     msgtimer("bounds computation");    pari_timer ti, ti2, TI;
   }  
   
   p = rep = polred = NULL; /* gcc -Wall */    if (DEBUGLEVEL>2) (void)TIMER(&ti);
   pr=NULL;  
   for (;;)    lP = absi(leading_term(P));
     if (is_pm1(lP)) lP = NULL;
   
     n0 = lg(famod) - 1;
    /* Lattice: (S PRK), small vector (vS vP). To find k bound for the image,
     * write S = S1 q + S0, P = P1 q + P0
     * |S1 vS + P1 vP|^2 <= Blow for all (vS,vP) assoc. to true factors */
     Btra = mulrr(ZC, mulsr(dP*dP, normlp(Br, 2, dnf)));
     Blow = get_Blow(n0, dnf, gceil(Btra));
     C = (long)ceil(sqrt(Blow/n0)) + 1; /* C^2 n0 ~ Blow */
     Bnorm = dbltor( n0 * C * C + Blow );
     ZERO = zeromat(n0, dnf);
   
     av = avma; lim = stack_lim(av, 1);
     TT = cgetg(n0+1, t_VEC);
     Tra  = cgetg(n0+1, t_MAT);
     for (i=1; i<=n0; i++) TT[i] = 0;
     CM_L = gscalsmat(C, n0);
     /* tmax = current number of traces used (and computed so far) */
     for(tmax = 0;; tmax++)
   {    {
     apr = choose_prime(nf,dki,minp, pr?30:0);      long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1;
     if (!apr) break;      GEN oldCM_L, M_L, q, S1, P1, VV;
       int first = 1;
   
     ap=(GEN)apr[1];      /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */
     aprh=prime_to_ideal(nf,apr);      Btra = mulrr(ZC, mulsr(dP*dP, normlp(Br, 2*tnew, dnf)));
       bmin = logint(ceil_safe(mpsqrt(Btra)), gdeux, NULL);
       if (DEBUGLEVEL>2)
         fprintferr("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
                    r, tmax, bmin);
   
     polred=gcopy(polbase);      /* compute Newton sums (possibly relifting first) */
     lt=(GEN)leading_term(polbase)[1];      if (gcmp(GSmin, Btra) < 0)
     lt=mpinvmod(lt,ap);      {
         nflift_t L1;
         GEN polred;
   
     for (i=2; i<d; i++)        bestlift_init(k<<1, nf, T->pr, Btra, &L1);
       polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), aprh)[1];        k      = L1.k;
         pk     = L1.pk;
         PRK    = L1.prk;
         PRKinv = L1.iprk;
         GSmin  = L1.GSmin;
   
     if (!divise(discsr(polred),ap))        polred = ZpX(Pbase, pk, L1.ZpProj);
         famod = hensel_lift_fact(polred,famod,NULL,p,pk,k);
         for (i=1; i<=n0; i++) TT[i] = 0;
       }
       for (i=1; i<=n0; i++)
     {      {
       rep=(GEN)simplefactmod(polred,ap)[1];        GEN p1, lPpow;
       anbf=lg(rep)-1;        GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], tnew, NULL, pk);
       ct--;        TT[i] = (long)p2;
       if (anbf < nbf)        p1 = (GEN)p2[tnew+1];
         /* make Newton sums integral */
         if (lP)
       {        {
         nbf=anbf;          lPpow = gpowgs(lP, tnew);
         pr=gcopy(apr);          p1 = mulii(p1, lPpow); /* assume pr has degree 1 */
         p=gcopy(ap);  
         if (DEBUGLEVEL>=4)  
         {  
           fprintferr("prime ideal considered: %Z\n", pr);  
           fprintferr("number of irreducible factors: %ld\n", nbf);  
         }  
         if (nbf == 1) break;  
       }        }
         if (dn) p1 = mulii(p1,dn);
         if (dn || lP) p1 = modii(p1, pk);
         Tra[i] = (long)nf_bestlift(p1, NULL, L); /* S_tnew(famod) */
     }      }
     if (pr && !ct) break;  
   
     minp = addis(ap,1);      /* compute truncation parameter */
   }      if (DEBUGLEVEL>2) { TIMERstart(&ti2); TIMERstart(&TI); }
       oldCM_L = CM_L;
       av2 = avma;
       b = delta = 0; /* -Wall */
   AGAIN:
       M_L = gdivexact(CM_L, stoi(C));
       T2 = gmul(Tra, M_L);
       VV = gdivround(gmul(PRKinv, T2), pk);
       T2 = gsub(T2, gmul(PRK, VV));
   
   k = itos(gceil(gdiv(k2,glog(p,BIGDEFAULTPREC))));      if (first)
       { /* initialize lattice, using few p-adic digits for traces */
         a = gexpo(T2);
         bgood = (long)(a - max(32, BitPerFactor * r));
         b = max(bmin, bgood);
         delta = a - b;
       }
       else
       { /* add more p-adic digits and continue reduction */
         long b0 = gexpo(T2);
         if (b0 < b) b = b0;
         b = max(b-delta, bmin);
         if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
       }
   
   if (DEBUGLEVEL>=4)      /* restart with truncated entries */
   {      q = shifti(gun, b);
     fprintferr("prime ideal chosen: %Z\n", pr);      P1 = gdivround(PRK, q);
     msgtimer("choice of the prime ideal");      S1 = gdivround(Tra, q);
   }      T2 = gsub(gmul(S1, M_L), gmul(P1, VV));
       m = vconcat( CM_L, T2 );
       if (first)
       {
         first = 0;
         m = concatsp( m, vconcat(ZERO, P1) );
         /*     [ C M_L   0  ]
          * m = [            ]   square matrix
          *     [  T2'   PRK ]   T2' = Tra * M_L  truncated
          */
       }
       CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti, &ti_LLL);
       if (DEBUGLEVEL>2)
         fprintferr("LLL_cmbf: b =%4ld; r =%3ld -->%3ld, time = %ld\n",
                    b, lg(m)-1, CM_L? lg(CM_L)-1: 1, TIMER(&TI));
       if (!CM_L) { list = _col(QXQ_normalize(P,nfT)); break; }
       if (b > bmin)
       {
         CM_L = gerepilecopy(av2, CM_L);
         goto AGAIN;
       }
       if (DEBUGLEVEL>2) msgTIMER(&ti2, "for this trace");
   
   if (lg(rep)==2)      i = lg(CM_L) - 1;
   {      if (i == r && gegal(CM_L, oldCM_L))
     if (fl) { avma=av; return cgetg(1,t_VEC); }      {
     rep=cgetg(2,t_VEC); rep[1]=lcopy(polmod);        CM_L = oldCM_L;
     nfcmbf.nfact=1; return gerepileupto(av, rep);        avma = av2; continue;
       }
   
       if (i <= r && i*rec < n0)
       {
         if (DEBUGLEVEL>2) (void)TIMER(&ti);
         list = nf_check_factors(T, P, M_L, famod, pk);
         if (DEBUGLEVEL>2) ti_CF += TIMER(&ti);
         if (list) break;
         CM_L = gerepilecopy(av2, CM_L);
       }
       if (low_stack(lim, stack_lim(av,1)))
       {
         if(DEBUGMEM>1) err(warnmem,"nf_LLL_cmbf");
         gerepileall(av, 8, &CM_L,&TT,&Tra,&famod,&pk,&GSmin,&PRK,&PRKinv);
       }
   }    }
     if (DEBUGLEVEL>2)
       fprintferr("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
     return list;
   }
   
   p2=cgetr(DEFAULTPREC);  static GEN
   affir(mulii(absi(dk),gpowgs(p,k)),p2);  nf_combine_factors(nfcmbf_t *T, GEN polred, GEN p, long a, long klim)
   p2=shifti(gceil(mplog(p2)),-1);  {
     GEN z, res, L, listmod, famod = T->fact, nf = T->nf;
     long i, m, l, maxK = 3, nft = lg(famod)-1;
   
   newprec = MEDDEFAULTPREC + (long)(itos(p2)*pariK1);    T->fact = hensel_lift_fact(polred,famod,NULL,p,T->pk,a);
   if (DEBUGLEVEL>=4)    if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */
     fprintferr("new precision: %ld\n",newprec);    if (DEBUGLEVEL>3) msgtimer("Hensel lift");
   
   prh = idealpows(nf,pr,k); m = k;    /* FIXME: neither nfcmbf nor LLL_cmbf can handle the non-nf case */
   h = T2_matrix_pow(nf,prh,p,C, &k, newprec);  
   if (m != k) prh=idealpows(nf,pr,k);  
   
   if (DEBUGLEVEL>=4)    T->res      = cgetg(nft+1,t_VEC);
     L = nfcmbf(T, p, a, maxK, klim);
   
     res     = (GEN)L[1];
     listmod = (GEN)L[2]; l = lg(listmod)-1;
     famod = (GEN)listmod[l];
     if (maxK >= 0 && lg(famod)-1 > 2*maxK)
   {    {
     fprintferr("a suitable exponent is: %ld\n",(long)k);      if (l > 1) T->polbase = unifpol(nf, (GEN)res[l], 0);
     msgtimer("computation of H");      L = nf_LLL_cmbf(T, p, a, maxK);
       /* remove last elt, possibly unfactored. Add all new ones. */
       setlg(res, l); res = concatsp(res, L);
   }    }
     if (DEBUGLEVEL>3) msgtimer("computation of the factors");
   
   pk = gcoeff(prh,1,1);    m = lg(res); z = cgetg(m, t_VEC);
   lt=(GEN)leading_term(polbase)[1];    for (i=1;i<m; i++) z[i] = (long)unifpol(nf,(GEN)res[i],1);
   lt=mpinvmod(lt,pk);    return z;
   }
   
   polred[1] = polbase[1];  /* return the factorization of the square-free polynomial x.
   for (i=2; i<d; i++)     The coeff of x are in Z_nf and its leading term is a rational integer.
     polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), prh)[1];     deg(x) > 1, deg(nfpol) > 1
      If fl = 1, return only the roots of x in nf
      If fl = 2, as fl=1 if pol splits, [] otherwise */
   static GEN
   nfsqff(GEN nf, GEN pol, long fl)
   {
     long i, m, n, nbf, ct, dpol = degpol(pol);
     gpmem_t av = avma;
     GEN pr, C0, C, dk, bad, polbase, pk;
     GEN N2, rep, p, ap, polmod, polred, lt, nfpol;
     byteptr pt=diffptr;
     nfcmbf_t T;
     nflift_t L;
   
   fact = lift_intern((GEN)factmod(polred,p)[1]);    if (DEBUGLEVEL>3) msgtimer("square-free");
   rep = hensel_lift_fact(polred,fact,NULL,p,pk,k);    nfpol = (GEN)nf[1]; n = degpol(nfpol);
     polbase = unifpol(nf,pol,0);
     polmod  = unifpol(nf,pol,1);
     /* heuristic */
   #if 1
     if (dpol*4 < n)
     {
       if (DEBUGLEVEL>2) fprintferr("Using Trager's method\n");
       return gerepilecopy(av, (GEN)polfnf(polmod, nfpol)[1]);
     }
   #endif
   
   if (DEBUGLEVEL >= 4) msgtimer("computation of the p-adic factorization");    pol = simplify_i(lift(polmod));
     lt  = (GEN)leading_term(polbase)[1]; /* t_INT */
   
   den = det(h); /* |den| = NP^k */    dk = absi((GEN)nf[3]);
   hinv= adj(h);    bad = mulii(mulii(dk,(GEN)nf[4]), lt);
   lt=(GEN)leading_term(polbase)[1];  
   
   if (fl)    p = polred = pr = NULL; /* gcc -Wall */
     nbf = 0; ap = NULL;
     for (ct = 5;;)
   {    {
     long x_a[] = { evaltyp(t_POL)|m_evallg(4), 0,0,0 };      GEN apr = choose_prime(nf, bad, &ap, &pt);
     GEN mlt = gneg_i(lt), rem;      GEN modpr = zkmodprinit(nf, apr);
     x_a[1] = polbase[1]; setlgef(x_a, 4);      long anbf;
     x_a[3] = un;  
     p1=cgetg(lg(rep)+1,t_VEC);  
     polbase = unifpol(nf,polbase,1);  
     for (m=1,i=1; i<lg(rep); i++)  
     {  
       p2=(GEN)rep[i]; if (lgef(p2)!=4) break;  
   
       p3 = algtobasis(nf, gmul(mlt,(GEN)p2[2]));      polred = modprX(polbase, nf, modpr);
       p3 = nf_bestlift(h,hinv,den,p3);      if (!FpX_is_squarefree(polred,ap)) continue;
       p3 = gdiv(p3,lt);  
       x_a[2] = lneg(p3); /* check P(p3) == 0 */      anbf = fl? FpX_nbroots(polred,ap): FpX_nbfact(polred,ap);
       p2 = poldivres(polbase, unifpol(nf,x_a,1), &rem);      if (!nbf || anbf < nbf)
       if (!signe(rem)) { polbase = p2; p1[m++] = (long)p3; }      {
         nbf = anbf; pr = apr; p = ap;
         if (DEBUGLEVEL>3)
           fprintferr("%3ld %s at prime ideal above %Z\n",
                      nbf, fl?"roots": "factors", p);
         if (fl == 2 && nbf < dpol) return cgetg(1,t_VEC);
         if (nbf <= 1)
         {
           if (!fl) /* irreducible */
             return gerepilecopy(av, _vec(QXQ_normalize(polmod, nfpol)));
           if (!nbf) return cgetg(1,t_VEC); /* no root */
         }
     }      }
     tetpil=avma; rep=cgetg(m,t_VEC);      if (--ct <= 0) break;
     for (i=1; i<m; i++) rep[i]=(long)basistoalg(nf,(GEN)p1[i]);  
     return gerepile(av,tetpil,rep);  
   }    }
     if (DEBUGLEVEL>3) {
       msgtimer("choice of a prime ideal");
       fprintferr("Prime ideal chosen: %Z\n", pr);
     }
   
   for (i=1; i<lg(rep); i++)    L.tozk = (GEN)nf[8];
     rep[i] = (long)unifpol(nf,(GEN)rep[i],0);    L.topow= (GEN)nf[7];
     T.ZC = L2_bound(nf, L.tozk, &(T.dn));
     T.Br = gmul(lt, nf_root_bounds(pol, nf));
   
   nfcmbf.pol      = polmod;    if (fl) C0 = normlp(T.Br, 2, n);
   nfcmbf.lt       = lt;    else    C0 = nf_factor_bound(nf, polbase); /* bound for T_2(Q_i), Q | P */
   nfcmbf.h        = h;    C = mulrr(T.ZC, C0); /* bound for |Q_i|^2 in Z^n on chosen Z-basis */
   nfcmbf.hinv     = hinv;    T.bound = C;
   nfcmbf.den      = den;  
   nfcmbf.fact     = rep;  
   nfcmbf.res      = cgetg(lg(rep)+1,t_VEC);  
   nfcmbf.nfact    = 0;  
   nfcmbf.nfactmod = lg(rep)-1;  
   nf_combine_factors(nf,1,NULL,d-3,1);  
   
   if (DEBUGLEVEL >= 4) msgtimer("computation of the factors");    N2 = mulsr(dpol*dpol, normlp(T.Br, 4, n)); /* bound for T_2(lt * S_2) */
     T.BS_2 = mulrr(T.ZC, N2); /* bound for |S_2|^2 on chosen Z-basis */
   
   i = nfcmbf.nfact;    if (DEBUGLEVEL>2) {
   if (lgef(nfcmbf.pol)>3)      msgtimer("bound computation");
   {      fprintferr("  1) T_2 bound for %s: %Z\n", fl?"root":"factor", C0);
     nfcmbf.res[++i] = (long) nf_pol_divres(nf,nfcmbf.pol,nfcmbf.lt,NULL);      fprintferr("  2) Conversion from T_2 --> | |^2 bound : %Z\n", T.ZC);
     nfcmbf.nfact = i;      fprintferr("  3) Final bound: %Z\n", C);
   }    }
   
   tetpil=avma; rep=cgetg(i+1,t_VEC);    if (fl)
   for (  ; i>=1; i--)      (void)logint(C, p, &pk);
     rep[i]=(long)unifpol(nf,(GEN)nfcmbf.res[i],1);  #if 0 /* overkill */
   return gerepile(av,tetpil,rep);    else
 }    { /* overlift needed for d-1/d-2 tests? */
       GEN pb; long b; /* junk */
 static int      if (cmbf_precs(p, C, T.BS_2, &a, &b, &pk, &pb))
 nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint)      { /* Rare */
 {        if (DEBUGLEVEL) err(warner,"nffactor: overlift for d-1/d-2 test");
   int val = 0; /* assume failure */        C = itor(pk, DEFAULTPREC);
   GEN newf, newpsf = NULL;      }
   long av,newd,ltop,i;  
   
   /* Assertion: fxn <= nfcmbf.nfactmod && dlim > 0 */  
   
   /* first, try deeper factors without considering the current one */  
   if (fxn != nfcmbf.nfactmod)  
   {  
     val = nf_combine_factors(nf,fxn+1,psf,dlim,hint);  
     if (val && psf) return 1;  
   }    }
   #endif
   
   /* second, try including the current modular factor in the product */    bestlift_init(0, nf, pr, C, &L);
   newf = (GEN)nfcmbf.fact[fxn];    T.pr = pr;
   if (!newf) return val; /* modular factor already used */    T.L  = &L;
   newd = degpol(newf);  
   if (newd > dlim) return val; /* degree of new factor is too large */  
   
   av = avma;    /* polred is monic */
   if (newd % hint == 0)    pk = L.pk;
   {    polred = ZpX(gmul(mpinvmod(lt,pk), polbase), pk, L.ZpProj);
     GEN p, quot,rem;  
   
     newpsf = nf_pol_mul(nf, psf? psf: nfcmbf.lt, newf);    if (fl)
     newpsf = nf_pol_lift(nfcmbf.h,nfcmbf.hinv,nfcmbf.den,newpsf);    { /* roots only */
     /* try out the new combination */      long x_r[] = { evaltyp(t_POL)|_evallg(4), 0,0,0 };
     ltop=avma;      rep = rootpadicfast(polred, p, L.k);
     quot=nf_pol_divres(nf,nfcmbf.pol,newpsf,&rem);      x_r[1] = evalsigne(1) | evalvarn(varn(pol)) | evallgef(4);
     if (gcmp0(rem))  /* found a factor */      x_r[3] = un;
       for (m=1,i=1; i<lg(rep); i++)
     {      {
       p = nf_pol_mul(nf,element_inv(nf,leading_term(newpsf)),newpsf);        GEN q, r = (GEN)rep[i];
       nfcmbf.res[++nfcmbf.nfact] = (long) p; /* store factor */  
       nfcmbf.fact[fxn]=0;                    /* remove used modular factor */  
   
       /* fix up target */        r = nf_bestlift_to_pol(gmul(lt,r), NULL, &L);
       p=gun; quot=unifpol(nf,quot,0);        r = gdiv(r,lt);
       for (i=2; i<lgef(quot); i++)        x_r[2] = lneg(r); /* check P(r) == 0 */
         p = mpppcm(p, denom((GEN)quot[i]));        q = RXQX_divrem(pol, x_r, nfpol, ONLY_DIVIDES);
         if (q) { pol = q; rep[m++] = (long)r; }
       nfcmbf.pol = nf_pol_mul(nf,p,quot);        else if (fl == 2) return cgetg(1, t_VEC);
       nfcmbf.lt  = leading_term(nfcmbf.pol);  
       return 1;  
     }      }
     avma=ltop;      rep[0] = evaltyp(t_VEC) | evallg(m);
       return gerepilecopy(av, rep);
   }    }
   
   /* If room in degree limit + more modular factors to try, add more factors to    T.polbase = polbase;
    * newpsf */    T.pol   = pol;
   if (newd < dlim && fxn < nfcmbf.nfactmod    T.nf    = nf;
                   && nf_combine_factors(nf,fxn+1,newpsf,dlim-newd,hint))    T.fact  = (GEN)factmod0(polred,p)[1];
   {    T.hint  = 1; /* useless */
     nfcmbf.fact[fxn]=0; /* remove used modular factor */  
     return 1;    rep = nf_combine_factors(&T, polred, p, L.k, dpol-1);
   }    return gerepileupto(av, rep);
   avma = av; return val;  
 }  }
   
 /* return the characteristic polynomial of alpha over nf, where alpha  /* return the characteristic polynomial of alpha over nf, where alpha
    is an element of the algebra nf[X]/(T) given as a polynomial in X */     is an element of the algebra nf[X]/(T) given as a polynomial in X */
 GEN  GEN
 rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)  rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
 {  {
   long av = avma, vnf, vT, lT;    long vnf, vT, lT;
     gpmem_t av = avma;
   GEN p1;    GEN p1;
   
   nf=checknf(nf); vnf = varn(nf[1]);    nf=checknf(nf); vnf = varn(nf[1]);
Line 1243  rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
Line 1414  rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
     return gerepileupto(av, gsub(polx[v], alpha));      return gerepileupto(av, gsub(polx[v], alpha));
   p1 = caract2(unifpol(nf,T,1), unifpol(nf,alpha,1), v);    p1 = caract2(unifpol(nf,T,1), unifpol(nf,alpha,1), v);
   return gerepileupto(av, unifpol(nf,p1,1));    return gerepileupto(av, unifpol(nf,p1,1));
 }  
   
 #if 0  
 /* return the minimal polynomial of alpha over nf, where alpha is  
    an element of the algebra nf[X]/(T) given as a polynomial in X */  
 GEN  
 rnfminpoly(GEN nf,GEN T,GEN alpha,int n)  
 {  
   long av=avma,tetpil;  
   GEN p1,p2;  
   
   nf=checknf(nf); p1=rnfcharpoly(nf,T,alpha,n);  
   tetpil=avma; p2=nf_pol_subres(nf,p1,deriv(p1,varn(T)));  
   if (lgef(p2)==3) { avma=tetpil; return p1; }  
   
   p1 = nf_pol_divres(nf,p1,p2,NULL);  
   p2 = element_inv(nf,leading_term(p1));  
   tetpil=avma; return gerepile(av,tetpil,unifpol(nf,nf_pol_mul(nf,p2,p1),1));  
 }  
 #endif  
   
 /* relative Dedekind criterion over nf, applied to the order defined by a  
  * root of irreducible polynomial T, 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 in pr of  
  * the order discriminant  
  */  
 GEN  
 rnfdedekind(GEN nf,GEN T,GEN pr)  
 {  
   long av=avma,vt,r,d,da,n,m,i,j;  
   GEN p1,p2,p,tau,g,vecun,veczero,matid;  
   GEN prhall,res,h,k,base,Ca;  
   
   nf=checknf(nf); Ca=unifpol(nf,T,0);  
   res=cgetg(4,t_VEC);  
   if (typ(pr)==t_VEC && lg(pr)==3)  
   { prhall = (GEN)pr[2]; pr = (GEN)pr[1]; }  
   else  
     prhall=nfmodprinit(nf,pr);  
   p=(GEN)pr[1]; tau=gdiv((GEN)pr[5],p);  
   n=degpol(nf[1]); m=degpol(T);  
   
   vecun = gscalcol_i(gun,n);  
   veczero = zerocol(n);  
   
   p1=(GEN)nffactormod(nf,Ca,pr)[1];  
   r=lg(p1); if (r < 2) err(constpoler,"rnfdedekind");  
   g=lift((GEN)p1[1]);  
   for (i=2; i<r; i++)  
     g = nf_pol_mul(nf,g, lift((GEN)p1[i]));  
   h=nfmod_pol_divres(nf,prhall,Ca,g,NULL);  
   k=nf_pol_mul(nf,tau,gsub(Ca, nf_pol_mul(nf,lift(g),lift(h))));  
   p2=nfmod_pol_gcd(nf,prhall,g,h);  
   k= nfmod_pol_gcd(nf,prhall,p2,k);  
   
   d=degpol(k);  /* <= m */  
   vt = idealval(nf,discsr(T),pr) - 2*d;  
   res[3]=lstoi(vt);  
   if (!d || vt<=1) res[1]=un; else res[1]=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++)  
   {  
     p2[j]=(long)matid;  
     p1[j]=lgetg(m+1,t_COL);  
     for (i=1; i<=m; i++)  
       coeff(p1,i,j) = (i==j)?(long)vecun:(long)veczero;  
   }  
   if (d)  
   {  
     GEN pal = lift(nfmod_pol_divres(nf,prhall,Ca,k,NULL));  
     GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */  
     GEN nfx = unifpol(nf,polx[varn(T)],0);  
   
     for (   ; j<=m+d; j++)  
     {  
       p1[j]=lgetg(m+1,t_COL);  
       da=degpol(pal);  
       for (i=1; i<=da+1; i++) coeff(p1,i,j)=pal[i+1];  
       for (   ; i<=m; i++) coeff(p1,i,j)=(long)veczero;  
       p2[j]=(long)prinvp;  
       nf_pol_divres(nf,nf_pol_mul(nf,pal,nfx),T,&pal);  
     }  
     /* 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);  
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>