=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/modules/Attic/nffactor.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/modules/Attic/nffactor.c 2001/10/02 11:17:11 1.1 +++ OpenXM_contrib/pari-2.2/src/modules/Attic/nffactor.c 2002/09/11 07:27:06 1.2 @@ -1,4 +1,4 @@ -/* $Id: nffactor.c,v 1.1 2001/10/02 11:17:11 noro Exp $ +/* $Id: nffactor.c,v 1.2 2002/09/11 07:27:06 noro Exp $ Copyright (C) 2000 The PARI group. @@ -21,39 +21,70 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, #include "pari.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 nf_get_T2(GEN base, GEN polr); -extern GEN nfreducemodpr_i(GEN x, GEN prh); +extern GEN initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis); +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 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 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; - long nfact, nfactmod; -} Nfcmbf; -static Nfcmbf nfcmbf; + nflift_t *L; + GEN nf, pol, fact, res, bound, pr, pk, Br, ZC, dn, polbase, BS_2; + long hint; +} nfcmbf_t; +GEN +RXQX_mul(GEN x, GEN y, GEN T) +{ + return RXQX_red(gmul(x,y), T); +} + static GEN unifpol0(GEN nf,GEN pol,long flag) { static long n = 0; static GEN vun = NULL; GEN f = (GEN) nf[1]; - long i = degpol(f), av; + long i = degpol(f); + gpmem_t av; if (i != n) { @@ -67,9 +98,11 @@ unifpol0(GEN nf,GEN pol,long flag) switch(typ(pol)) { case t_INT: case t_FRAC: case t_RFRAC: + if (flag) return gcopy(pol); pol = gmul(pol,vun); break; case t_POL: + if (flag && !degpol(pol)) return gcopy(constant_term(pol)); pol = gmodulcp(pol,f); /* fall through */ case t_POLMOD: pol = algtobasis(nf,pol); @@ -100,1132 +133,1270 @@ unifpol(GEN nf,GEN pol,long flag) return unifpol0(nf,(GEN) pol, flag); } -#if 0 -/* return a monic polynomial of degree d with random coefficients in Z_nf */ -static GEN -random_pol(GEN nf,long d) +/* factorization of x modulo pr. Assume x already reduced */ +GEN +FqX_factor(GEN x, GEN T, GEN p) { - long i,j, n = degpol(nf[1]); - GEN pl,p; - - pl=cgetg(d+3,t_POL); - for (i=2; i= 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 -nf_pol_sqr(GEN nf,GEN x) +choose_prime(GEN nf, GEN bad, GEN *p, byteptr *PT) { - long tetpil,av=avma; - GEN res = gsqr(unifpol(nf,x,1)); - - tetpil = avma; - return gerepile(av,tetpil,unifpol(nf,res,0)); + GEN q = icopy(gun), r, x = (GEN)nf[1]; + ulong pp = *p? itou(*p): 0; + byteptr pt = *PT; + gpmem_t av = avma; + 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 -nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol) +QXQ_normalize(GEN P, GEN T) { - long av=avma,tetpil,i; - GEN p1; - - if (typ(pol)!=t_POL) return nfreducemodpr(nf,pol,prhall); - pol=unifpol(nf,pol,0); - - tetpil=avma; i=lgef(pol); - p1=cgetg(i,t_POL); p1[1]=pol[1]; - for (i--; i>=2; i--) - p1[i] = (long) nfreducemodpr(nf,(GEN)pol[i],prhall); - return gerepile(av,tetpil, normalizepol(p1)); + GEN t = leading_term(P); + if (!gcmp1(t)) + { + if (is_rational_t(typ(t))) + P = gdiv(P, t); + else + P = RXQX_mul(QX_invmod(t,T), P, T); + } + return P; } -/* x^2 modulo prhall */ -static GEN -nfmod_pol_sqr(GEN nf,GEN prhall,GEN x) +/* return the roots of pol in nf */ +GEN +nfroots(GEN nf,GEN pol) { - long av=avma,tetpil; - GEN px; + gpmem_t av = avma; + int d = degpol(pol); + GEN A,g, T; - px = nfmod_pol_reduce(nf,prhall,x); - px = unifpol(nf,lift(px),1); - px = unifpol(nf,nf_pol_sqr(nf,px),0); - tetpil=avma; - return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px)); + nf = checknf(nf); T = (GEN)nf[1]; + if (typ(pol) != t_POL) err(notpoler,"nfroots"); + if (varn(pol) >= varn(T)) + err(talker,"polynomial variable must have highest priority in nfroots"); + 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 */ -static GEN -nfmod_pol_mul(GEN nf,GEN prhall,GEN x,GEN y) +int +nfissplit(GEN nf, GEN x) { - long av=avma,tetpil; - GEN px,py; - - 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)); + gpmem_t av = avma; + long l = lg(nfsqff(checknf(nf), x, 2)); + avma = av; return l != 1; } -/* Euclidan division of x by y */ -static GEN -nf_pol_divres(GEN nf,GEN x,GEN y,GEN *pr) +/* nf = K[y] / (P). Galois over K ? */ +int +nfisgalois(GEN nf, GEN P) { - long av = avma,tetpil; - 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; + return degpol(P) <= 2 || nfissplit(nf, P); } -/* Euclidan division of x by y modulo prhall */ +/* return a minimal lift of elt modulo id */ 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 z,p1,p3,px,py; - - 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 (dxprk), t = typ(elt); + if (t != t_INT) { - GEN vzero; - - if (pr) *pr = gerepile(av,tetpil,px); - 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; + if (t == t_POL) elt = mulmat_pol(T->tozk, elt); + u = gmul(T->iprk,elt); + for (i=1; ipk); } - p1 = NULL; /* gcc -Wall */ - - 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) + else { - l=avma; p1=(GEN)px[i+2]; - for (j=i-dy+1; j<=i && j<=dz; j++) - p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2])); - 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); + u = gmul(elt, (GEN)T->iprk[1]); + for (i=1; ipk); + elt = gscalcol(elt, l-1); } - l=avma; - for (i=dy-1; i>=0; --i) - { - 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; + u = gsub(elt, gmul(T->prk, u)); + if (bound && gcmp(QuickNormL2(u,DEFAULTPREC), bound) > 0) u = NULL; + return u; } -/* GCD of x and y */ 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 s = srgcd(unifpol(nf,x,1), unifpol(nf,y,1)); - - tetpil=avma; return gerepile(av,tetpil,unifpol(nf,s,1)); + GEN u = nf_bestlift(elt,bound,T); + if (u) u = gmul(T->topow, u); + return u; } -/* GCD of x and y modulo prhall */ +/* return the lift of pol with coefficients of T2-norm <= C (if possible) */ 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; - GEN p1,p2; + long i, d = lgef(pol); + GEN x = cgetg(d,t_POL); + nflift_t *L = T->L; - if (lgef(x)3) (void)timer2(); - vun=cgetg(n+1,t_COL); vun[1]=un; for (i=2; i<=n; i++) vun[i]=zero; - p1=gcopy(polun[varn(pol)]); p1[2]=(long)vun; - if (gcmp0(e)) return p1; + nf = checknf(nf); T = (GEN)nf[1]; + if (typ(pol) != t_POL) err(notpoler,"nffactor"); + if (varn(pol) >= varn(T)) + err(talker,"polynomial variable must have highest priority in nffactor"); - p2=nfmod_pol_reduce(nf,prhall,pol); - for(;;) + if (d == 0) return trivfact(); + rep = cgetg(3, t_MAT); av = avma; + if (d == 1) { - if (!vali(e)) - { - p1=nfmod_pol_mul(nf,prhall,p1,p2); - 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); + rep[1] = (long)_col( gcopy(pol) ); + rep[2] = (long)_col( gun ); + return rep; } - return gerepileupto(av,p1); -} -static long -isdivbyprime(GEN nf, GEN x, GEN pr) -{ - GEN elt, p = (GEN)pr[1], tau = (GEN)pr[5]; + A = fix_relative_pol(nf,pol,0); + if (degpol(nf[1]) == 1) + return gerepileupto(av, factpol(simplify(pol), 0)); - elt = element_mul(nf, x, tau); - if (divise(content(elt), p)) return 1; + A = primpart( lift_intern(A) ); + 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 */ -static GEN -localpol(GEN nf, GEN pr) -{ - long i, l; - GEN fct, pol = (GEN)nf[1], p = (GEN)pr[1]; + if (degpol(g)) + { /* not squarefree */ + gpmem_t av1; + GEN ex; + g = QXQ_normalize(g, T); + A = RXQX_div(A,g, T); - fct = lift((GEN)factmod(pol, p)[1]); - l = lg(fct) - 1; - for (i = 1; i <= l; i++) - if (isdivbyprime(nf, (GEN)fct[i], pr)) return (GEN)fct[i]; - - err(talker, "cannot find a suitable factor in localpol"); - return NULL; /* not reached */ + y = nfsqff(nf,A,0); av1 = avma; + l = lg(y); + ex=(GEN)gpmalloc(l * sizeof(long)); + for (j=l-1; j>=1; j--) + { + GEN fact = lift((GEN)y[j]), quo = g, q; + 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 -nffactormod0(GEN nf, GEN x, GEN pr) +nf_Mignotte_bound(GEN nf, GEN polbase) { - long av = avma, j, l, vx = varn(x), vn; - GEN rep, pol, xrd, prh, p1; + GEN G = gmael(nf,5,2), lS = leading_term(polbase); /* t_INT */ + 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); - 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"); + if (typ(lS) == t_COL) lS = (GEN)lS[1]; + binlS = bin = vecbinome(d-1); + if (!gcmp1(lS)) binlS = gmul(lS, bin); - 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])) - return gerepileupto(av, nffactormod2(nf,x,pr)); + matGS = cgetg(d+2, t_MAT); + 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); - if (gcmp1((GEN)pr[4])) + /* Take sup over 0 <= i <= d of + * 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 */ - for (j = 2; j < lg(xrd); 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++) + GEN s = gzero; + for (j = 1; j <= n; j++) { - xrd[j] = lmod((GEN)xrd[j], pol); - p1 = mpppcm(p1, denom(content((GEN)xrd[j]))); + p1 = mpadd( mpmul((GEN)bin[i], (GEN)N2[j]), (GEN)binlS[i+1] ); + s = mpadd(s, gsqr(p1)); } - rep = factmod9(gmul(xrd, p1), (GEN)pr[1], pol); - rep = lift(lift(rep)); + if (gcmp(C, s) < 0) C = s; } + return C; +} - l = lg((GEN)rep[1]); - for (j = 1; j < l; j++) - mael(rep, 1, j) = (long)unifpol(nf, gmael(rep, 1, j), 1); +/* return a bound for T_2(P), P | polbase + * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2, + * 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 -nffactormod(GEN nf, GEN x, GEN pr) +static GEN +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); - -/* factorization of x modulo pr */ -GEN -nffactormod2(GEN nf,GEN pol,GEN pr) +static long +ZXY_get_prec(GEN P) { - long av = avma, tetpil,lb,nbfact,psim,N,n,i,j,k,d,e,vf,r,kk; - GEN y,ex,*t,f1,f2,f3,df1,g1,polb,pold,polu,vker; - 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 + long i, j, z, prec = 0; + for (i=2; i3) + GEN p = (GEN)P[i]; + if (typ(p) == t_INT) { - df1=deriv(f1,vf); f2=nfmod_pol_gcd(nf,prhall,f1,df1); - g1=nfmod_pol_divres(nf,prhall,f1,f2,NULL); k=0; - while (lgef(g1)>3) + z = lgefint(p); + if (z > prec) prec = z; + } + else + { + for (j=2; j3) - { - 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 (j1) - { - 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>8, q); - vz=nfreducemodpr(nf,vz,prhall); - v=gadd(v,nfmod_pol_mul(nf,prhall,vz,(GEN)vker[i])); - } - for (i=1; i<=kk && kk4) - { - 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) prec) prec = z; } - e*=psim; j=(degpol(f2))/psim+3; f1=cgetg(j,t_POL); - f1[1] = evalsigne(1) | evallgef(j) | evalvarn(vf); - for (i=2; i k) k = l; } } - - 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 k) k = l; } - 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 -p_ok(GEN nf, GEN p, GEN a) +complex_bound(GEN P) { - long av,m,i; - GEN dec; - - if (divise(a,p)) return NULL; - av = avma; dec = primedec(nf,p); m=lg(dec); - for (i=1; i= 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 polbase is square-free) - and 0 otherwise (=> polbase may or may not be square-free) */ -static int -is_sqf(GEN nf, GEN polbase) +/* return B such that if x in O_K, K = Z[X]/(T), then the L2-norm of the + * coordinates of the numerator of x [on the power, resp. integral, basis if T + * is a polynomial, resp. an nf] is <= B T_2(x) + * *ptden = multiplicative bound for denom(x) */ +static GEN +L2_bound(GEN T, GEN tozk, GEN *ptden) { - GEN lt, pr, prh, p2, p; - long i, d = lgef(polbase), ct = 5; + GEN nf, M, L, prep, den, u; + long prec; - lt = (GEN)leading_term(polbase)[1]; - p = stoi(101); + T = get_nfpol(T, &nf); + u = NULL; /* gcc -Wall */ - while (ct > 0) - { - /* small primes tend to divide discriminants more often - than large ones so we look at primes >= 101 */ - pr = choose_prime(nf,lt,p,30); - if (!pr) break; + prec = ZX_get_prec(T) + ZM_get_prec(tozk); + den = nf? gun: NULL; + den = initgaloisborne(T, den, prec, &L, &prep, NULL); + M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep); + if (nf) M = gmul(tozk, M); + if (gcmp1(den)) den = NULL; + *ptden = den; + return QuickNormL2(M, DEFAULTPREC); +} - p=(GEN)pr[1]; - prh=prime_to_ideal(nf,pr); +/* || L ||_p^p in dimension n (L may be a scalar) */ +static GEN +normlp(GEN L, long p, long n) +{ + long i,l, t = typ(L); + GEN z; - p2=gcopy(polbase); - lt=mpinvmod(lt,p); + if (!is_vec_t(t)) return gmulsg(n, gpowgs(L, p)); - for (i=2; i polynomial is square-free */ - if (!gcmp0(p2) && !divise(discsr(p2),p)) { return 1; } - - ct--; - p=addis(p,1); - } - - return 0; + l = lg(L); z = gzero; + /* assert(n == l-1); */ + for (i=1; i primitive in O_K[X] */ -GEN -nf_pol_to_int(GEN p, GEN *den) +/* S = S0 + qS1, P = P0 + qP1 (Euclidean div. by q integer). For a true + * factor (vS, vP), we have: + * | 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); - GEN d = gun; - for (i=2; i30) + 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 */ -GEN -nfroots(GEN nf,GEN pol) +typedef struct { + GEN d; + 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; - GEN p1,p2,polbase,polmod,den; + long i, j, l, K = lg(ind)-1; + GEN z, s, v; - p2=NULL; /* gcc -Wall */ - nf=checknf(nf); - 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"); + s = (GEN)T->S1[ ind[1] ]; + if (K == 1) return s; - 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; iPinvSdbl[ 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) - { - tetpil=avma; p1=cgetg(2,t_VEC); - p1[1] = (long)basistoalg(nf,gneg_i( - element_div(nf,(GEN)polbase[2],(GEN)polbase[3]))); - return gerepile(av,tetpil,p1); - } +static trace_data * +init_trace(trace_data *T, GEN S, nflift_t *L, GEN q) +{ + long e = gexpo((GEN)S), i,j, l,h; + GEN qgood, S1, invd; - p1=element_inv(nf,leading_term(polbase)); - polbase=nf_pol_mul(nf,p1,polbase); + if (e < 0) return NULL; /* S = 0 */ - polbase = nf_pol_to_int(polbase, &den); - polmod=unifpol(nf,polbase,1); + qgood = shifti(gun, e - 32); /* single precision check */ + if (cmpii(qgood, q) > 0) q = qgood; - if (DEBUGLEVEL>=4) - fprintferr("test if the polynomial is square-free\n"); + S1 = gdivround(S, q); + if (gcmp0(S1)) return NULL; - fl = is_sqf(nf, polbase); + invd = ginv(itor(L->den, DEFAULTPREC)); - /* the polynomial may not be square-free ... */ - if (!fl) + T->dPinvS = gmul(L->iprk, S); + 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); - p2=nf_pol_subres(nf,polmod,p1); - if (degpol(p2) == 0) fl = 1; + double *t = dalloc(h * sizeof(double)); + GEN c = (GEN)T->dPinvS[j]; + T->PinvSdbl[j] = t; + for (i=1; i < h; i++) t[i] = rtodbl(gmul(invd, (GEN)c[i])); } - if (!fl) - { - p1=element_inv(nf,leading_term(p2)); - 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)); + T->d = L->den; + T->P1 = gdivround(L->prk, q); + T->S1 = S1; return T; } -/* return a minimal lift of elt modulo id */ -static GEN -nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt) +static void +update_trace(trace_data *T, long k, long i) { - GEN u = gmul(idinv,elt); - long i, l = lg(u); - for (i=1; iS1[k] = T->S1[i]; + T->dPinvS[k] = T->dPinvS[i]; + 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 -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); - GEN x = cgetg(d,t_POL); - x[1] = pol[1]; - for (i=2; ipol, nf = T->nf, famod = T->fact; + GEN bound = T->bound; + GEN nfpol = (GEN)nf[1]; + long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol); + GEN lc, lcpol; + GEN pk = gpowgs(p,a), pas2 = shifti(pk,-1); -#if 0 -/* return pol(elt) */ -static GEN -nf_pol_eval(GEN nf,GEN pol,GEN elt) -{ - long av=avma,tetpil,i; - GEN p1; + GEN trace1 = cgetg(lfamod+1, t_MAT); + GEN trace2 = cgetg(lfamod+1, t_MAT); + GEN ind = cgetg(lfamod+1, t_VECSMALL); + GEN degpol = cgetg(lfamod+1, t_VECSMALL); + GEN degsofar = cgetg(lfamod+1, t_VECSMALL); + GEN listmod = cgetg(lfamod+1, t_COL); + 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); - for (i--; i>=3; i--) - p1=element_mul(nf,elt,gadd((GEN)pol[i],p1)); - tetpil=avma; return gerepile(av,tetpil,gadd(p1,(GEN)pol[2])); -} -#endif + lc = absi(leading_term(pol)); + if (gcmp1(lc)) lc = NULL; + lcpol = lc? gmul(lc,pol): pol; -/* 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); - 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; - } + GEN t1,t2, lc2 = lc? sqri(lc): NULL; - p1=element_inv(nf,leading_term(pol)); - pol=nf_pol_mul(nf,p1,pol); + for (i=1; i <= lfamod; i++) + { + GEN P = (GEN)famod[i]; + long d = degpol(P); - p1=unifpol(nf,pol,0); - p1 = nf_pol_to_int(p1, &den); + degpol[i] = d; P += 2; + 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) - fprintferr("test if the polynomial is square-free\n"); - - 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; + T1 = init_trace(&_T1, trace1, T->L, q); + T2 = init_trace(&_T2, trace2, T->L, q); } + degsofar[0] = 0; /* sentinel */ - if (!fl) - { - p4=element_inv(nf,leading_term(p3)); - p3=nf_pol_mul(nf,p4,p3); + /* ind runs through strictly increasing sequences of length K, + * 1 <= ind[i] <= lfamod */ +nextK: + 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); - p4=element_inv(nf,leading_term(p2)); - p2=nf_pol_mul(nf,p4,p2); + av = avma; + /* d - 1 test */ + 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); - tetpil = avma; y = nfsqff(nf,p2,0); - i = nfcmbf.nfact; + y = primpart(y); + fa[cnt++] = (long)QXQ_normalize(y, nfpol); + /* 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); - p3=(GEN)gpmalloc((i+1) * sizeof(long)); - for (j=i; j>=1; j--) +NEXT: + for (i = K+1;;) { - GEN fact=(GEN)y[j], quo = quot, rem; - long e = 0; - - do + if (--i == 0) { K++; goto nextK; } + if (++ind[i] <= lfamod - K + i) { - quo = nf_pol_divres(nf,quo,fact,&rem); - e++; + curdeg = degsofar[i-1] + degpol[ind[i]]; + 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 - { - tetpil=avma; y = gerepile(av,tetpil,nfsqff(nf,p1,0)); - i = nfcmbf.nfact; - 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); -} +END: + if (degpol(pol) > 0) + { /* leftover factor */ + if (signe(leading_term(pol)) < 0) pol = gneg_i(pol); -/* test if the matrix M is suitable */ -static long -test_mat(GEN M, GEN p, GEN C2, long k) -{ - 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; + if (lc && lfamod < 2*K) pol = QXQ_normalize(primpart(pol), nfpol); + setlg(famod, lfamod+1); + listmod[cnt] = (long)dummycopy(famod); + fa[cnt++] = (long)pol; } - R = mpmul(min, gpowgs(p, k<<1)); - i = mpcmp(mpmul(C2,prod), R); avma = av; - return (i < 0); + if (DEBUGLEVEL>6) fprintferr("\n"); + setlg(listmod, cnt); setlg(fa, cnt); + res[1] = (long)fa; + res[2] = (long)listmod; return res; } -/* return the matrix corresponding to pr^e with R(pr^e) > C */ 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]); - int tot_real = !signe(gmael(nf,2,2)); - GEN p1,p2,p3,u,C2,T2, x = (GEN)nf[1]; + GEN nf = T->nf, bound = T->bound; + GEN nfT = (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])); - p1 = gmul(pre,lllintpartial(pre)); av1 = avma; - T2 = tot_real? gmael(nf,5,4) - : nf_get_T2((GEN) nf[7], roots(x,prec)); - p3 = qf_base_change(T2,p1,1); + piv = special_pivot(M_L); + if (!piv) return NULL; + if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv); - 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; i2) fprintferr("exponent: %ld\n",k); + y = lc; + for (j=1; j<=n0; j++) + if (signe(c[j])) + { + 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(;;) - { - u = tot_real? lllgramall(p3,2,lll_IM) : lllgramintern(p3,2,1,prec); - if (u) break; + /* y is the candidate factor */ + pol = RXQX_divrem(lcpol,y,nfT, ONLY_DIVIDES); + if (!pol) return NULL; - prec=(prec<<1)-2; - if (DEBUGLEVEL > 1) err(warnprec,"nffactor[1]",prec); - 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)) + y = primpart(y); + if (lc) { - *kmax = k; - return gerepileupto(av,gmul(p1,u)); + pol = primpart(pol); + lc = absi(leading_term(pol)); } - - /* we also need to increase the precision */ - 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); + lcpol = lc? gmul(lc, pol): pol; + list[i] = (long)QXQ_normalize(y, nfT); } + 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 -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; - 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; + return centermodii(gmul(ffproj, x), pk, pas2); +} - 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]); - dki=mulii(dk,(GEN)nf[4]); - n=degpol(nf[1]); + if (typ(P)!=t_POL) return nf_to_Zp(P,pk,pas2,ffproj); + l = lgef(P); + z = cgetg(l,t_POL); z[1] = P[1]; + for (i=2; i= C. Let alpha the LLL parameter, + * 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; - for (;;) - { - if (prec <= gprecision(nf)) - T2 = gprec_w(gmael(nf,5,3), prec); - else - T2 = nf_get_T2((GEN)nf[7], roots((GEN)nf[1], prec)); +static void +bestlift_init(long a, GEN nf, GEN pr, GEN C, nflift_t *T) +{ + const int D = 4; + const double alpha = ((double)D-1) / D; /* LLL parameter */ + gpmem_t av = avma; + GEN prk, PRK, B, GSmin, pk; - run=realun(prec); - p1=realzero(prec); - for (i=2; i 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; - if (DEBUGLEVEL > 1) err(warnprec, "nfsqff", prec); + PRK = lllint_i(PRK, D, 0, NULL, NULL, &B); + 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]; - p1 = gmul(p1, mulis(sqri(lt), n)); - C = gpow(stoi(3), gadd(gmulsg(3, ghalf), stoi(d)), prec); - C = gdiv(gmul(C, p1), gmulsg(d, mppi(prec))); - - if (DEBUGLEVEL>=4) - fprintferr("the bound on the T2-norm of the coeff. is: %Z\n", C); +/* assume pr has degree 1 */ +static GEN +nf_LLL_cmbf(nfcmbf_t *T, GEN p, long k, long rec) +{ + nflift_t *L = T->L; + GEN pk = L->pk, PRK = L->prk, PRKinv = L->iprk, GSmin = L->GSmin; - /* the theoretical bound for the exponent should be: - k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.347))); */ - k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.15))); - k2=gmul2n(gmulgs(k2,n),-1); + GEN famod = T->fact, nf = T->nf, ZC = T->ZC, Br = T->Br; + GEN Pbase = T->polbase, P = T->pol, dn = T->dn; + GEN nfT = (GEN)nf[1]; + GEN Btra; + long dnf = degpol(nfT), dP = degpol(P); - minp=gmin(gexp(gmul2n(k2,-6),BIGDEFAULTPREC), maxp); - minp=gceil(minp); - - if (DEBUGLEVEL>=4) - { - fprintferr("lower bound for the prime numbers: %Z\n", minp); - msgtimer("bounds computation"); - } + double BitPerFactor = 0.5; /* nb bits / modular factor */ + long i, C, tmax, n0; + GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO; + double Blow; + gpmem_t av, av2, lim; + long ti_LLL = 0, ti_CF = 0; + pari_timer ti, ti2, TI; - p = rep = polred = NULL; /* gcc -Wall */ - pr=NULL; - for (;;) + if (DEBUGLEVEL>2) (void)TIMER(&ti); + + 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); - if (!apr) break; + long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1; + GEN oldCM_L, M_L, q, S1, P1, VV; + int first = 1; - ap=(GEN)apr[1]; - aprh=prime_to_ideal(nf,apr); + /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */ + 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); - lt=(GEN)leading_term(polbase)[1]; - lt=mpinvmod(lt,ap); + /* compute Newton sums (possibly relifting first) */ + if (gcmp(GSmin, Btra) < 0) + { + nflift_t L1; + GEN polred; - for (i=2; ipr, Btra, &L1); + 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]; - anbf=lg(rep)-1; - ct--; - if (anbf < nbf) + GEN p1, lPpow; + GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], tnew, NULL, pk); + TT[i] = (long)p2; + p1 = (GEN)p2[tnew+1]; + /* make Newton sums integral */ + if (lP) { - nbf=anbf; - pr=gcopy(apr); - 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; + lPpow = gpowgs(lP, tnew); + p1 = mulii(p1, lPpow); /* assume pr has degree 1 */ } + 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) - { - fprintferr("prime ideal chosen: %Z\n", pr); - msgtimer("choice of the prime ideal"); - } + /* restart with truncated entries */ + q = shifti(gun, b); + P1 = gdivround(PRK, q); + 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) - { - if (fl) { avma=av; return cgetg(1,t_VEC); } - rep=cgetg(2,t_VEC); rep[1]=lcopy(polmod); - nfcmbf.nfact=1; return gerepileupto(av, rep); + i = lg(CM_L) - 1; + if (i == r && gegal(CM_L, oldCM_L)) + { + CM_L = oldCM_L; + 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); - affir(mulii(absi(dk),gpowgs(p,k)),p2); - p2=shifti(gceil(mplog(p2)),-1); +static GEN +nf_combine_factors(nfcmbf_t *T, GEN polred, GEN p, long a, long klim) +{ + 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); - if (DEBUGLEVEL>=4) - fprintferr("new precision: %ld\n",newprec); + T->fact = hensel_lift_fact(polred,famod,NULL,p,T->pk,a); + if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */ + if (DEBUGLEVEL>3) msgtimer("Hensel lift"); - prh = idealpows(nf,pr,k); m = k; - h = T2_matrix_pow(nf,prh,p,C, &k, newprec); - if (m != k) prh=idealpows(nf,pr,k); + /* FIXME: neither nfcmbf nor LLL_cmbf can handle the non-nf case */ - 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); - msgtimer("computation of H"); + if (l > 1) T->polbase = unifpol(nf, (GEN)res[l], 0); + 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); - lt=(GEN)leading_term(polbase)[1]; - lt=mpinvmod(lt,pk); + m = lg(res); z = cgetg(m, t_VEC); + for (i=1;i 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]); - rep = hensel_lift_fact(polred,fact,NULL,p,pk,k); + if (DEBUGLEVEL>3) msgtimer("square-free"); + 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 */ - hinv= adj(h); - lt=(GEN)leading_term(polbase)[1]; + dk = absi((GEN)nf[3]); + bad = mulii(mulii(dk,(GEN)nf[4]), lt); - 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 mlt = gneg_i(lt), rem; - x_a[1] = polbase[1]; setlgef(x_a, 4); - x_a[3] = un; - p1=cgetg(lg(rep)+1,t_VEC); - polbase = unifpol(nf,polbase,1); - for (m=1,i=1; i3) + 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); - for (i=1; i3) { + msgtimer("choice of a prime ideal"); + fprintferr("Prime ideal chosen: %Z\n", pr); + } - for (i=1; i= 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 (lgef(nfcmbf.pol)>3) - { - nfcmbf.res[++i] = (long) nf_pol_divres(nf,nfcmbf.pol,nfcmbf.lt,NULL); - nfcmbf.nfact = i; + if (DEBUGLEVEL>2) { + msgtimer("bound computation"); + fprintferr(" 1) T_2 bound for %s: %Z\n", fl?"root":"factor", C0); + fprintferr(" 2) Conversion from T_2 --> | |^2 bound : %Z\n", T.ZC); + fprintferr(" 3) Final bound: %Z\n", C); } - tetpil=avma; rep=cgetg(i+1,t_VEC); - for ( ; i>=1; i--) - rep[i]=(long)unifpol(nf,(GEN)nfcmbf.res[i],1); - return gerepile(av,tetpil,rep); -} - -static int -nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint) -{ - int val = 0; /* assume failure */ - 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; + if (fl) + (void)logint(C, p, &pk); +#if 0 /* overkill */ + else + { /* overlift needed for d-1/d-2 tests? */ + GEN pb; long b; /* junk */ + if (cmbf_precs(p, C, T.BS_2, &a, &b, &pk, &pb)) + { /* Rare */ + if (DEBUGLEVEL) err(warner,"nffactor: overlift for d-1/d-2 test"); + C = itor(pk, DEFAULTPREC); + } } +#endif - /* second, try including the current modular factor in the product */ - newf = (GEN)nfcmbf.fact[fxn]; - if (!newf) return val; /* modular factor already used */ - newd = degpol(newf); - if (newd > dlim) return val; /* degree of new factor is too large */ + bestlift_init(0, nf, pr, C, &L); + T.pr = pr; + T.L = &L; - av = avma; - if (newd % hint == 0) - { - GEN p, quot,rem; + /* polred is monic */ + pk = L.pk; + polred = ZpX(gmul(mpinvmod(lt,pk), polbase), pk, L.ZpProj); - newpsf = nf_pol_mul(nf, psf? psf: nfcmbf.lt, newf); - newpsf = nf_pol_lift(nfcmbf.h,nfcmbf.hinv,nfcmbf.den,newpsf); - /* try out the new combination */ - ltop=avma; - quot=nf_pol_divres(nf,nfcmbf.pol,newpsf,&rem); - if (gcmp0(rem)) /* found a factor */ + if (fl) + { /* roots only */ + long x_r[] = { evaltyp(t_POL)|_evallg(4), 0,0,0 }; + rep = rootpadicfast(polred, p, L.k); + x_r[1] = evalsigne(1) | evalvarn(varn(pol)) | evallgef(4); + x_r[3] = un; + for (m=1,i=1; i 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); }