Annotation of OpenXM_contrib/pari-2.2/src/modules/nffactor.c, Revision 1.2
1.2 ! noro 1: /* $Id: nffactor.c,v 1.104 2002/09/07 13:37:18 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* POLYNOMIAL FACTORIZATION IN A NUMBER FIELD */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
23:
1.2 ! noro 24: extern void init_dalloc();
! 25: extern double *dalloc(size_t n);
! 26: extern GEN gmul_mati_smallvec(GEN x, GEN y);
! 27: extern GEN GS_norms(GEN B, long prec);
! 28: extern GEN RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr);
! 29: extern GEN RXQX_red(GEN P, GEN T);
! 30: extern GEN apply_kummer(GEN nf,GEN pol,GEN e,GEN p);
! 31: extern GEN centermod_i(GEN x, GEN p, GEN ps2);
! 32: extern GEN dim1proj(GEN prh);
1.1 noro 33: extern GEN hensel_lift_fact(GEN pol, GEN fact, GEN T, GEN p, GEN pev, long e);
1.2 ! noro 34: extern GEN initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis);
! 35: extern GEN max_modulus(GEN p, double tau);
! 36: extern GEN mulmat_pol(GEN A, GEN x);
! 37: extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
! 38: extern GEN polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N);
1.1 noro 39: extern GEN sort_factor(GEN y, int (*cmp)(GEN,GEN));
1.2 ! noro 40: extern GEN special_pivot(GEN x);
! 41: extern GEN trivfact(void);
! 42: extern GEN vandermondeinverse(GEN L, GEN T, GEN den, GEN prep);
! 43: extern GEN vconcat(GEN A, GEN B);
! 44: extern int cmbf_precs(GEN q, GEN A, GEN B, long *a, long *b, GEN *qa, GEN *qb);
! 45: extern int isrational(GEN x);
! 46: extern GEN LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, pari_timer *T, long *ti_LLL);
! 47: extern void remake_GM(GEN nf, nffp_t *F, long prec);
! 48: #define RXQX_div(x,y,T) RXQX_divrem((x),(y),(T),NULL)
! 49: #define RXQX_rem(x,y,T) RXQX_divrem((x),(y),(T),ONLY_REM)
1.1 noro 50:
51: static GEN nfsqff(GEN nf,GEN pol,long fl);
52:
1.2 ! noro 53: /* for nf_bestlift: reconstruction of algebraic integers known mod \wp^k */
! 54: /* FIXME: assume \wp has degree 1 for now */
! 55: typedef struct {
! 56: long k; /* input known mod \wp^k */
! 57: GEN pk; /* p^k = Norm(\wp^k)*/
! 58: GEN den; /* denom(prk^-1) = p^k */
! 59: GEN prk; /* |.|^2 LLL-reduced basis (b_i) of \wp^k (NOT T2-reduced) */
! 60: GEN iprk; /* den * prk^-1 */
! 61: GEN GSmin; /* min |b_i^*|^2 */
! 62: GEN ZpProj;/* projector to Zp / \wp^k = Z/p^k (\wp unramified, degree 1) */
! 63: GEN tozk;
! 64: GEN topow;
! 65: } nflift_t;
! 66:
! 67: typedef struct /* for use in nfsqff */
! 68: {
! 69: nflift_t *L;
! 70: GEN nf, pol, fact, res, bound, pr, pk, Br, ZC, dn, polbase, BS_2;
! 71: long hint;
! 72: } nfcmbf_t;
! 73:
! 74: GEN
! 75: RXQX_mul(GEN x, GEN y, GEN T)
1.1 noro 76: {
1.2 ! noro 77: return RXQX_red(gmul(x,y), T);
! 78: }
1.1 noro 79:
80: static GEN
81: unifpol0(GEN nf,GEN pol,long flag)
82: {
83: static long n = 0;
84: static GEN vun = NULL;
85: GEN f = (GEN) nf[1];
1.2 ! noro 86: long i = degpol(f);
! 87: gpmem_t av;
1.1 noro 88:
89: if (i != n)
90: {
91: n=i; if (vun) free(vun);
92: vun = (GEN) gpmalloc((n+1)*sizeof(long));
93: vun[0] = evaltyp(t_COL) | evallg(n+1); vun[1]=un;
94: for ( ; i>=2; i--) vun[i]=zero;
95: }
96:
97: av = avma;
98: switch(typ(pol))
99: {
100: case t_INT: case t_FRAC: case t_RFRAC:
1.2 ! noro 101: if (flag) return gcopy(pol);
1.1 noro 102: pol = gmul(pol,vun); break;
103:
104: case t_POL:
1.2 ! noro 105: if (flag && !degpol(pol)) return gcopy(constant_term(pol));
1.1 noro 106: pol = gmodulcp(pol,f); /* fall through */
107: case t_POLMOD:
108: pol = algtobasis(nf,pol);
109: }
110: if (flag) pol = basistoalg(nf, lift(pol));
111: return gerepileupto(av,pol);
112: }
113:
114: /* Let pol be a polynomial with coefficients in Z or nf (vectors or polymods)
115: * return the same polynomial with coefficients expressed:
116: * if flag=0: as vectors (on the integral basis).
117: * if flag=1: as polymods.
118: */
119: GEN
120: unifpol(GEN nf,GEN pol,long flag)
121: {
122: if (typ(pol)==t_POL && varn(pol) < varn(nf[1]))
123: {
124: long i, d = lgef(pol);
125: GEN p1 = pol;
126:
127: pol=cgetg(d,t_POL); pol[1]=p1[1];
128: for (i=2; i<d; i++)
129: pol[i] = (long) unifpol0(nf,(GEN) p1[i],flag);
130:
131: return pol;
132: }
133: return unifpol0(nf,(GEN) pol, flag);
134: }
135:
1.2 ! noro 136: /* factorization of x modulo pr. Assume x already reduced */
! 137: GEN
! 138: FqX_factor(GEN x, GEN T, GEN p)
1.1 noro 139: {
1.2 ! noro 140: GEN rep;
! 141: if (!T)
! 142: {
! 143: rep = factmod0(x, p);
! 144: rep[2] = (long)small_to_vec((GEN)rep[2]);
! 145: }
! 146: else
1.1 noro 147: {
1.2 ! noro 148: rep = factmod9(x, p, T);
! 149: rep = lift_intern(lift_intern(rep));
1.1 noro 150: }
1.2 ! noro 151: return rep;
1.1 noro 152: }
153:
1.2 ! noro 154: GEN
! 155: nffactormod(GEN nf, GEN x, GEN pr)
1.1 noro 156: {
1.2 ! noro 157: long j, l, vx = varn(x), vn;
! 158: gpmem_t av = avma;
! 159: GEN z, rep, xrd, modpr, T, p;
! 160:
! 161: nf = checknf(nf);
! 162: vn = varn((GEN)nf[1]);
! 163: if (typ(x)!=t_POL) err(typeer,"nffactormod");
! 164: if (vx >= vn)
! 165: err(talker,"polynomial variable must have highest priority in nffactormod");
1.1 noro 166:
1.2 ! noro 167: modpr = nf_to_ff_init(nf, &pr, &T, &p);
! 168: xrd = modprX(x, nf, modpr);
! 169: rep = FqX_factor(xrd,T,p);
! 170: z = (GEN)rep[1]; l = lg(z);
! 171: for (j = 1; j < l; j++) z[j] = (long)modprX_lift((GEN)z[j], modpr);
! 172: return gerepilecopy(av, rep);
1.1 noro 173: }
174:
1.2 ! noro 175: /* If p doesn't divide bad and has a divisor of degree 1, return the latter. */
1.1 noro 176: static GEN
1.2 ! noro 177: choose_prime(GEN nf, GEN bad, GEN *p, byteptr *PT)
1.1 noro 178: {
1.2 ! noro 179: GEN q = icopy(gun), r, x = (GEN)nf[1];
! 180: ulong pp = *p? itou(*p): 0;
! 181: byteptr pt = *PT;
! 182: gpmem_t av = avma;
! 183: for (;;)
! 184: {
! 185: NEXT_PRIME_VIADIFF_CHECK(pp, pt);
! 186: if (! smodis(bad,pp)) continue;
! 187: affui(pp, q);
! 188: r = rootmod(x, q); if (lg(r) > 1) break;
! 189: avma = av;
! 190: }
! 191: r = gsub(polx[varn(x)], lift_intern((GEN)r[1]));
! 192: *p = q;
! 193: *PT = pt; return apply_kummer(nf,r,gun,q);
1.1 noro 194: }
195:
196: static GEN
1.2 ! noro 197: QXQ_normalize(GEN P, GEN T)
1.1 noro 198: {
1.2 ! noro 199: GEN t = leading_term(P);
! 200: if (!gcmp1(t))
! 201: {
! 202: if (is_rational_t(typ(t)))
! 203: P = gdiv(P, t);
! 204: else
! 205: P = RXQX_mul(QX_invmod(t,T), P, T);
! 206: }
! 207: return P;
1.1 noro 208: }
209:
1.2 ! noro 210: /* return the roots of pol in nf */
! 211: GEN
! 212: nfroots(GEN nf,GEN pol)
1.1 noro 213: {
1.2 ! noro 214: gpmem_t av = avma;
! 215: int d = degpol(pol);
! 216: GEN A,g, T;
! 217:
! 218: nf = checknf(nf); T = (GEN)nf[1];
! 219: if (typ(pol) != t_POL) err(notpoler,"nfroots");
! 220: if (varn(pol) >= varn(T))
! 221: err(talker,"polynomial variable must have highest priority in nfroots");
! 222: if (d == 0) return cgetg(1,t_VEC);
! 223: if (d == 1)
! 224: {
! 225: A = gneg_i(gdiv((GEN)pol[2],(GEN)pol[3]));
! 226: return gerepilecopy(av, _vec( basistoalg(nf,A) ));
! 227: }
! 228: A = fix_relative_pol(nf,pol,0);
! 229: A = primpart( lift_intern(A) );
! 230: if (DEBUGLEVEL>3) fprintferr("test if polynomial is square-free\n");
! 231: g = nfgcd(A, derivpol(A), T, NULL);
1.1 noro 232:
1.2 ! noro 233: if (degpol(g))
! 234: { /* not squarefree */
! 235: g = QXQ_normalize(g, T);
! 236: A = RXQX_div(A,g,T);
! 237: }
! 238: A = QXQ_normalize(A, T);
! 239: A = primpart(A);
! 240: A = nfsqff(nf,A,1);
! 241: return gerepileupto(av, gen_sort(A, 0, cmp_pol));
1.1 noro 242: }
243:
1.2 ! noro 244: int
! 245: nfissplit(GEN nf, GEN x)
1.1 noro 246: {
1.2 ! noro 247: gpmem_t av = avma;
! 248: long l = lg(nfsqff(checknf(nf), x, 2));
! 249: avma = av; return l != 1;
1.1 noro 250: }
251:
1.2 ! noro 252: /* nf = K[y] / (P). Galois over K ? */
! 253: int
! 254: nfisgalois(GEN nf, GEN P)
1.1 noro 255: {
1.2 ! noro 256: return degpol(P) <= 2 || nfissplit(nf, P);
1.1 noro 257: }
258:
1.2 ! noro 259: /* return a minimal lift of elt modulo id */
1.1 noro 260: static GEN
1.2 ! noro 261: nf_bestlift(GEN elt, GEN bound, nflift_t *T)
1.1 noro 262: {
1.2 ! noro 263: GEN u;
! 264: long i,l = lg(T->prk), t = typ(elt);
! 265: if (t != t_INT)
! 266: {
! 267: if (t == t_POL) elt = mulmat_pol(T->tozk, elt);
! 268: u = gmul(T->iprk,elt);
! 269: for (i=1; i<l; i++) u[i] = (long)diviiround((GEN)u[i], T->pk);
1.1 noro 270: }
1.2 ! noro 271: else
1.1 noro 272: {
1.2 ! noro 273: u = gmul(elt, (GEN)T->iprk[1]);
! 274: for (i=1; i<l; i++) u[i] = (long)diviiround((GEN)u[i], T->pk);
! 275: elt = gscalcol(elt, l-1);
1.1 noro 276: }
1.2 ! noro 277: u = gsub(elt, gmul(T->prk, u));
! 278: if (bound && gcmp(QuickNormL2(u,DEFAULTPREC), bound) > 0) u = NULL;
! 279: return u;
1.1 noro 280: }
281:
282: static GEN
1.2 ! noro 283: nf_bestlift_to_pol(GEN elt, GEN bound, nflift_t *T)
1.1 noro 284: {
1.2 ! noro 285: GEN u = nf_bestlift(elt,bound,T);
! 286: if (u) u = gmul(T->topow, u);
! 287: return u;
1.1 noro 288: }
289:
1.2 ! noro 290: /* return the lift of pol with coefficients of T2-norm <= C (if possible) */
1.1 noro 291: static GEN
1.2 ! noro 292: nf_pol_lift(GEN pol, GEN bound, nfcmbf_t *T)
1.1 noro 293: {
1.2 ! noro 294: long i, d = lgef(pol);
! 295: GEN x = cgetg(d,t_POL);
! 296: nflift_t *L = T->L;
1.1 noro 297:
1.2 ! noro 298: x[1] = pol[1];
! 299: for (i=2; i<d; i++)
1.1 noro 300: {
1.2 ! noro 301: x[i] = (long)nf_bestlift_to_pol((GEN)pol[i], bound, L);
! 302: if (!x[i]) return NULL;
1.1 noro 303: }
1.2 ! noro 304: return x;
1.1 noro 305: }
306:
1.2 ! noro 307: /* return the factorization of x in nf */
! 308: GEN
! 309: nffactor(GEN nf,GEN pol)
1.1 noro 310: {
1.2 ! noro 311: GEN A,g,y,p1,rep,T;
! 312: long l, j, d = degpol(pol);
! 313: gpmem_t av;
! 314: if (DEBUGLEVEL>3) (void)timer2();
! 315:
! 316: nf = checknf(nf); T = (GEN)nf[1];
! 317: if (typ(pol) != t_POL) err(notpoler,"nffactor");
! 318: if (varn(pol) >= varn(T))
! 319: err(talker,"polynomial variable must have highest priority in nffactor");
1.1 noro 320:
1.2 ! noro 321: if (d == 0) return trivfact();
! 322: rep = cgetg(3, t_MAT); av = avma;
! 323: if (d == 1)
! 324: {
! 325: rep[1] = (long)_col( gcopy(pol) );
! 326: rep[2] = (long)_col( gun );
! 327: return rep;
! 328: }
1.1 noro 329:
1.2 ! noro 330: A = fix_relative_pol(nf,pol,0);
! 331: if (degpol(nf[1]) == 1)
! 332: return gerepileupto(av, factpol(simplify(pol), 0));
! 333:
! 334: A = primpart( lift_intern(A) );
! 335: if (DEBUGLEVEL>3) fprintferr("test if polynomial is square-free\n");
! 336: g = nfgcd(A, derivpol(A), T, NULL);
! 337:
! 338: A = QXQ_normalize(A, T);
! 339: A = primpart(A);
! 340:
! 341: if (degpol(g))
! 342: { /* not squarefree */
! 343: gpmem_t av1;
! 344: GEN ex;
! 345: g = QXQ_normalize(g, T);
! 346: A = RXQX_div(A,g, T);
! 347:
! 348: y = nfsqff(nf,A,0); av1 = avma;
! 349: l = lg(y);
! 350: ex=(GEN)gpmalloc(l * sizeof(long));
! 351: for (j=l-1; j>=1; j--)
1.1 noro 352: {
1.2 ! noro 353: GEN fact = lift((GEN)y[j]), quo = g, q;
! 354: long e = 0;
! 355: for(e = 1;; e++)
! 356: {
! 357: q = RXQX_divrem(quo,fact,T, ONLY_DIVIDES);
! 358: if (!q) break;
! 359: quo = q;
! 360: }
! 361: ex[j] = e;
1.1 noro 362: }
1.2 ! noro 363: avma = av1; y = gerepileupto(av, y);
! 364: p1 = cgetg(l, t_COL); for (j=l-1; j>=1; j--) p1[j] = lstoi(ex[j]);
! 365: free(ex);
1.1 noro 366: }
1.2 ! noro 367: else
! 368: {
! 369: y = gerepileupto(av, nfsqff(nf,A,0));
! 370: l = lg(y);
! 371: p1 = cgetg(l, t_COL); for (j=l-1; j>=1; j--) p1[j] = un;
! 372: }
! 373: if (DEBUGLEVEL>3)
! 374: fprintferr("number of factor(s) found: %ld\n", lg(y)-1);
! 375: rep[1] = (long)y;
! 376: rep[2] = (long)p1; return sort_factor(rep, cmp_pol);
! 377: }
! 378:
! 379: /* return a bound for T_2(P), P | polbase in C[X]
! 380: * NB: Mignotte bound: A | S ==>
! 381: * |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S)
! 382: *
! 383: * Apply to sigma(S) for all embeddings sigma, then take the L_2 norm over
! 384: * sigma, then take the sup over i.
! 385: **/
! 386: static GEN
! 387: nf_Mignotte_bound(GEN nf, GEN polbase)
! 388: {
! 389: GEN G = gmael(nf,5,2), lS = leading_term(polbase); /* t_INT */
! 390: GEN p1, C, N2, matGS, binlS, bin;
! 391: long prec, i, j, d = degpol(polbase), n = degpol(nf[1]), r1 = nf_get_r1(nf);
! 392:
! 393: if (typ(lS) == t_COL) lS = (GEN)lS[1];
! 394: binlS = bin = vecbinome(d-1);
! 395: if (!gcmp1(lS)) binlS = gmul(lS, bin);
1.1 noro 396:
1.2 ! noro 397: N2 = cgetg(n+1, t_VEC);
! 398: prec = gprecision(G);
! 399: for (;;)
! 400: {
! 401: nffp_t F;
1.1 noro 402:
1.2 ! noro 403: matGS = cgetg(d+2, t_MAT);
! 404: for (j=0; j<=d; j++) matGS[j+1] = lmul(G, (GEN)polbase[j+2]);
! 405: matGS = gtrans_i(matGS);
! 406: for (j=1; j <= r1; j++) /* N2[j] = || sigma_j(S) ||_2 */
! 407: {
! 408: N2[j] = lsqrt( QuickNormL2((GEN)matGS[j], DEFAULTPREC), DEFAULTPREC );
! 409: if (lg(N2[j]) < DEFAULTPREC) goto PRECPB;
! 410: }
! 411: for ( ; j <= n; j+=2)
! 412: {
! 413: GEN q1 = QuickNormL2((GEN)matGS[j ], DEFAULTPREC);
! 414: GEN q2 = QuickNormL2((GEN)matGS[j+1], DEFAULTPREC);
! 415: p1 = gmul2n(mpadd(q1, q2), -1);
! 416: N2[j] = N2[j+1] = lsqrt( p1, DEFAULTPREC );
! 417: if (lg(N2[j]) < DEFAULTPREC) goto PRECPB;
! 418: }
! 419: if (j > n) break; /* done */
! 420: PRECPB:
! 421: prec = (prec<<1)-2;
! 422: remake_GM(nf, &F, prec); G = F.G;
! 423: if (DEBUGLEVEL>1) err(warnprec, "nf_factor_bound", prec);
! 424: }
1.1 noro 425:
1.2 ! noro 426: /* Take sup over 0 <= i <= d of
! 427: * sum_sigma | binom(d-1, i-1) ||sigma(S)||_2 + binom(d-1,i) lc(S) |^2 */
1.1 noro 428:
1.2 ! noro 429: /* i = 0: n lc(S)^2 */
! 430: C = mulsi(n, sqri(lS));
! 431: /* i = d: sum_sigma ||sigma(S)||_2^2 */
! 432: p1 = gnorml2(N2); if (gcmp(C, p1) < 0) C = p1;
! 433: for (i = 1; i < d; i++)
1.1 noro 434: {
1.2 ! noro 435: GEN s = gzero;
! 436: for (j = 1; j <= n; j++)
1.1 noro 437: {
1.2 ! noro 438: p1 = mpadd( mpmul((GEN)bin[i], (GEN)N2[j]), (GEN)binlS[i+1] );
! 439: s = mpadd(s, gsqr(p1));
1.1 noro 440: }
1.2 ! noro 441: if (gcmp(C, s) < 0) C = s;
1.1 noro 442: }
1.2 ! noro 443: return C;
1.1 noro 444: }
445:
1.2 ! noro 446: /* return a bound for T_2(P), P | polbase
! 447: * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2,
! 448: * where [P]_2 is Bombieri's 2-norm
! 449: * Sum over conjugates
! 450: */
! 451: static GEN
! 452: nf_Beauzamy_bound(GEN nf, GEN polbase)
1.1 noro 453: {
1.2 ! noro 454: GEN lt,C,run,s, G = gmael(nf,5,2), POL, bin;
! 455: long i,prec,precnf, d = degpol(polbase), n = degpol(nf[1]);
1.1 noro 456:
1.2 ! noro 457: precnf = gprecision(G);
! 458: prec = MEDDEFAULTPREC;
! 459: bin = vecbinome(d);
! 460: POL = polbase + 2;
! 461: /* compute [POL]_2 */
! 462: for (;;)
1.1 noro 463: {
1.2 ! noro 464: run= realun(prec);
! 465: s = realzero(prec);
! 466: for (i=0; i<=d; i++)
! 467: {
! 468: GEN p1 = gnorml2(gmul(G, gmul(run, (GEN)POL[i]))); /* T2(POL[i]) */
! 469: if (!signe(p1)) continue;
! 470: if (lg(p1) == 3) break;
! 471: /* s += T2(POL[i]) / binomial(d,i) */
! 472: s = addrr(s, gdiv(p1, (GEN)bin[i+1]));
1.1 noro 473: }
1.2 ! noro 474: if (i > d) break;
1.1 noro 475:
1.2 ! noro 476: prec = (prec<<1)-2;
! 477: if (prec > precnf)
1.1 noro 478: {
1.2 ! noro 479: nffp_t F; remake_GM(nf, &F, prec); G = F.G;
! 480: if (DEBUGLEVEL>1) err(warnprec, "nf_factor_bound", prec);
1.1 noro 481: }
1.2 ! noro 482: }
! 483: lt = (GEN)leading_term(polbase)[1];
! 484: s = gmul(s, mulis(sqri(lt), n));
! 485: C = gpow(stoi(3), dbltor(1.5 + d), DEFAULTPREC); /* 3^{3/2 + d} */
! 486: return gdiv(gmul(C, s), gmulsg(d, mppi(DEFAULTPREC)));
1.1 noro 487: }
488:
489: static GEN
1.2 ! noro 490: nf_factor_bound(GEN nf, GEN polbase)
1.1 noro 491: {
1.2 ! noro 492: gpmem_t av = avma;
! 493: GEN a = nf_Mignotte_bound(nf, polbase);
! 494: GEN b = nf_Beauzamy_bound(nf, polbase);
! 495: if (DEBUGLEVEL>2)
1.1 noro 496: {
1.2 ! noro 497: fprintferr("Mignotte bound: %Z\n",a);
! 498: fprintferr("Beauzamy bound: %Z\n",b);
1.1 noro 499: }
1.2 ! noro 500: return gerepileupto(av, gmin(a, b));
1.1 noro 501: }
502:
1.2 ! noro 503: static long
! 504: ZXY_get_prec(GEN P)
1.1 noro 505: {
1.2 ! noro 506: long i, j, z, prec = 0;
! 507: for (i=2; i<lgef(P); i++)
1.1 noro 508: {
1.2 ! noro 509: GEN p = (GEN)P[i];
! 510: if (typ(p) == t_INT)
! 511: {
! 512: z = lgefint(p);
! 513: if (z > prec) prec = z;
! 514: }
! 515: else
! 516: {
! 517: for (j=2; j<lgef(p); j++)
! 518: {
! 519: z = lgefint(p[j]);
! 520: if (z > prec) prec = z;
! 521: }
! 522: }
! 523: }
! 524: return prec + 1;
1.1 noro 525: }
526:
1.2 ! noro 527: long
! 528: ZM_get_prec(GEN x)
1.1 noro 529: {
1.2 ! noro 530: long i, j, l, k = 2, lx = lg(x);
1.1 noro 531:
1.2 ! noro 532: for (j=1; j<lx; j++)
1.1 noro 533: {
1.2 ! noro 534: GEN c = (GEN)x[j];
! 535: for (i=1; i<lx; i++) { l = lgefint(c[i]); if (l > k) k = l; }
1.1 noro 536: }
1.2 ! noro 537: return k;
1.1 noro 538: }
1.2 ! noro 539: long
! 540: ZX_get_prec(GEN x)
1.1 noro 541: {
1.2 ! noro 542: long j, l, k = 2, lx = lgef(x);
1.1 noro 543:
1.2 ! noro 544: for (j=2; j<lx; j++)
1.1 noro 545: {
1.2 ! noro 546: l = lgefint(x[j]); if (l > k) k = l;
1.1 noro 547: }
1.2 ! noro 548: return k;
1.1 noro 549: }
550:
1.2 ! noro 551: static GEN
! 552: complex_bound(GEN P)
1.1 noro 553: {
1.2 ! noro 554: return gmul(max_modulus(P, 0.01), dbltor(1.0101)); /* exp(0.01) = 1.0100 */
1.1 noro 555: }
556:
1.2 ! noro 557: /* return Bs: if r a root of sigma_i(P), |r| < Bs[i] */
! 558: static GEN
! 559: nf_root_bounds(GEN P, GEN T)
1.1 noro 560: {
1.2 ! noro 561: long lR, i, j, l, prec;
! 562: GEN Ps, R, V, nf;
1.1 noro 563:
1.2 ! noro 564: if (isrational(P)) return complex_bound(P);
1.1 noro 565:
1.2 ! noro 566: T = get_nfpol(T, &nf);
1.1 noro 567:
1.2 ! noro 568: prec = ZXY_get_prec(P);
! 569: l = lgef(P);
! 570: if (nf && nfgetprec(nf) >= prec)
! 571: R = (GEN)nf[6];
! 572: else
! 573: R = roots(T, prec);
! 574: lR = lg(R);
! 575: V = cgetg(lR, t_VEC);
! 576: Ps = cgetg(l, t_POL); /* sigma (P) */
! 577: Ps[1] = P[1];
! 578: for (j=1; j<lg(R); j++)
1.1 noro 579: {
1.2 ! noro 580: GEN r = (GEN)R[j];
! 581: for (i=2; i<l; i++) Ps[i] = (long)poleval((GEN)P[i], r);
! 582: V[j] = (long)complex_bound(Ps);
1.1 noro 583: }
1.2 ! noro 584: return V;
1.1 noro 585: }
586:
1.2 ! noro 587: /* return B such that if x in O_K, K = Z[X]/(T), then the L2-norm of the
! 588: * coordinates of the numerator of x [on the power, resp. integral, basis if T
! 589: * is a polynomial, resp. an nf] is <= B T_2(x)
! 590: * *ptden = multiplicative bound for denom(x) */
1.1 noro 591: static GEN
1.2 ! noro 592: L2_bound(GEN T, GEN tozk, GEN *ptden)
1.1 noro 593: {
1.2 ! noro 594: GEN nf, M, L, prep, den, u;
! 595: long prec;
! 596:
! 597: T = get_nfpol(T, &nf);
! 598: u = NULL; /* gcc -Wall */
1.1 noro 599:
1.2 ! noro 600: prec = ZX_get_prec(T) + ZM_get_prec(tozk);
! 601: den = nf? gun: NULL;
! 602: den = initgaloisborne(T, den, prec, &L, &prep, NULL);
! 603: M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
! 604: if (nf) M = gmul(tozk, M);
! 605: if (gcmp1(den)) den = NULL;
! 606: *ptden = den;
! 607: return QuickNormL2(M, DEFAULTPREC);
1.1 noro 608: }
609:
1.2 ! noro 610: /* || L ||_p^p in dimension n (L may be a scalar) */
1.1 noro 611: static GEN
1.2 ! noro 612: normlp(GEN L, long p, long n)
1.1 noro 613: {
1.2 ! noro 614: long i,l, t = typ(L);
! 615: GEN z;
1.1 noro 616:
1.2 ! noro 617: if (!is_vec_t(t)) return gmulsg(n, gpowgs(L, p));
1.1 noro 618:
1.2 ! noro 619: l = lg(L); z = gzero;
! 620: /* assert(n == l-1); */
! 621: for (i=1; i<l; i++)
! 622: z = gadd(z, gpowgs((GEN)L[i], p));
! 623: return z;
1.1 noro 624: }
625:
1.2 ! noro 626: /* S = S0 + qS1, P = P0 + qP1 (Euclidean div. by q integer). For a true
! 627: * factor (vS, vP), we have:
! 628: * | S vS + P vP |^2 < Btra
! 629: * This implies | S1 vS + P1 vP |^2 < Blow, assuming q > sqrt(Btra).
! 630: * d = dimension of low part (= [nf:Q])
! 631: * n0 = bound for |vS|^2
! 632: * */
! 633: static double
! 634: get_Blow(long n0, long d, GEN q)
! 635: {
! 636: #if 0
! 637: double sqrtn0 = sqrt((double)n0);
! 638: double t = sqrtn0 + sqrt((double)d) * (sqrtn0 + 1);
! 639: #else
! 640: double sqrtd = sqrt((double)d);
! 641: double t, aux;
! 642:
! 643: if (gexpo(q)>30)
! 644: aux = 1.0001;
! 645: else
1.1 noro 646: {
1.2 ! noro 647: aux = gtodouble(q);
! 648: aux = sqrt(1 + 4/(aux*aux));
1.1 noro 649: }
1.2 ! noro 650: t = n0*sqrtd + sqrtd/2 * aux * (sqrtd * (n0+1)); /* assume pr degree 1 */
! 651: #endif
! 652: t = 1. + 0.5 * t;
! 653: return t * t;
! 654: }
! 655:
! 656: typedef struct {
! 657: GEN d;
! 658: GEN dPinvS; /* d P^(-1) S [ integral ] */
! 659: double **PinvSdbl; /* P^(-1) S as double */
! 660: GEN S1, P1; /* S = S0 + S1 q, idem P */
! 661: } trace_data;
! 662:
! 663: /* S1 * u - P1 * round(P^-1 S u). K non-zero coords in u given by ind */
! 664: static GEN
! 665: get_trace(GEN ind, trace_data *T)
! 666: {
! 667: long i, j, l, K = lg(ind)-1;
! 668: GEN z, s, v;
! 669:
! 670: s = (GEN)T->S1[ ind[1] ];
! 671: if (K == 1) return s;
! 672:
! 673: /* compute s = S1 u */
! 674: for (j=2; j<=K; j++) s = gadd(s, (GEN)T->S1[ ind[j] ]);
! 675:
! 676: /* compute v := - round(P^1 S u) */
! 677: l = lg(s);
! 678: v = cgetg(l, t_VECSMALL);
! 679: for (i=1; i<l; i++)
! 680: {
! 681: double r, t = 0.;
! 682: /* quick approximate computation */
! 683: for (j=1; j<=K; j++) t += T->PinvSdbl[ ind[j] ][i];
! 684: r = floor(t + 0.5);
! 685: if (fabs(t + 0.5 - r) < 0.01)
! 686: { /* dubious, compute exactly */
! 687: z = gzero;
! 688: for (j=1; j<=K; j++) z = addii(z, ((GEN**)T->dPinvS)[ ind[j] ][i]);
! 689: v[i] = - itos( diviiround(z, T->d) );
! 690: }
! 691: else
! 692: v[i] = - (long)r;
1.1 noro 693: }
1.2 ! noro 694: return gadd(s, gmul_mati_smallvec(T->P1, v));
! 695: }
1.1 noro 696:
1.2 ! noro 697: static trace_data *
! 698: init_trace(trace_data *T, GEN S, nflift_t *L, GEN q)
! 699: {
! 700: long e = gexpo((GEN)S), i,j, l,h;
! 701: GEN qgood, S1, invd;
! 702:
! 703: if (e < 0) return NULL; /* S = 0 */
1.1 noro 704:
1.2 ! noro 705: qgood = shifti(gun, e - 32); /* single precision check */
! 706: if (cmpii(qgood, q) > 0) q = qgood;
1.1 noro 707:
1.2 ! noro 708: S1 = gdivround(S, q);
! 709: if (gcmp0(S1)) return NULL;
1.1 noro 710:
1.2 ! noro 711: invd = ginv(itor(L->den, DEFAULTPREC));
1.1 noro 712:
1.2 ! noro 713: T->dPinvS = gmul(L->iprk, S);
! 714: l = lg(S);
! 715: h = lg(T->dPinvS[1]);
! 716: T->PinvSdbl = (double**)cgetg(l, t_MAT);
! 717: init_dalloc();
! 718: for (j = 1; j < l; j++)
1.1 noro 719: {
1.2 ! noro 720: double *t = dalloc(h * sizeof(double));
! 721: GEN c = (GEN)T->dPinvS[j];
! 722: T->PinvSdbl[j] = t;
! 723: for (i=1; i < h; i++) t[i] = rtodbl(gmul(invd, (GEN)c[i]));
! 724: }
! 725:
! 726: T->d = L->den;
! 727: T->P1 = gdivround(L->prk, q);
! 728: T->S1 = S1; return T;
! 729: }
! 730:
! 731: static void
! 732: update_trace(trace_data *T, long k, long i)
! 733: {
! 734: if (!T) return;
! 735: T->S1[k] = T->S1[i];
! 736: T->dPinvS[k] = T->dPinvS[i];
! 737: T->PinvSdbl[k] = T->PinvSdbl[i];
! 738: }
! 739:
! 740: /* Naive recombination of modular factors: combine up to maxK modular
! 741: * factors, degree <= klim and divisible by hint
! 742: *
! 743: * target = polynomial we want to factor
! 744: * famod = array of modular factors. Product should be congruent to
! 745: * target/lc(target) modulo p^a
! 746: * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
! 747: static GEN
! 748: nfcmbf(nfcmbf_t *T, GEN p, long a, long maxK, long klim)
! 749: {
! 750: long Sbound;
! 751: GEN pol = T->pol, nf = T->nf, famod = T->fact;
! 752: GEN bound = T->bound;
! 753: GEN nfpol = (GEN)nf[1];
! 754: long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol);
! 755: GEN lc, lcpol;
! 756: GEN pk = gpowgs(p,a), pas2 = shifti(pk,-1);
! 757:
! 758: GEN trace1 = cgetg(lfamod+1, t_MAT);
! 759: GEN trace2 = cgetg(lfamod+1, t_MAT);
! 760: GEN ind = cgetg(lfamod+1, t_VECSMALL);
! 761: GEN degpol = cgetg(lfamod+1, t_VECSMALL);
! 762: GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
! 763: GEN listmod = cgetg(lfamod+1, t_COL);
! 764: GEN fa = cgetg(lfamod+1, t_COL);
! 765: GEN res = cgetg(3, t_VEC);
! 766: GEN q = ceil_safe(mpsqrt(T->BS_2));
! 767: const double Blow = get_Blow(lfamod, dnf, q);
! 768: trace_data _T1, _T2, *T1, *T2;
! 769:
! 770: if (maxK < 0) maxK = lfamod-1;
! 771:
! 772: lc = absi(leading_term(pol));
! 773: if (gcmp1(lc)) lc = NULL;
! 774: lcpol = lc? gmul(lc,pol): pol;
1.1 noro 775:
776: {
1.2 ! noro 777: GEN t1,t2, lc2 = lc? sqri(lc): NULL;
1.1 noro 778:
1.2 ! noro 779: for (i=1; i <= lfamod; i++)
! 780: {
! 781: GEN P = (GEN)famod[i];
! 782: long d = degpol(P);
1.1 noro 783:
1.2 ! noro 784: degpol[i] = d; P += 2;
! 785: t1 = (GEN)P[d-1];/* = - S_1 */
! 786: t2 = sqri(t1);
! 787: if (d > 1) t2 = subii(t2, shifti((GEN)P[d-2],1));
! 788: t2 = modii(t2, pk); /* = S_2 Newton sum */
! 789: if (lc)
! 790: {
! 791: t1 = modii(mulii(lc, t1), pk);
! 792: t2 = modii(mulii(lc2,t2), pk);
! 793: }
! 794: trace1[i] = (long)nf_bestlift(t1, NULL, T->L);
! 795: trace2[i] = (long)nf_bestlift(t2, NULL, T->L);
! 796: }
1.1 noro 797:
1.2 ! noro 798: T1 = init_trace(&_T1, trace1, T->L, q);
! 799: T2 = init_trace(&_T2, trace2, T->L, q);
! 800: }
! 801: degsofar[0] = 0; /* sentinel */
1.1 noro 802:
1.2 ! noro 803: /* ind runs through strictly increasing sequences of length K,
! 804: * 1 <= ind[i] <= lfamod */
! 805: nextK:
! 806: if (K > maxK || 2*K > lfamod) goto END;
! 807: if (DEBUGLEVEL > 3)
! 808: fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K));
! 809: setlg(ind, K+1); ind[1] = 1;
! 810: Sbound = ((K+1)>>1);
! 811: i = 1; curdeg = degpol[ind[1]];
! 812: for(;;)
! 813: { /* try all combinations of K factors */
! 814: for (j = i; j < K; j++)
1.1 noro 815: {
1.2 ! noro 816: degsofar[j] = curdeg;
! 817: ind[j+1] = ind[j]+1; curdeg += degpol[ind[j+1]];
! 818: }
! 819: if (curdeg <= klim && curdeg % T->hint == 0) /* trial divide */
! 820: {
! 821: GEN t, y, q, list;
! 822: gpmem_t av;
! 823:
! 824: av = avma;
! 825: /* d - 1 test */
! 826: if (T1)
! 827: {
! 828: t = get_trace(ind, T1);
! 829: if (rtodbl(QuickNormL2(t,DEFAULTPREC)) > Blow)
! 830: {
! 831: if (DEBUGLEVEL>6) fprintferr(".");
! 832: avma = av; goto NEXT;
! 833: }
! 834: }
! 835: /* d - 2 test */
! 836: if (T2)
! 837: {
! 838: t = get_trace(ind, T2);
! 839: if (rtodbl(QuickNormL2(t,DEFAULTPREC)) > Blow)
! 840: {
! 841: if (DEBUGLEVEL>3) fprintferr("|");
! 842: avma = av; goto NEXT;
! 843: }
! 844: }
! 845: avma = av;
! 846: y = lc; /* full computation */
! 847: for (i=1; i<=K; i++)
! 848: {
! 849: GEN q = (GEN)famod[ind[i]];
! 850: if (y) q = gmul(y, q);
! 851: y = centermod_i(q, pk, pas2);
! 852: }
! 853: y = nf_pol_lift(y, bound, T);
! 854: if (!y)
! 855: {
! 856: if (DEBUGLEVEL>3) fprintferr("@");
! 857: avma = av; goto NEXT;
! 858: }
! 859: /* try out the new combination: y is the candidate factor */
! 860: q = RXQX_divrem(lcpol,y, nfpol, ONLY_DIVIDES);
! 861: if (!q)
! 862: {
! 863: if (DEBUGLEVEL>3) fprintferr("*");
! 864: avma = av; goto NEXT;
! 865: }
1.1 noro 866:
1.2 ! noro 867: /* found a factor */
! 868: list = cgetg(K+1, t_VEC);
! 869: listmod[cnt] = (long)list;
! 870: for (i=1; i<=K; i++) list[i] = famod[ind[i]];
! 871:
! 872: y = primpart(y);
! 873: fa[cnt++] = (long)QXQ_normalize(y, nfpol);
! 874: /* fix up pol */
! 875: pol = q;
! 876: if (lc) pol = primpart(pol);
! 877: for (i=j=k=1; i <= lfamod; i++)
! 878: { /* remove used factors */
! 879: if (j <= K && i == ind[j]) j++;
! 880: else
! 881: {
! 882: famod[k] = famod[i];
! 883: update_trace(T1, k, i);
! 884: update_trace(T2, k, i);
! 885: degpol[k] = degpol[i]; k++;
! 886: }
! 887: }
! 888: lfamod -= K;
! 889: if (lfamod < 2*K) goto END;
! 890: i = 1; curdeg = degpol[ind[1]];
! 891: if (lc) lc = absi(leading_term(pol));
! 892: lcpol = lc? gmul(lc,pol): pol;
! 893: if (DEBUGLEVEL > 2)
1.1 noro 894: {
1.2 ! noro 895: fprintferr("\n"); msgtimer("to find factor %Z",y);
! 896: fprintferr("remaining modular factor(s): %ld\n", lfamod);
! 897: }
! 898: continue;
! 899: }
! 900:
! 901: NEXT:
! 902: for (i = K+1;;)
! 903: {
! 904: if (--i == 0) { K++; goto nextK; }
! 905: if (++ind[i] <= lfamod - K + i)
! 906: {
! 907: curdeg = degsofar[i-1] + degpol[ind[i]];
! 908: if (curdeg <= klim) break;
1.1 noro 909: }
910: }
911: }
1.2 ! noro 912: END:
! 913: if (degpol(pol) > 0)
! 914: { /* leftover factor */
! 915: if (signe(leading_term(pol)) < 0) pol = gneg_i(pol);
! 916:
! 917: if (lc && lfamod < 2*K) pol = QXQ_normalize(primpart(pol), nfpol);
! 918: setlg(famod, lfamod+1);
! 919: listmod[cnt] = (long)dummycopy(famod);
! 920: fa[cnt++] = (long)pol;
! 921: }
! 922: if (DEBUGLEVEL>6) fprintferr("\n");
! 923: setlg(listmod, cnt); setlg(fa, cnt);
! 924: res[1] = (long)fa;
! 925: res[2] = (long)listmod; return res;
1.1 noro 926: }
927:
1.2 ! noro 928: static GEN
! 929: nf_check_factors(nfcmbf_t *T, GEN P, GEN M_L, GEN famod, GEN pk)
! 930: {
! 931: GEN nf = T->nf, bound = T->bound;
! 932: GEN nfT = (GEN)nf[1];
! 933: long i, j, r, n0;
! 934: GEN pol = P, lcpol, lc, list, piv, y, pas2;
! 935:
! 936: piv = special_pivot(M_L);
! 937: if (!piv) return NULL;
! 938: if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv);
! 939:
! 940: pas2 = shifti(pk,-1);
! 941: r = lg(piv)-1;
! 942: n0 = lg(piv[1])-1;
! 943: list = cgetg(r+1, t_COL);
! 944: lc = absi(leading_term(pol));
! 945: if (is_pm1(lc)) lc = NULL;
! 946: lcpol = lc? gmul(lc, pol): pol;
! 947: for (i=1; i<r; i++)
! 948: {
! 949: GEN c = (GEN)piv[i];
! 950: if (DEBUGLEVEL) fprintferr("nf_LLL_cmbf: checking factor %ld\n",i);
! 951:
! 952: y = lc;
! 953: for (j=1; j<=n0; j++)
! 954: if (signe(c[j]))
! 955: {
! 956: GEN q = (GEN)famod[j];
! 957: if (y) q = gmul(y, q);
! 958: y = centermod_i(q, pk, pas2);
! 959: }
! 960: y = nf_pol_lift(y, bound, T);
! 961: if (!y) return NULL;
! 962:
! 963: /* y is the candidate factor */
! 964: pol = RXQX_divrem(lcpol,y,nfT, ONLY_DIVIDES);
! 965: if (!pol) return NULL;
1.1 noro 966:
1.2 ! noro 967: y = primpart(y);
! 968: if (lc)
! 969: {
! 970: pol = primpart(pol);
! 971: lc = absi(leading_term(pol));
! 972: }
! 973: lcpol = lc? gmul(lc, pol): pol;
! 974: list[i] = (long)QXQ_normalize(y, nfT);
1.1 noro 975: }
1.2 ! noro 976: y = primpart(pol);
! 977: list[i] = (long)QXQ_normalize(y, nfT); return list;
1.1 noro 978: }
979:
980: static GEN
1.2 ! noro 981: nf_to_Zp(GEN x, GEN pk, GEN pas2, GEN ffproj)
1.1 noro 982: {
1.2 ! noro 983: return centermodii(gmul(ffproj, x), pk, pas2);
! 984: }
1.1 noro 985:
1.2 ! noro 986: /* assume lc(pol) != 0 mod prh. Reduce P to Zp[X] mod p^a */
! 987: static GEN
! 988: ZpX(GEN P, GEN pk, GEN ffproj)
! 989: {
! 990: long i, l;
! 991: GEN z, pas2 = shifti(pk,-1);
1.1 noro 992:
1.2 ! noro 993: if (typ(P)!=t_POL) return nf_to_Zp(P,pk,pas2,ffproj);
! 994: l = lgef(P);
! 995: z = cgetg(l,t_POL); z[1] = P[1];
! 996: for (i=2; i<l; i++) z[i] = (long)nf_to_Zp((GEN)P[i],pk,pas2,ffproj);
! 997: return normalizepol(z);
! 998: }
! 999:
! 1000: /* We want to be able to reconstruct x, |x|^2 < C, from x mod pr^k
! 1001: * We want B := min B_i >= C. Let alpha the LLL parameter,
! 1002: * y = 1/(alpha-1/4) > 4/3, the theoretical guaranteed bound follows from
! 1003: * (Npr)^(2k/n) = det(L)^(2/n) <= 1/n sum B_i
! 1004: * <= 1/n B sum y^i
! 1005: * <= 1/n B y^(n-1) / (y-1)
! 1006: *
! 1007: * Hence log(B/n) >= 2k/n log(Npr) - (n-1)log(y) + log(y-1)
! 1008: * >= log(C/n), provided
! 1009: * k >= [ log(C/n) + (n-1)log(y) - log(y-1) ] * n/(2log(Npr))
! 1010: */
! 1011: static double
! 1012: bestlift_bound(GEN C, long n, double alpha, GEN Npr)
! 1013: {
! 1014: const double y = 1 / (alpha - 0.25); /* = 2 if alpha = 3/4 */
! 1015: double t;
! 1016: if (typ(C) != t_REAL) C = gmul(C, realun(DEFAULTPREC));
! 1017: setlg(C, DEFAULTPREC);
! 1018: t = rtodbl(mplog(divrs(C,n))) + (n-1) * log(y) - log(y-1);
! 1019: return ceil((t * n) / (2.* log(gtodouble(Npr))));
! 1020: }
! 1021:
! 1022: static void
! 1023: bestlift_init(long a, GEN nf, GEN pr, GEN C, nflift_t *T)
! 1024: {
! 1025: const int D = 4;
! 1026: const double alpha = ((double)D-1) / D; /* LLL parameter */
! 1027: gpmem_t av = avma;
! 1028: GEN prk, PRK, B, GSmin, pk;
! 1029:
! 1030: if (!a) a = (long)bestlift_bound(C, degpol(nf[1]), alpha, idealnorm(nf,pr));
! 1031:
! 1032: for (;; avma = av, a<<=1)
! 1033: {
! 1034: if (DEBUGLEVEL>2) fprintferr("exponent: %ld\n",a);
! 1035: prk = idealpows(nf, pr, a);
! 1036: pk = gcoeff(prk,1,1);
! 1037: PRK = hnfmodid(prk, pk);
! 1038:
! 1039: PRK = lllint_i(PRK, D, 0, NULL, NULL, &B);
! 1040: GSmin = vecmin(GS_norms(B, DEFAULTPREC));
! 1041: if (gcmp(GSmin, C) >= 0) break;
! 1042: }
! 1043: if (DEBUGLEVEL>2) fprintferr("for this exponent, GSmin = %Z\n",GSmin);
! 1044: T->k = a;
! 1045: T->pk = T->den = pk;
! 1046: T->prk = PRK;
! 1047: T->iprk = ZM_inv(PRK, pk);
! 1048: T->GSmin= GSmin;
! 1049: T->ZpProj = dim1proj(prk); /* nf --> Zp */
! 1050: }
! 1051:
! 1052: /* assume pr has degree 1 */
! 1053: static GEN
! 1054: nf_LLL_cmbf(nfcmbf_t *T, GEN p, long k, long rec)
! 1055: {
! 1056: nflift_t *L = T->L;
! 1057: GEN pk = L->pk, PRK = L->prk, PRKinv = L->iprk, GSmin = L->GSmin;
! 1058:
! 1059: GEN famod = T->fact, nf = T->nf, ZC = T->ZC, Br = T->Br;
! 1060: GEN Pbase = T->polbase, P = T->pol, dn = T->dn;
! 1061: GEN nfT = (GEN)nf[1];
! 1062: GEN Btra;
! 1063: long dnf = degpol(nfT), dP = degpol(P);
! 1064:
! 1065: double BitPerFactor = 0.5; /* nb bits / modular factor */
! 1066: long i, C, tmax, n0;
! 1067: GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
! 1068: double Blow;
! 1069: gpmem_t av, av2, lim;
! 1070: long ti_LLL = 0, ti_CF = 0;
! 1071: pari_timer ti, ti2, TI;
! 1072:
! 1073: if (DEBUGLEVEL>2) (void)TIMER(&ti);
! 1074:
! 1075: lP = absi(leading_term(P));
! 1076: if (is_pm1(lP)) lP = NULL;
! 1077:
! 1078: n0 = lg(famod) - 1;
! 1079: /* Lattice: (S PRK), small vector (vS vP). To find k bound for the image,
! 1080: * write S = S1 q + S0, P = P1 q + P0
! 1081: * |S1 vS + P1 vP|^2 <= Blow for all (vS,vP) assoc. to true factors */
! 1082: Btra = mulrr(ZC, mulsr(dP*dP, normlp(Br, 2, dnf)));
! 1083: Blow = get_Blow(n0, dnf, gceil(Btra));
! 1084: C = (long)ceil(sqrt(Blow/n0)) + 1; /* C^2 n0 ~ Blow */
! 1085: Bnorm = dbltor( n0 * C * C + Blow );
! 1086: ZERO = zeromat(n0, dnf);
! 1087:
! 1088: av = avma; lim = stack_lim(av, 1);
! 1089: TT = cgetg(n0+1, t_VEC);
! 1090: Tra = cgetg(n0+1, t_MAT);
! 1091: for (i=1; i<=n0; i++) TT[i] = 0;
! 1092: CM_L = gscalsmat(C, n0);
! 1093: /* tmax = current number of traces used (and computed so far) */
! 1094: for(tmax = 0;; tmax++)
! 1095: {
! 1096: long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1;
! 1097: GEN oldCM_L, M_L, q, S1, P1, VV;
! 1098: int first = 1;
! 1099:
! 1100: /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */
! 1101: Btra = mulrr(ZC, mulsr(dP*dP, normlp(Br, 2*tnew, dnf)));
! 1102: bmin = logint(ceil_safe(mpsqrt(Btra)), gdeux, NULL);
! 1103: if (DEBUGLEVEL>2)
! 1104: fprintferr("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
! 1105: r, tmax, bmin);
1.1 noro 1106:
1.2 ! noro 1107: /* compute Newton sums (possibly relifting first) */
! 1108: if (gcmp(GSmin, Btra) < 0)
! 1109: {
! 1110: nflift_t L1;
! 1111: GEN polred;
1.1 noro 1112:
1.2 ! noro 1113: bestlift_init(k<<1, nf, T->pr, Btra, &L1);
! 1114: k = L1.k;
! 1115: pk = L1.pk;
! 1116: PRK = L1.prk;
! 1117: PRKinv = L1.iprk;
! 1118: GSmin = L1.GSmin;
! 1119:
! 1120: polred = ZpX(Pbase, pk, L1.ZpProj);
! 1121: famod = hensel_lift_fact(polred,famod,NULL,p,pk,k);
! 1122: for (i=1; i<=n0; i++) TT[i] = 0;
! 1123: }
! 1124: for (i=1; i<=n0; i++)
1.1 noro 1125: {
1.2 ! noro 1126: GEN p1, lPpow;
! 1127: GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], tnew, NULL, pk);
! 1128: TT[i] = (long)p2;
! 1129: p1 = (GEN)p2[tnew+1];
! 1130: /* make Newton sums integral */
! 1131: if (lP)
! 1132: {
! 1133: lPpow = gpowgs(lP, tnew);
! 1134: p1 = mulii(p1, lPpow); /* assume pr has degree 1 */
! 1135: }
! 1136: if (dn) p1 = mulii(p1,dn);
! 1137: if (dn || lP) p1 = modii(p1, pk);
! 1138: Tra[i] = (long)nf_bestlift(p1, NULL, L); /* S_tnew(famod) */
! 1139: }
1.1 noro 1140:
1.2 ! noro 1141: /* compute truncation parameter */
! 1142: if (DEBUGLEVEL>2) { TIMERstart(&ti2); TIMERstart(&TI); }
! 1143: oldCM_L = CM_L;
! 1144: av2 = avma;
! 1145: b = delta = 0; /* -Wall */
! 1146: AGAIN:
! 1147: M_L = gdivexact(CM_L, stoi(C));
! 1148: T2 = gmul(Tra, M_L);
! 1149: VV = gdivround(gmul(PRKinv, T2), pk);
! 1150: T2 = gsub(T2, gmul(PRK, VV));
! 1151:
! 1152: if (first)
! 1153: { /* initialize lattice, using few p-adic digits for traces */
! 1154: a = gexpo(T2);
! 1155: bgood = (long)(a - max(32, BitPerFactor * r));
! 1156: b = max(bmin, bgood);
! 1157: delta = a - b;
1.1 noro 1158: }
1.2 ! noro 1159: else
! 1160: { /* add more p-adic digits and continue reduction */
! 1161: long b0 = gexpo(T2);
! 1162: if (b0 < b) b = b0;
! 1163: b = max(b-delta, bmin);
! 1164: if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
! 1165: }
! 1166:
! 1167: /* restart with truncated entries */
! 1168: q = shifti(gun, b);
! 1169: P1 = gdivround(PRK, q);
! 1170: S1 = gdivround(Tra, q);
! 1171: T2 = gsub(gmul(S1, M_L), gmul(P1, VV));
! 1172: m = vconcat( CM_L, T2 );
! 1173: if (first)
! 1174: {
! 1175: first = 0;
! 1176: m = concatsp( m, vconcat(ZERO, P1) );
! 1177: /* [ C M_L 0 ]
! 1178: * m = [ ] square matrix
! 1179: * [ T2' PRK ] T2' = Tra * M_L truncated
! 1180: */
! 1181: }
! 1182: CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti, &ti_LLL);
! 1183: if (DEBUGLEVEL>2)
! 1184: fprintferr("LLL_cmbf: b =%4ld; r =%3ld -->%3ld, time = %ld\n",
! 1185: b, lg(m)-1, CM_L? lg(CM_L)-1: 1, TIMER(&TI));
! 1186: if (!CM_L) { list = _col(QXQ_normalize(P,nfT)); break; }
! 1187: if (b > bmin)
! 1188: {
! 1189: CM_L = gerepilecopy(av2, CM_L);
! 1190: goto AGAIN;
! 1191: }
! 1192: if (DEBUGLEVEL>2) msgTIMER(&ti2, "for this trace");
! 1193:
! 1194: i = lg(CM_L) - 1;
! 1195: if (i == r && gegal(CM_L, oldCM_L))
1.1 noro 1196: {
1.2 ! noro 1197: CM_L = oldCM_L;
! 1198: avma = av2; continue;
1.1 noro 1199: }
1200:
1.2 ! noro 1201: if (i <= r && i*rec < n0)
! 1202: {
! 1203: if (DEBUGLEVEL>2) (void)TIMER(&ti);
! 1204: list = nf_check_factors(T, P, M_L, famod, pk);
! 1205: if (DEBUGLEVEL>2) ti_CF += TIMER(&ti);
! 1206: if (list) break;
! 1207: CM_L = gerepilecopy(av2, CM_L);
! 1208: }
! 1209: if (low_stack(lim, stack_lim(av,1)))
1.1 noro 1210: {
1.2 ! noro 1211: if(DEBUGMEM>1) err(warnmem,"nf_LLL_cmbf");
! 1212: gerepileall(av, 8, &CM_L,&TT,&Tra,&famod,&pk,&GSmin,&PRK,&PRKinv);
1.1 noro 1213: }
1214: }
1.2 ! noro 1215: if (DEBUGLEVEL>2)
! 1216: fprintferr("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
! 1217: return list;
1.1 noro 1218: }
1219:
1220: static GEN
1.2 ! noro 1221: nf_combine_factors(nfcmbf_t *T, GEN polred, GEN p, long a, long klim)
1.1 noro 1222: {
1.2 ! noro 1223: GEN z, res, L, listmod, famod = T->fact, nf = T->nf;
! 1224: long i, m, l, maxK = 3, nft = lg(famod)-1;
1.1 noro 1225:
1.2 ! noro 1226: T->fact = hensel_lift_fact(polred,famod,NULL,p,T->pk,a);
! 1227: if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */
! 1228: if (DEBUGLEVEL>3) msgtimer("Hensel lift");
1.1 noro 1229:
1.2 ! noro 1230: /* FIXME: neither nfcmbf nor LLL_cmbf can handle the non-nf case */
1.1 noro 1231:
1.2 ! noro 1232: T->res = cgetg(nft+1,t_VEC);
! 1233: L = nfcmbf(T, p, a, maxK, klim);
1.1 noro 1234:
1.2 ! noro 1235: res = (GEN)L[1];
! 1236: listmod = (GEN)L[2]; l = lg(listmod)-1;
! 1237: famod = (GEN)listmod[l];
! 1238: if (maxK >= 0 && lg(famod)-1 > 2*maxK)
1.1 noro 1239: {
1.2 ! noro 1240: if (l > 1) T->polbase = unifpol(nf, (GEN)res[l], 0);
! 1241: L = nf_LLL_cmbf(T, p, a, maxK);
! 1242: /* remove last elt, possibly unfactored. Add all new ones. */
! 1243: setlg(res, l); res = concatsp(res, L);
1.1 noro 1244: }
1.2 ! noro 1245: if (DEBUGLEVEL>3) msgtimer("computation of the factors");
1.1 noro 1246:
1.2 ! noro 1247: m = lg(res); z = cgetg(m, t_VEC);
! 1248: for (i=1;i<m; i++) z[i] = (long)unifpol(nf,(GEN)res[i],1);
! 1249: return z;
! 1250: }
1.1 noro 1251:
1.2 ! noro 1252: /* return the factorization of the square-free polynomial x.
! 1253: The coeff of x are in Z_nf and its leading term is a rational integer.
! 1254: deg(x) > 1, deg(nfpol) > 1
! 1255: If fl = 1, return only the roots of x in nf
! 1256: If fl = 2, as fl=1 if pol splits, [] otherwise */
! 1257: static GEN
! 1258: nfsqff(GEN nf, GEN pol, long fl)
! 1259: {
! 1260: long i, m, n, nbf, ct, dpol = degpol(pol);
! 1261: gpmem_t av = avma;
! 1262: GEN pr, C0, C, dk, bad, polbase, pk;
! 1263: GEN N2, rep, p, ap, polmod, polred, lt, nfpol;
! 1264: byteptr pt=diffptr;
! 1265: nfcmbf_t T;
! 1266: nflift_t L;
1.1 noro 1267:
1.2 ! noro 1268: if (DEBUGLEVEL>3) msgtimer("square-free");
! 1269: nfpol = (GEN)nf[1]; n = degpol(nfpol);
! 1270: polbase = unifpol(nf,pol,0);
! 1271: polmod = unifpol(nf,pol,1);
! 1272: /* heuristic */
! 1273: #if 1
! 1274: if (dpol*4 < n)
1.1 noro 1275: {
1.2 ! noro 1276: if (DEBUGLEVEL>2) fprintferr("Using Trager's method\n");
! 1277: return gerepilecopy(av, (GEN)polfnf(polmod, nfpol)[1]);
1.1 noro 1278: }
1.2 ! noro 1279: #endif
1.1 noro 1280:
1.2 ! noro 1281: pol = simplify_i(lift(polmod));
! 1282: lt = (GEN)leading_term(polbase)[1]; /* t_INT */
1.1 noro 1283:
1.2 ! noro 1284: dk = absi((GEN)nf[3]);
! 1285: bad = mulii(mulii(dk,(GEN)nf[4]), lt);
1.1 noro 1286:
1.2 ! noro 1287: p = polred = pr = NULL; /* gcc -Wall */
! 1288: nbf = 0; ap = NULL;
! 1289: for (ct = 5;;)
! 1290: {
! 1291: GEN apr = choose_prime(nf, bad, &ap, &pt);
! 1292: GEN modpr = zkmodprinit(nf, apr);
! 1293: long anbf;
1.1 noro 1294:
1.2 ! noro 1295: polred = modprX(polbase, nf, modpr);
! 1296: if (!FpX_is_squarefree(polred,ap)) continue;
1.1 noro 1297:
1.2 ! noro 1298: anbf = fl? FpX_nbroots(polred,ap): FpX_nbfact(polred,ap);
! 1299: if (!nbf || anbf < nbf)
1.1 noro 1300: {
1.2 ! noro 1301: nbf = anbf; pr = apr; p = ap;
! 1302: if (DEBUGLEVEL>3)
! 1303: fprintferr("%3ld %s at prime ideal above %Z\n",
! 1304: nbf, fl?"roots": "factors", p);
! 1305: if (fl == 2 && nbf < dpol) return cgetg(1,t_VEC);
! 1306: if (nbf <= 1)
1.1 noro 1307: {
1.2 ! noro 1308: if (!fl) /* irreducible */
! 1309: return gerepilecopy(av, _vec(QXQ_normalize(polmod, nfpol)));
! 1310: if (!nbf) return cgetg(1,t_VEC); /* no root */
1.1 noro 1311: }
1312: }
1.2 ! noro 1313: if (--ct <= 0) break;
1.1 noro 1314: }
1.2 ! noro 1315: if (DEBUGLEVEL>3) {
! 1316: msgtimer("choice of a prime ideal");
! 1317: fprintferr("Prime ideal chosen: %Z\n", pr);
! 1318: }
! 1319:
! 1320: L.tozk = (GEN)nf[8];
! 1321: L.topow= (GEN)nf[7];
! 1322: T.ZC = L2_bound(nf, L.tozk, &(T.dn));
! 1323: T.Br = gmul(lt, nf_root_bounds(pol, nf));
! 1324:
! 1325: if (fl) C0 = normlp(T.Br, 2, n);
! 1326: else C0 = nf_factor_bound(nf, polbase); /* bound for T_2(Q_i), Q | P */
! 1327: C = mulrr(T.ZC, C0); /* bound for |Q_i|^2 in Z^n on chosen Z-basis */
! 1328: T.bound = C;
! 1329:
! 1330: N2 = mulsr(dpol*dpol, normlp(T.Br, 4, n)); /* bound for T_2(lt * S_2) */
! 1331: T.BS_2 = mulrr(T.ZC, N2); /* bound for |S_2|^2 on chosen Z-basis */
! 1332:
! 1333: if (DEBUGLEVEL>2) {
! 1334: msgtimer("bound computation");
! 1335: fprintferr(" 1) T_2 bound for %s: %Z\n", fl?"root":"factor", C0);
! 1336: fprintferr(" 2) Conversion from T_2 --> | |^2 bound : %Z\n", T.ZC);
! 1337: fprintferr(" 3) Final bound: %Z\n", C);
1.1 noro 1338: }
1339:
1.2 ! noro 1340: if (fl)
! 1341: (void)logint(C, p, &pk);
! 1342: #if 0 /* overkill */
! 1343: else
! 1344: { /* overlift needed for d-1/d-2 tests? */
! 1345: GEN pb; long b; /* junk */
! 1346: if (cmbf_precs(p, C, T.BS_2, &a, &b, &pk, &pb))
! 1347: { /* Rare */
! 1348: if (DEBUGLEVEL) err(warner,"nffactor: overlift for d-1/d-2 test");
! 1349: C = itor(pk, DEFAULTPREC);
! 1350: }
1.1 noro 1351: }
1.2 ! noro 1352: #endif
1.1 noro 1353:
1.2 ! noro 1354: bestlift_init(0, nf, pr, C, &L);
! 1355: T.pr = pr;
! 1356: T.L = &L;
! 1357:
! 1358: /* polred is monic */
! 1359: pk = L.pk;
! 1360: polred = ZpX(gmul(mpinvmod(lt,pk), polbase), pk, L.ZpProj);
1.1 noro 1361:
1362: if (fl)
1.2 ! noro 1363: { /* roots only */
! 1364: long x_r[] = { evaltyp(t_POL)|_evallg(4), 0,0,0 };
! 1365: rep = rootpadicfast(polred, p, L.k);
! 1366: x_r[1] = evalsigne(1) | evalvarn(varn(pol)) | evallgef(4);
! 1367: x_r[3] = un;
1.1 noro 1368: for (m=1,i=1; i<lg(rep); i++)
1369: {
1.2 ! noro 1370: GEN q, r = (GEN)rep[i];
1.1 noro 1371:
1.2 ! noro 1372: r = nf_bestlift_to_pol(gmul(lt,r), NULL, &L);
! 1373: r = gdiv(r,lt);
! 1374: x_r[2] = lneg(r); /* check P(r) == 0 */
! 1375: q = RXQX_divrem(pol, x_r, nfpol, ONLY_DIVIDES);
! 1376: if (q) { pol = q; rep[m++] = (long)r; }
! 1377: else if (fl == 2) return cgetg(1, t_VEC);
1.1 noro 1378: }
1.2 ! noro 1379: rep[0] = evaltyp(t_VEC) | evallg(m);
! 1380: return gerepilecopy(av, rep);
1.1 noro 1381: }
1382:
1.2 ! noro 1383: T.polbase = polbase;
! 1384: T.pol = pol;
! 1385: T.nf = nf;
! 1386: T.fact = (GEN)factmod0(polred,p)[1];
! 1387: T.hint = 1; /* useless */
1.1 noro 1388:
1.2 ! noro 1389: rep = nf_combine_factors(&T, polred, p, L.k, dpol-1);
! 1390: return gerepileupto(av, rep);
1.1 noro 1391: }
1392:
1.2 ! noro 1393: /* return the characteristic polynomial of alpha over nf, where alpha
1.1 noro 1394: is an element of the algebra nf[X]/(T) given as a polynomial in X */
1395: GEN
1396: rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
1397: {
1.2 ! noro 1398: long vnf, vT, lT;
! 1399: gpmem_t av = avma;
1.1 noro 1400: GEN p1;
1401:
1402: nf=checknf(nf); vnf = varn(nf[1]);
1403: if (v<0) v = 0;
1404: T = fix_relative_pol(nf,T,1);
1405: if (typ(alpha) == t_POLMOD) alpha = lift_to_pol(alpha);
1406: lT = lgef(T);
1407: if (typ(alpha) != t_POL || varn(alpha) == vnf)
1408: return gerepileupto(av, gpowgs(gsub(polx[v], alpha), lT - 3));
1409: vT = varn(T);
1410: if (varn(alpha) != vT || v >= vnf)
1411: err(talker,"incorrect variables in rnfcharpoly");
1412: if (lgef(alpha) >= lT) alpha = gmod(alpha,T);
1413: if (lT <= 4)
1414: return gerepileupto(av, gsub(polx[v], alpha));
1415: p1 = caract2(unifpol(nf,T,1), unifpol(nf,alpha,1), v);
1416: return gerepileupto(av, unifpol(nf,p1,1));
1417: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to nffactor.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / modules