Annotation of OpenXM_contrib/pari-2.2/src/modules/nffactor.c, Revision 1.1
1.1 ! noro 1: /* $Id: nffactor.c,v 1.29 2001/10/01 12:11:33 karim Exp $
! 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:
! 24: extern GEN hensel_lift_fact(GEN pol, GEN fact, GEN T, GEN p, GEN pev, long e);
! 25: extern GEN nf_get_T2(GEN base, GEN polr);
! 26: extern GEN nfreducemodpr_i(GEN x, GEN prh);
! 27: extern GEN sort_factor(GEN y, int (*cmp)(GEN,GEN));
! 28: extern GEN pidealprimeinv(GEN nf, GEN x);
! 29:
! 30: static GEN nffactormod2(GEN nf,GEN pol,GEN pr);
! 31: static GEN nfmod_split2(GEN nf, GEN prhall, GEN polb, GEN v, GEN q);
! 32: static GEN nf_pol_mul(GEN nf,GEN pol1,GEN pol2);
! 33: static GEN nf_pol_divres(GEN nf,GEN pol1,GEN pol2, GEN *pr);
! 34: static GEN nf_pol_subres(GEN nf,GEN pol1,GEN pol2);
! 35: static GEN nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol);
! 36: static GEN nfmod_pol_divres(GEN nf,GEN prhall,GEN pol1,GEN pol2, GEN *pr);
! 37: static GEN nfmod_pol_gcd(GEN nf,GEN prhall,GEN pol1,GEN pol2);
! 38: static GEN nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt);
! 39: static GEN nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol);
! 40: static GEN nfsqff(GEN nf,GEN pol,long fl);
! 41: static int nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint);
! 42:
! 43: typedef struct Nfcmbf /* for use in nfsqff */
! 44: {
! 45: GEN pol, h, hinv, fact, res, lt, den;
! 46: long nfact, nfactmod;
! 47: } Nfcmbf;
! 48: static Nfcmbf nfcmbf;
! 49:
! 50: static GEN
! 51: unifpol0(GEN nf,GEN pol,long flag)
! 52: {
! 53: static long n = 0;
! 54: static GEN vun = NULL;
! 55: GEN f = (GEN) nf[1];
! 56: long i = degpol(f), av;
! 57:
! 58: if (i != n)
! 59: {
! 60: n=i; if (vun) free(vun);
! 61: vun = (GEN) gpmalloc((n+1)*sizeof(long));
! 62: vun[0] = evaltyp(t_COL) | evallg(n+1); vun[1]=un;
! 63: for ( ; i>=2; i--) vun[i]=zero;
! 64: }
! 65:
! 66: av = avma;
! 67: switch(typ(pol))
! 68: {
! 69: case t_INT: case t_FRAC: case t_RFRAC:
! 70: pol = gmul(pol,vun); break;
! 71:
! 72: case t_POL:
! 73: pol = gmodulcp(pol,f); /* fall through */
! 74: case t_POLMOD:
! 75: pol = algtobasis(nf,pol);
! 76: }
! 77: if (flag) pol = basistoalg(nf, lift(pol));
! 78: return gerepileupto(av,pol);
! 79: }
! 80:
! 81: /* Let pol be a polynomial with coefficients in Z or nf (vectors or polymods)
! 82: * return the same polynomial with coefficients expressed:
! 83: * if flag=0: as vectors (on the integral basis).
! 84: * if flag=1: as polymods.
! 85: */
! 86: GEN
! 87: unifpol(GEN nf,GEN pol,long flag)
! 88: {
! 89: if (typ(pol)==t_POL && varn(pol) < varn(nf[1]))
! 90: {
! 91: long i, d = lgef(pol);
! 92: GEN p1 = pol;
! 93:
! 94: pol=cgetg(d,t_POL); pol[1]=p1[1];
! 95: for (i=2; i<d; i++)
! 96: pol[i] = (long) unifpol0(nf,(GEN) p1[i],flag);
! 97:
! 98: return pol;
! 99: }
! 100: return unifpol0(nf,(GEN) pol, flag);
! 101: }
! 102:
! 103: #if 0
! 104: /* return a monic polynomial of degree d with random coefficients in Z_nf */
! 105: static GEN
! 106: random_pol(GEN nf,long d)
! 107: {
! 108: long i,j, n = degpol(nf[1]);
! 109: GEN pl,p;
! 110:
! 111: pl=cgetg(d+3,t_POL);
! 112: for (i=2; i<d+2; i++)
! 113: {
! 114: p=cgetg(n+1,t_COL); pl[i]=(long)p;
! 115: for (j=1; j<=n; j++)
! 116: p[j] = lstoi(mymyrand()%101 - 50);
! 117: }
! 118: p=cgetg(n+1,t_COL); pl[i]=(long)p;
! 119: p[1]=un; for (i=2; i<=n; i++) p[i]=zero;
! 120:
! 121: pl[1] = evalsigne(1) | evallgef(d+3) | evalvarn(0);
! 122: return pl;
! 123: }
! 124: #endif
! 125:
! 126: /* multiplication of x by y */
! 127: static GEN
! 128: nf_pol_mul(GEN nf,GEN x,GEN y)
! 129: {
! 130: long tetpil,av=avma;
! 131: GEN res = gmul(unifpol(nf,x,1), unifpol(nf,y,1));
! 132:
! 133: tetpil = avma;
! 134: return gerepile(av,tetpil,unifpol(nf,res,0));
! 135: }
! 136:
! 137: /* compute x^2 in nf */
! 138: static GEN
! 139: nf_pol_sqr(GEN nf,GEN x)
! 140: {
! 141: long tetpil,av=avma;
! 142: GEN res = gsqr(unifpol(nf,x,1));
! 143:
! 144: tetpil = avma;
! 145: return gerepile(av,tetpil,unifpol(nf,res,0));
! 146: }
! 147:
! 148: /* reduce the coefficients of pol modulo prhall */
! 149: static GEN
! 150: nfmod_pol_reduce(GEN nf,GEN prhall,GEN pol)
! 151: {
! 152: long av=avma,tetpil,i;
! 153: GEN p1;
! 154:
! 155: if (typ(pol)!=t_POL) return nfreducemodpr(nf,pol,prhall);
! 156: pol=unifpol(nf,pol,0);
! 157:
! 158: tetpil=avma; i=lgef(pol);
! 159: p1=cgetg(i,t_POL); p1[1]=pol[1];
! 160: for (i--; i>=2; i--)
! 161: p1[i] = (long) nfreducemodpr(nf,(GEN)pol[i],prhall);
! 162: return gerepile(av,tetpil, normalizepol(p1));
! 163: }
! 164:
! 165: /* x^2 modulo prhall */
! 166: static GEN
! 167: nfmod_pol_sqr(GEN nf,GEN prhall,GEN x)
! 168: {
! 169: long av=avma,tetpil;
! 170: GEN px;
! 171:
! 172: px = nfmod_pol_reduce(nf,prhall,x);
! 173: px = unifpol(nf,lift(px),1);
! 174: px = unifpol(nf,nf_pol_sqr(nf,px),0);
! 175: tetpil=avma;
! 176: return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));
! 177: }
! 178:
! 179: /* multiplication of x by y modulo prhall */
! 180: static GEN
! 181: nfmod_pol_mul(GEN nf,GEN prhall,GEN x,GEN y)
! 182: {
! 183: long av=avma,tetpil;
! 184: GEN px,py;
! 185:
! 186: px = nfmod_pol_reduce(nf,prhall,x); px = unifpol(nf,lift(px),1);
! 187: py = nfmod_pol_reduce(nf,prhall,y); py = unifpol(nf,lift(py),1);
! 188: px = unifpol(nf,nf_pol_mul(nf,px,py),0);
! 189: tetpil=avma;
! 190: return gerepile(av,tetpil,nfmod_pol_reduce(nf,prhall,px));
! 191: }
! 192:
! 193: /* Euclidan division of x by y */
! 194: static GEN
! 195: nf_pol_divres(GEN nf,GEN x,GEN y,GEN *pr)
! 196: {
! 197: long av = avma,tetpil;
! 198: GEN nq = poldivres(unifpol(nf,x,1),unifpol(nf,y,1),pr);
! 199: GEN *gptr[2];
! 200:
! 201: tetpil=avma; nq=unifpol(nf,nq,0);
! 202: if (pr) *pr = unifpol(nf,*pr,0);
! 203: gptr[0]=&nq; gptr[1]=pr;
! 204: gerepilemanysp(av,tetpil,gptr,pr ? 2:1);
! 205: return nq;
! 206: }
! 207:
! 208: /* Euclidan division of x by y modulo prhall */
! 209: static GEN
! 210: nfmod_pol_divres(GEN nf,GEN prhall,GEN x,GEN y, GEN *pr)
! 211: {
! 212: long av=avma,dx,dy,dz,i,j,k,l,n,tetpil;
! 213: GEN z,p1,p3,px,py;
! 214:
! 215: py = nfmod_pol_reduce(nf,prhall,y);
! 216: if (gcmp0(py))
! 217: err(talker, "division by zero in nfmod_pol_divres");
! 218:
! 219: tetpil=avma;
! 220: px=nfmod_pol_reduce(nf,prhall,x);
! 221: dx=degpol(px); dy=degpol(py); dz=dx-dy;
! 222: if (dx<dy)
! 223: {
! 224: GEN vzero;
! 225:
! 226: if (pr) *pr = gerepile(av,tetpil,px);
! 227: else avma = av;
! 228:
! 229: n=degpol(nf[1]);
! 230: vzero = cgetg(n+1,t_COL);
! 231: n=degpol(nf[1]);
! 232: for (i=1; i<=n; i++) vzero[i]=zero;
! 233:
! 234: z=cgetg(3,t_POL); z[2]=(long)vzero;
! 235: z[1]=evallgef(2) | evalvarn(varn(px));
! 236: return z;
! 237: }
! 238: p1 = NULL; /* gcc -Wall */
! 239:
! 240: z=cgetg(dz+3,t_POL); z[1]=evalsigne(1) | evallgef(3+dz);
! 241: setvarn(z,varn(px));
! 242: z[dz+2] = (long) element_divmodpr(nf,(GEN)px[dx+2],(GEN)py[dy+2],prhall);
! 243: for (i=dx-1; i>=dy; --i)
! 244: {
! 245: l=avma; p1=(GEN)px[i+2];
! 246: for (j=i-dy+1; j<=i && j<=dz; j++)
! 247: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));
! 248: p1 = nfreducemodpr(nf,p1,prhall);
! 249: tetpil=avma; p3=element_divmodpr(nf,p1,(GEN)py[dy+2],prhall);
! 250: z[i-dy+2]=lpile(l,tetpil,p3);
! 251: z[i-dy+2]=(long)nfreducemodpr(nf,(GEN)z[i-dy+2],prhall);
! 252: }
! 253: l=avma;
! 254: for (i=dy-1; i>=0; --i)
! 255: {
! 256: l=avma; p1=((GEN)px[i+2]);
! 257: for (j=0; j<=i && j<=dz; j++)
! 258: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[i-j+2]));
! 259: p1 = gerepileupto(l, nfreducemodpr(nf,p1,prhall));
! 260: if (!gcmp0(p1)) break;
! 261: }
! 262:
! 263: if (!pr) { avma = l; return z; }
! 264:
! 265: if (i<0)
! 266: {
! 267: avma=l;
! 268: p3 = cgetg(3,t_POL); p3[2]=zero;
! 269: p3[1] = evallgef(2) | evalvarn(varn(px));
! 270: *pr=p3; return z;
! 271: }
! 272:
! 273: p3=cgetg(i+3,t_POL);
! 274: p3[1]=evalsigne(1) | evallgef(i+3) | evalvarn(varn(px));
! 275: p3[i+2]=(long)nfreducemodpr(nf,p1,prhall);
! 276: for (k=i-1; k>=0; --k)
! 277: {
! 278: l=avma; p1=((GEN)px[k+2]);
! 279: for (j=0; j<=k && j<=dz; j++)
! 280: p1 = gsub(p1, element_mul(nf,(GEN)z[j+2],(GEN)py[k-j+2]));
! 281: p3[k+2]=lpileupto(l,nfreducemodpr(nf,p1,prhall));
! 282: }
! 283: *pr=p3; return z;
! 284: }
! 285:
! 286: /* GCD of x and y */
! 287: static GEN
! 288: nf_pol_subres(GEN nf,GEN x,GEN y)
! 289: {
! 290: long av=avma,tetpil;
! 291: GEN s = srgcd(unifpol(nf,x,1), unifpol(nf,y,1));
! 292:
! 293: tetpil=avma; return gerepile(av,tetpil,unifpol(nf,s,1));
! 294: }
! 295:
! 296: /* GCD of x and y modulo prhall */
! 297: static GEN
! 298: nfmod_pol_gcd(GEN nf,GEN prhall,GEN x,GEN y)
! 299: {
! 300: long av=avma;
! 301: GEN p1,p2;
! 302:
! 303: if (lgef(x)<lgef(y)) { p1=y; y=x; x=p1; }
! 304: p1=nfmod_pol_reduce(nf,prhall,x);
! 305: p2=nfmod_pol_reduce(nf,prhall,y);
! 306: while (!isexactzero(p2))
! 307: {
! 308: GEN p3;
! 309:
! 310: nfmod_pol_divres(nf,prhall,p1,p2,&p3);
! 311: p1=p2; p2=p3;
! 312: }
! 313: return gerepileupto(av,p1);
! 314: }
! 315:
! 316: /* return pol^e modulo prhall and the polynomial pmod */
! 317: static GEN
! 318: nfmod_pol_pow(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN e)
! 319: {
! 320: long i, av = avma, n = degpol(nf[1]);
! 321: GEN p1,p2,vun;
! 322:
! 323: vun=cgetg(n+1,t_COL); vun[1]=un; for (i=2; i<=n; i++) vun[i]=zero;
! 324: p1=gcopy(polun[varn(pol)]); p1[2]=(long)vun;
! 325: if (gcmp0(e)) return p1;
! 326:
! 327: p2=nfmod_pol_reduce(nf,prhall,pol);
! 328: for(;;)
! 329: {
! 330: if (!vali(e))
! 331: {
! 332: p1=nfmod_pol_mul(nf,prhall,p1,p2);
! 333: nfmod_pol_divres(nf,prhall,p1,pmod,&p1);
! 334: }
! 335: if (gcmp1(e)) break;
! 336:
! 337: e=shifti(e,-1);
! 338: p2=nfmod_pol_sqr(nf,prhall,p2);
! 339: nfmod_pol_divres(nf,prhall,p2,pmod,&p2);
! 340: }
! 341: return gerepileupto(av,p1);
! 342: }
! 343:
! 344: static long
! 345: isdivbyprime(GEN nf, GEN x, GEN pr)
! 346: {
! 347: GEN elt, p = (GEN)pr[1], tau = (GEN)pr[5];
! 348:
! 349: elt = element_mul(nf, x, tau);
! 350: if (divise(content(elt), p)) return 1;
! 351:
! 352: return 0;
! 353: }
! 354:
! 355: /* return the factor of nf.pol modulo p corresponding to pr */
! 356: static GEN
! 357: localpol(GEN nf, GEN pr)
! 358: {
! 359: long i, l;
! 360: GEN fct, pol = (GEN)nf[1], p = (GEN)pr[1];
! 361:
! 362: fct = lift((GEN)factmod(pol, p)[1]);
! 363: l = lg(fct) - 1;
! 364: for (i = 1; i <= l; i++)
! 365: if (isdivbyprime(nf, (GEN)fct[i], pr)) return (GEN)fct[i];
! 366:
! 367: err(talker, "cannot find a suitable factor in localpol");
! 368: return NULL; /* not reached */
! 369: }
! 370:
! 371: /* factorization of x modulo pr */
! 372: static GEN
! 373: nffactormod0(GEN nf, GEN x, GEN pr)
! 374: {
! 375: long av = avma, j, l, vx = varn(x), vn;
! 376: GEN rep, pol, xrd, prh, p1;
! 377:
! 378: nf=checknf(nf);
! 379: vn = varn((GEN)nf[1]);
! 380: if (typ(x)!=t_POL) err(typeer,"nffactormod");
! 381: if (vx >= vn)
! 382: err(talker,"polynomial variable must have highest priority in nffactormod");
! 383:
! 384: prh = nfmodprinit(nf, pr);
! 385:
! 386: if (divise((GEN)nf[4], (GEN)pr[1]))
! 387: return gerepileupto(av, nffactormod2(nf,x,pr));
! 388:
! 389: xrd = nfmod_pol_reduce(nf, prh, x);
! 390: if (gcmp1((GEN)pr[4]))
! 391: {
! 392: pol = gun; /* dummy */
! 393: for (j = 2; j < lg(xrd); j++)
! 394: xrd[j] = mael(xrd, j, 1);
! 395: rep = factmod(xrd, (GEN)pr[1]);
! 396: rep = lift(rep);
! 397: }
! 398: else
! 399: {
! 400: pol = localpol(nf, pr);
! 401: xrd = lift(unifpol(nf, xrd, 1));
! 402: p1 = gun;
! 403: for (j = 2; j < lg(xrd); j++)
! 404: {
! 405: xrd[j] = lmod((GEN)xrd[j], pol);
! 406: p1 = mpppcm(p1, denom(content((GEN)xrd[j])));
! 407: }
! 408: rep = factmod9(gmul(xrd, p1), (GEN)pr[1], pol);
! 409: rep = lift(lift(rep));
! 410: }
! 411:
! 412: l = lg((GEN)rep[1]);
! 413: for (j = 1; j < l; j++)
! 414: mael(rep, 1, j) = (long)unifpol(nf, gmael(rep, 1, j), 1);
! 415:
! 416: return gerepilecopy(av, rep);
! 417: }
! 418:
! 419: GEN
! 420: nffactormod(GEN nf, GEN x, GEN pr)
! 421: {
! 422: return nffactormod0(nf, x, pr);
! 423: }
! 424:
! 425: extern GEN trivfact(void);
! 426:
! 427: /* factorization of x modulo pr */
! 428: GEN
! 429: nffactormod2(GEN nf,GEN pol,GEN pr)
! 430: {
! 431: long av = avma, tetpil,lb,nbfact,psim,N,n,i,j,k,d,e,vf,r,kk;
! 432: GEN y,ex,*t,f1,f2,f3,df1,g1,polb,pold,polu,vker;
! 433: GEN Q,f,x,u,v,v2,v3,vz,q,vun,vzero,prhall;
! 434:
! 435: nf=checknf(nf);
! 436: if (typ(pol)!=t_POL) err(typeer,"nffactormod");
! 437: if (varn(pol) >= varn(nf[1]))
! 438: err(talker,"polynomial variable must have highest priority in nffactormod");
! 439:
! 440: prhall=nfmodprinit(nf,pr); n=degpol(nf[1]);
! 441: vun = gscalcol_i(gun, n);
! 442: vzero = gscalcol_i(gzero, n);
! 443: q = vker = NULL; /* gcc -Wall */
! 444:
! 445: f=unifpol(nf,pol,0); f=nfmod_pol_reduce(nf,prhall,f);
! 446: if (!signe(f)) err(zeropoler,"nffactormod");
! 447: d=degpol(f); vf=varn(f);
! 448: if (d == 0) return trivfact();
! 449: t = (GEN*)new_chunk(d+1); ex = cgetg(d+1, t_VECSMALL);
! 450: x=gcopy(polx[vf]); x[3]=(long)vun; x[2]=(long)vzero;
! 451: if (d<=1) { nbfact=2; t[1] = f; ex[1]=1; }
! 452: else
! 453: {
! 454: q = (GEN)pr[1]; psim = VERYBIGINT;
! 455: if (!is_bigint(q)) psim = itos(q);
! 456: /* psim has an effect only when p is small. If too big, set it to a huge
! 457: * number (i.e ignore it) to avoid an error in itos on next line.
! 458: */
! 459: q=gpuigs(q, itos((GEN)pr[4]));
! 460: f1=f; e=1; nbfact=1;
! 461: while (lgef(f1)>3)
! 462: {
! 463: df1=deriv(f1,vf); f2=nfmod_pol_gcd(nf,prhall,f1,df1);
! 464: g1=nfmod_pol_divres(nf,prhall,f1,f2,NULL); k=0;
! 465: while (lgef(g1)>3)
! 466: {
! 467: k++;
! 468: if (k%psim == 0)
! 469: {
! 470: k++; f2=nfmod_pol_divres(nf,prhall,f2,g1,NULL);
! 471: }
! 472: f3=nfmod_pol_gcd(nf,prhall,f2,g1);
! 473: u = nfmod_pol_divres(nf,prhall,g1,f3,NULL);
! 474: f2= nfmod_pol_divres(nf,prhall,f2,f3,NULL);
! 475: g1=f3;
! 476: if (lgef(u)>3)
! 477: {
! 478: N=degpol(u); Q=cgetg(N+1,t_MAT);
! 479: v3=cgetg(N+1,t_COL); Q[1]=(long)v3;
! 480: v3[1]=(long)vun; for (i=2; i<=N; i++) v3[i]=(long)vzero;
! 481:
! 482: v2 = v = nfmod_pol_pow(nf,prhall,u,x,q);
! 483: for (j=2; j<=N; j++)
! 484: {
! 485: v3=cgetg(N+1,t_COL); Q[j]=(long)v3;
! 486: for (i=1; i<=lgef(v2)-2; i++) v3[i]=v2[i+1];
! 487: for (; i<=N; i++) v3[i]=(long)vzero;
! 488: if (j<N)
! 489: {
! 490: v2=nfmod_pol_mul(nf,prhall,v2,v);
! 491: nfmod_pol_divres(nf,prhall,v2,u,&v2);
! 492: }
! 493: }
! 494: for (i=1; i<=N; i++)
! 495: coeff(Q,i,i)=lsub((GEN)coeff(Q,i,i),vun);
! 496: v2=nfkermodpr(nf,Q,prhall); r=lg(v2)-1; t[nbfact]=gcopy(u); kk=1;
! 497: if (r>1)
! 498: {
! 499: vker=cgetg(r+1,t_COL);
! 500: for (i=1; i<=r; i++)
! 501: {
! 502: v3=cgetg(N+2,t_POL);
! 503: v3[1]=evalsigne(1)+evallgef(2+N); setvarn(v3,vf);
! 504: vker[i]=(long)v3; for (j=1; j<=N; j++) v3[j+1]=coeff(v2,j,i);
! 505: normalizepol(v3);
! 506: }
! 507: }
! 508: while (kk<r)
! 509: {
! 510: v=gcopy(polun[vf]); v[2]=(long)vzero;
! 511: for (i=1; i<=r; i++)
! 512: {
! 513: vz=cgetg(n+1,t_COL);
! 514: for (j=1; j<=n; j++)
! 515: vz[j] = lmodsi(mymyrand()>>8, q);
! 516: vz=nfreducemodpr(nf,vz,prhall);
! 517: v=gadd(v,nfmod_pol_mul(nf,prhall,vz,(GEN)vker[i]));
! 518: }
! 519: for (i=1; i<=kk && kk<r; i++)
! 520: {
! 521: polb=t[nbfact+i-1]; lb=lgef(polb);
! 522: if (lb>4)
! 523: {
! 524: if(psim==2)
! 525: polu=nfmod_split2(nf,prhall,polb,v,q);
! 526: else
! 527: {
! 528: polu=nfmod_pol_pow(nf,prhall,polb,v,shifti(q,-1));
! 529: polu=gsub(polu,vun);
! 530: }
! 531: pold=nfmod_pol_gcd(nf,prhall,polb,polu);
! 532: if (lgef(pold)>3 && lgef(pold)<lb)
! 533: {
! 534: t[nbfact+i-1]=pold; kk++;
! 535: t[nbfact+kk-1]=nfmod_pol_divres(nf,prhall,polb,pold,NULL);
! 536: }
! 537: }
! 538: }
! 539: }
! 540: for (i=nbfact; i<nbfact+r; i++) ex[i]=e*k;
! 541: nbfact+=r;
! 542: }
! 543: }
! 544: e*=psim; j=(degpol(f2))/psim+3; f1=cgetg(j,t_POL);
! 545: f1[1] = evalsigne(1) | evallgef(j) | evalvarn(vf);
! 546: for (i=2; i<j; i++)
! 547: f1[i]=(long)element_powmodpr(nf,(GEN)f2[psim*(i-2)+2],
! 548: gdiv(q,(GEN)pr[1]),prhall);
! 549: }
! 550: }
! 551: if (nbfact < 2)
! 552: err(talker, "%Z not a prime (nffactormod)", q);
! 553: for (j=1; j<nbfact; j++)
! 554: {
! 555: v = element_divmodpr(nf,vun,gmael(t,j,lgef(t[j])-1),prhall);
! 556: t[j] = unifpol(nf,nfmod_pol_mul(nf,prhall,v,(GEN)t[j]),1);
! 557: }
! 558:
! 559: tetpil=avma; y=cgetg(3,t_MAT);
! 560: u=cgetg(nbfact,t_COL); y[1]=(long)u;
! 561: v=cgetg(nbfact,t_COL); y[2]=(long)v;
! 562: for (j=1,k=0; j<nbfact; j++)
! 563: if (ex[j])
! 564: {
! 565: k++;
! 566: u[k]=lcopy((GEN)t[j]);
! 567: v[k]=lstoi(ex[j]);
! 568: }
! 569: return gerepile(av,tetpil,y);
! 570: }
! 571:
! 572: /* return pol + pol^2 + ... + pol^(q/2) modulo prhall and
! 573: the polynomial pmod */
! 574: static GEN
! 575: nfmod_split2(GEN nf,GEN prhall,GEN pmod,GEN pol,GEN exp)
! 576: {
! 577: long av = avma;
! 578: GEN p1,p2,q;
! 579:
! 580: if (cmpis(exp,2)<=0) return pol;
! 581: p2=p1=pol; q=shifti(exp,-1);
! 582: while (!gcmp1(q))
! 583: {
! 584: p2=nfmod_pol_sqr(nf,prhall,p2);
! 585: nfmod_pol_divres(nf,prhall,p2,pmod,&p2);
! 586: q=shifti(q,-1); p1=gadd(p1,p2);
! 587: }
! 588: return gerepileupto(av,p1);
! 589: }
! 590:
! 591: /* If p doesn't divide either a or b and has a divisor of degree 1, return it.
! 592: * Return NULL otherwise.
! 593: */
! 594: static GEN
! 595: p_ok(GEN nf, GEN p, GEN a)
! 596: {
! 597: long av,m,i;
! 598: GEN dec;
! 599:
! 600: if (divise(a,p)) return NULL;
! 601: av = avma; dec = primedec(nf,p); m=lg(dec);
! 602: for (i=1; i<m; i++)
! 603: {
! 604: GEN pr = (GEN)dec[i];
! 605: if (is_pm1(pr[4]))
! 606: return pr;
! 607: }
! 608: avma = av; return NULL;
! 609: }
! 610:
! 611: /* for each new prime ct--, if ct = 0, return NULL */
! 612: static GEN
! 613: choose_prime(GEN nf, GEN dk, GEN lim, long ct)
! 614: {
! 615: GEN p, pr;
! 616:
! 617: p = nextprime(lim);
! 618: for (;;)
! 619: {
! 620: if ((pr = p_ok(nf,p,dk))) break;
! 621: ct--;
! 622: if (!ct) return NULL;
! 623: p = nextprime(addis(p,2));
! 624: }
! 625:
! 626: return pr;
! 627: }
! 628:
! 629: /* test if the discriminant of polbase modulo some few primes
! 630: is non-zero. Return 1 if it is so (=> polbase is square-free)
! 631: and 0 otherwise (=> polbase may or may not be square-free) */
! 632: static int
! 633: is_sqf(GEN nf, GEN polbase)
! 634: {
! 635: GEN lt, pr, prh, p2, p;
! 636: long i, d = lgef(polbase), ct = 5;
! 637:
! 638: lt = (GEN)leading_term(polbase)[1];
! 639: p = stoi(101);
! 640:
! 641: while (ct > 0)
! 642: {
! 643: /* small primes tend to divide discriminants more often
! 644: than large ones so we look at primes >= 101 */
! 645: pr = choose_prime(nf,lt,p,30);
! 646: if (!pr) break;
! 647:
! 648: p=(GEN)pr[1];
! 649: prh=prime_to_ideal(nf,pr);
! 650:
! 651: p2=gcopy(polbase);
! 652: lt=mpinvmod(lt,p);
! 653:
! 654: for (i=2; i<d; i++)
! 655: p2[i] = nfreducemodpr_i(gmul(lt,(GEN)p2[i]), prh)[1];
! 656: p2 = normalizepol(p2);
! 657:
! 658: /* discriminant is non-zero => polynomial is square-free */
! 659: if (!gcmp0(p2) && !divise(discsr(p2),p)) { return 1; }
! 660:
! 661: ct--;
! 662: p=addis(p,1);
! 663: }
! 664:
! 665: return 0;
! 666: }
! 667:
! 668: /* rescale p in K[X] (coeffs in algtobasis form) --> primitive in O_K[X] */
! 669: GEN
! 670: nf_pol_to_int(GEN p, GEN *den)
! 671: {
! 672: long i, l = lgef(p);
! 673: GEN d = gun;
! 674: for (i=2; i<l; i++)
! 675: d = mpppcm(d,denom((GEN)p[i]));
! 676: if (! gcmp1(d)) p = gmul(p, d);
! 677: *den = d; return p;
! 678: }
! 679:
! 680: /* return the roots of pol in nf */
! 681: GEN
! 682: nfroots(GEN nf,GEN pol)
! 683: {
! 684: long av=avma,tetpil,d=lgef(pol),fl;
! 685: GEN p1,p2,polbase,polmod,den;
! 686:
! 687: p2=NULL; /* gcc -Wall */
! 688: nf=checknf(nf);
! 689: if (typ(pol)!=t_POL) err(talker,"not a polynomial in nfroots");
! 690: if (varn(pol) >= varn(nf[1]))
! 691: err(talker,"polynomial variable must have highest priority in nfroots");
! 692:
! 693: polbase=unifpol(nf,pol,0);
! 694:
! 695: if (d==3)
! 696: {
! 697: tetpil=avma; p1=cgetg(1,t_VEC);
! 698: return gerepile(av,tetpil,p1);
! 699: }
! 700:
! 701: if (d==4)
! 702: {
! 703: tetpil=avma; p1=cgetg(2,t_VEC);
! 704: p1[1] = (long)basistoalg(nf,gneg_i(
! 705: element_div(nf,(GEN)polbase[2],(GEN)polbase[3])));
! 706: return gerepile(av,tetpil,p1);
! 707: }
! 708:
! 709: p1=element_inv(nf,leading_term(polbase));
! 710: polbase=nf_pol_mul(nf,p1,polbase);
! 711:
! 712: polbase = nf_pol_to_int(polbase, &den);
! 713: polmod=unifpol(nf,polbase,1);
! 714:
! 715: if (DEBUGLEVEL>=4)
! 716: fprintferr("test if the polynomial is square-free\n");
! 717:
! 718: fl = is_sqf(nf, polbase);
! 719:
! 720: /* the polynomial may not be square-free ... */
! 721: if (!fl)
! 722: {
! 723: p1=derivpol(polmod);
! 724: p2=nf_pol_subres(nf,polmod,p1);
! 725: if (degpol(p2) == 0) fl = 1;
! 726: }
! 727:
! 728: if (!fl)
! 729: {
! 730: p1=element_inv(nf,leading_term(p2));
! 731: p2=nf_pol_mul(nf,p1,p2);
! 732: polmod=nf_pol_divres(nf,polmod,p2,NULL);
! 733:
! 734: p1=element_inv(nf,leading_term(polmod));
! 735: polmod=nf_pol_mul(nf,p1,polmod);
! 736:
! 737: polmod = nf_pol_to_int(polmod, &den);
! 738: polmod=unifpol(nf,polmod,1);
! 739: }
! 740:
! 741: p1 = nfsqff(nf,polmod,1);
! 742: tetpil=avma; return gerepile(av, tetpil, gen_sort(p1, 0, cmp_pol));
! 743: }
! 744:
! 745: /* return a minimal lift of elt modulo id */
! 746: static GEN
! 747: nf_bestlift(GEN id,GEN idinv,GEN den,GEN elt)
! 748: {
! 749: GEN u = gmul(idinv,elt);
! 750: long i, l = lg(u);
! 751: for (i=1; i<l; i++) u[i] = (long)gdivround((GEN)u[i], den);
! 752: return gsub(elt, gmul(id, u));
! 753: }
! 754:
! 755: /* return the lift of pol with coefficients of T2-norm <= C (if possible) */
! 756: static GEN
! 757: nf_pol_lift(GEN id,GEN idinv,GEN den,GEN pol)
! 758: {
! 759: long i, d = lgef(pol);
! 760: GEN x = cgetg(d,t_POL);
! 761: x[1] = pol[1];
! 762: for (i=2; i<d; i++)
! 763: x[i] = (long) nf_bestlift(id,idinv,den,(GEN)pol[i]);
! 764: return x;
! 765: }
! 766:
! 767: #if 0
! 768: /* return pol(elt) */
! 769: static GEN
! 770: nf_pol_eval(GEN nf,GEN pol,GEN elt)
! 771: {
! 772: long av=avma,tetpil,i;
! 773: GEN p1;
! 774:
! 775: i=lgef(pol)-1; if (i==2) return gcopy((GEN)pol[2]);
! 776:
! 777: p1=element_mul(nf,(GEN)pol[i],elt);
! 778: for (i--; i>=3; i--)
! 779: p1=element_mul(nf,elt,gadd((GEN)pol[i],p1));
! 780: tetpil=avma; return gerepile(av,tetpil,gadd(p1,(GEN)pol[2]));
! 781: }
! 782: #endif
! 783:
! 784: /* return the factorization of x in nf */
! 785: GEN
! 786: nffactor(GEN nf,GEN pol)
! 787: {
! 788: GEN y,p1,p2,den,p3,p4,quot,rep=cgetg(3,t_MAT);
! 789: long av = avma,tetpil,i,j,d,fl;
! 790:
! 791: if (DEBUGLEVEL >= 4) timer2();
! 792:
! 793: p3=NULL; /* gcc -Wall */
! 794: nf=checknf(nf);
! 795: if (typ(pol)!=t_POL) err(typeer,"nffactor");
! 796: if (varn(pol) >= varn(nf[1]))
! 797: err(talker,"polynomial variable must have highest priority in nffactor");
! 798:
! 799: d=lgef(pol);
! 800: if (d==3)
! 801: {
! 802: rep[1]=lgetg(1,t_COL);
! 803: rep[2]=lgetg(1,t_COL);
! 804: return rep;
! 805: }
! 806: if (d==4)
! 807: {
! 808: p1=cgetg(2,t_COL); rep[1]=(long)p1; p1[1]=lcopy(pol);
! 809: p1=cgetg(2,t_COL); rep[2]=(long)p1; p1[1]=un;
! 810: return rep;
! 811: }
! 812:
! 813: p1=element_inv(nf,leading_term(pol));
! 814: pol=nf_pol_mul(nf,p1,pol);
! 815:
! 816: p1=unifpol(nf,pol,0);
! 817: p1 = nf_pol_to_int(p1, &den);
! 818:
! 819: if (DEBUGLEVEL>=4)
! 820: fprintferr("test if the polynomial is square-free\n");
! 821:
! 822: fl = is_sqf(nf, p1);
! 823:
! 824: /* polynomial may not be square-free ... */
! 825: if (!fl)
! 826: {
! 827: p2=derivpol(p1);
! 828: p3=nf_pol_subres(nf,p1,p2);
! 829: if (degpol(p3) == 0) fl = 1;
! 830: }
! 831:
! 832: if (!fl)
! 833: {
! 834: p4=element_inv(nf,leading_term(p3));
! 835: p3=nf_pol_mul(nf,p4,p3);
! 836:
! 837: p2=nf_pol_divres(nf,p1,p3,NULL);
! 838: p4=element_inv(nf,leading_term(p2));
! 839: p2=nf_pol_mul(nf,p4,p2);
! 840:
! 841: p2 = nf_pol_to_int(p2, &den);
! 842:
! 843: p2=unifpol(nf,p2,1);
! 844: tetpil = avma; y = nfsqff(nf,p2,0);
! 845: i = nfcmbf.nfact;
! 846:
! 847: quot=nf_pol_divres(nf,p1,p2,NULL);
! 848: p3=(GEN)gpmalloc((i+1) * sizeof(long));
! 849: for (j=i; j>=1; j--)
! 850: {
! 851: GEN fact=(GEN)y[j], quo = quot, rem;
! 852: long e = 0;
! 853:
! 854: do
! 855: {
! 856: quo = nf_pol_divres(nf,quo,fact,&rem);
! 857: e++;
! 858: }
! 859: while (gcmp0(rem));
! 860: p3[j]=lstoi(e);
! 861: }
! 862: avma = (long)y; y = gerepile(av, tetpil, y);
! 863: p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=lcopy((GEN)p3[i]);
! 864: free(p3);
! 865: }
! 866: else
! 867: {
! 868: tetpil=avma; y = gerepile(av,tetpil,nfsqff(nf,p1,0));
! 869: i = nfcmbf.nfact;
! 870: p2=cgetg(i+1, t_COL); for (; i>=1; i--) p2[i]=un;
! 871: }
! 872: if (DEBUGLEVEL>=4)
! 873: fprintferr("number of factor(s) found: %ld\n", nfcmbf.nfact);
! 874: rep[1]=(long)y;
! 875: rep[2]=(long)p2; return sort_factor(rep, cmp_pol);
! 876: }
! 877:
! 878: /* test if the matrix M is suitable */
! 879: static long
! 880: test_mat(GEN M, GEN p, GEN C2, long k)
! 881: {
! 882: long av = avma, i, N = lg(M);
! 883: GEN min, prod, L2, R;
! 884:
! 885: min = prod = gcoeff(M,1,1);
! 886: for (i = 2; i < N; i++)
! 887: {
! 888: L2 = gcoeff(M,i,i); prod = mpmul(prod,L2);
! 889: if (mpcmp(L2,min) < 0) min = L2;
! 890: }
! 891: R = mpmul(min, gpowgs(p, k<<1));
! 892: i = mpcmp(mpmul(C2,prod), R); avma = av;
! 893: return (i < 0);
! 894: }
! 895:
! 896: /* return the matrix corresponding to pr^e with R(pr^e) > C */
! 897: static GEN
! 898: T2_matrix_pow(GEN nf, GEN pre, GEN p, GEN C, long *kmax, long prec)
! 899: {
! 900: long av=avma,av1,lim, k = *kmax, N = degpol((GEN)nf[1]);
! 901: int tot_real = !signe(gmael(nf,2,2));
! 902: GEN p1,p2,p3,u,C2,T2, x = (GEN)nf[1];
! 903:
! 904: C2 = gdiv(gmul2n(C,2), absi((GEN)nf[3]));
! 905: p1 = gmul(pre,lllintpartial(pre)); av1 = avma;
! 906: T2 = tot_real? gmael(nf,5,4)
! 907: : nf_get_T2((GEN) nf[7], roots(x,prec));
! 908: p3 = qf_base_change(T2,p1,1);
! 909:
! 910: if (N <= 6 && test_mat(p3,p,C2,k))
! 911: {
! 912: avma = av1; return gerepileupto(av,p1);
! 913: }
! 914:
! 915: av1=avma; lim = stack_lim(av1,2);
! 916: for (;;)
! 917: {
! 918: if (DEBUGLEVEL>2) fprintferr("exponent: %ld\n",k);
! 919:
! 920: for(;;)
! 921: {
! 922: u = tot_real? lllgramall(p3,2,lll_IM) : lllgramintern(p3,2,1,prec);
! 923: if (u) break;
! 924:
! 925: prec=(prec<<1)-2;
! 926: if (DEBUGLEVEL > 1) err(warnprec,"nffactor[1]",prec);
! 927: T2 = nf_get_T2((GEN) nf[7], roots(x,prec));
! 928: p3 = qf_base_change(T2,p1,1);
! 929: }
! 930: if (DEBUGLEVEL>2) msgtimer("lllgram + base change");
! 931: p3 = qf_base_change(p3,u,1);
! 932: if (test_mat(p3,p,C2,k))
! 933: {
! 934: *kmax = k;
! 935: return gerepileupto(av,gmul(p1,u));
! 936: }
! 937:
! 938: /* we also need to increase the precision */
! 939: p2=shifti(gceil(mulsr(k, glog(p,DEFAULTPREC))),-1);
! 940: prec += (long)(itos(p2)*pariK1);
! 941: if (DEBUGLEVEL > 1) err(warnprec,"nffactor[2]",prec);
! 942: k = k<<1; p1 = idealmullll(nf,p1,p1);
! 943: if (low_stack(lim, stack_lim(av1,2)))
! 944: {
! 945: if (DEBUGMEM>1) err(warnmem,"T2_matrix_pow");
! 946: p1 = gerepileupto(av1,p1);
! 947: }
! 948: if (!tot_real) T2 = nf_get_T2((GEN) nf[7], roots(x,prec));
! 949: p3 = qf_base_change(T2,p1,1);
! 950: }
! 951: }
! 952:
! 953: /* return the factorization of the square-free polynomial x.
! 954: The coeff of x are in Z_nf and its leading term is a rational
! 955: integer. If fl = 1,return only the roots of x in nf */
! 956: static GEN
! 957: nfsqff(GEN nf,GEN pol, long fl)
! 958: {
! 959: long d=lgef(pol),i,k,m,n,av=avma,tetpil,newprec,prec,nbf=BIGINT,anbf,ct=5;
! 960: GEN p1,pr,p2,rep,k2,C,h,dk,dki,p,prh,p3,T2,polbase,fact,pk,ap,apr;
! 961: GEN polmod,polred,hinv,lt,minp,den,maxp=shifti(gun,32),run,aprh;
! 962:
! 963: if (DEBUGLEVEL>=4) msgtimer("square-free");
! 964:
! 965: dk=absi((GEN)nf[3]);
! 966: dki=mulii(dk,(GEN)nf[4]);
! 967: n=degpol(nf[1]);
! 968:
! 969: polbase = unifpol(nf,pol,0);
! 970: polmod = unifpol(nf,pol,1);
! 971: dki=mulii(dki,gnorm((GEN)polmod[d-1]));
! 972:
! 973: prec = DEFAULTPREC;
! 974: for (;;)
! 975: {
! 976: if (prec <= gprecision(nf))
! 977: T2 = gprec_w(gmael(nf,5,3), prec);
! 978: else
! 979: T2 = nf_get_T2((GEN)nf[7], roots((GEN)nf[1], prec));
! 980:
! 981: run=realun(prec);
! 982: p1=realzero(prec);
! 983: for (i=2; i<d; i++)
! 984: {
! 985: p2 = gmul(run, (GEN)polbase[i]);
! 986: p2 = qfeval(T2, p2);
! 987: p1 = addrr(p1, gdiv(p2, binome(stoi(d), i-2)));
! 988: if (signe(p1) < 0) break;
! 989: }
! 990:
! 991: if (signe(p1) > 0) break;
! 992:
! 993: prec = (prec<<1)-2;
! 994: if (DEBUGLEVEL > 1) err(warnprec, "nfsqff", prec);
! 995: }
! 996:
! 997: lt = (GEN)leading_term(polbase)[1];
! 998: p1 = gmul(p1, mulis(sqri(lt), n));
! 999: C = gpow(stoi(3), gadd(gmulsg(3, ghalf), stoi(d)), prec);
! 1000: C = gdiv(gmul(C, p1), gmulsg(d, mppi(prec)));
! 1001:
! 1002: if (DEBUGLEVEL>=4)
! 1003: fprintferr("the bound on the T2-norm of the coeff. is: %Z\n", C);
! 1004:
! 1005: /* the theoretical bound for the exponent should be:
! 1006: k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.347))); */
! 1007: k2=gadd(glog(gdivgs(C,n),DEFAULTPREC), mulsr(n*(n-1), dbltor(0.15)));
! 1008: k2=gmul2n(gmulgs(k2,n),-1);
! 1009:
! 1010: minp=gmin(gexp(gmul2n(k2,-6),BIGDEFAULTPREC), maxp);
! 1011: minp=gceil(minp);
! 1012:
! 1013: if (DEBUGLEVEL>=4)
! 1014: {
! 1015: fprintferr("lower bound for the prime numbers: %Z\n", minp);
! 1016: msgtimer("bounds computation");
! 1017: }
! 1018:
! 1019: p = rep = polred = NULL; /* gcc -Wall */
! 1020: pr=NULL;
! 1021: for (;;)
! 1022: {
! 1023: apr = choose_prime(nf,dki,minp, pr?30:0);
! 1024: if (!apr) break;
! 1025:
! 1026: ap=(GEN)apr[1];
! 1027: aprh=prime_to_ideal(nf,apr);
! 1028:
! 1029: polred=gcopy(polbase);
! 1030: lt=(GEN)leading_term(polbase)[1];
! 1031: lt=mpinvmod(lt,ap);
! 1032:
! 1033: for (i=2; i<d; i++)
! 1034: polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), aprh)[1];
! 1035:
! 1036: if (!divise(discsr(polred),ap))
! 1037: {
! 1038: rep=(GEN)simplefactmod(polred,ap)[1];
! 1039: anbf=lg(rep)-1;
! 1040: ct--;
! 1041: if (anbf < nbf)
! 1042: {
! 1043: nbf=anbf;
! 1044: pr=gcopy(apr);
! 1045: p=gcopy(ap);
! 1046: if (DEBUGLEVEL>=4)
! 1047: {
! 1048: fprintferr("prime ideal considered: %Z\n", pr);
! 1049: fprintferr("number of irreducible factors: %ld\n", nbf);
! 1050: }
! 1051: if (nbf == 1) break;
! 1052: }
! 1053: }
! 1054: if (pr && !ct) break;
! 1055:
! 1056: minp = addis(ap,1);
! 1057: }
! 1058:
! 1059: k = itos(gceil(gdiv(k2,glog(p,BIGDEFAULTPREC))));
! 1060:
! 1061: if (DEBUGLEVEL>=4)
! 1062: {
! 1063: fprintferr("prime ideal chosen: %Z\n", pr);
! 1064: msgtimer("choice of the prime ideal");
! 1065: }
! 1066:
! 1067: if (lg(rep)==2)
! 1068: {
! 1069: if (fl) { avma=av; return cgetg(1,t_VEC); }
! 1070: rep=cgetg(2,t_VEC); rep[1]=lcopy(polmod);
! 1071: nfcmbf.nfact=1; return gerepileupto(av, rep);
! 1072: }
! 1073:
! 1074: p2=cgetr(DEFAULTPREC);
! 1075: affir(mulii(absi(dk),gpowgs(p,k)),p2);
! 1076: p2=shifti(gceil(mplog(p2)),-1);
! 1077:
! 1078: newprec = MEDDEFAULTPREC + (long)(itos(p2)*pariK1);
! 1079: if (DEBUGLEVEL>=4)
! 1080: fprintferr("new precision: %ld\n",newprec);
! 1081:
! 1082: prh = idealpows(nf,pr,k); m = k;
! 1083: h = T2_matrix_pow(nf,prh,p,C, &k, newprec);
! 1084: if (m != k) prh=idealpows(nf,pr,k);
! 1085:
! 1086: if (DEBUGLEVEL>=4)
! 1087: {
! 1088: fprintferr("a suitable exponent is: %ld\n",(long)k);
! 1089: msgtimer("computation of H");
! 1090: }
! 1091:
! 1092: pk = gcoeff(prh,1,1);
! 1093: lt=(GEN)leading_term(polbase)[1];
! 1094: lt=mpinvmod(lt,pk);
! 1095:
! 1096: polred[1] = polbase[1];
! 1097: for (i=2; i<d; i++)
! 1098: polred[i] = nfreducemodpr_i(gmul(lt,(GEN)polbase[i]), prh)[1];
! 1099:
! 1100: fact = lift_intern((GEN)factmod(polred,p)[1]);
! 1101: rep = hensel_lift_fact(polred,fact,NULL,p,pk,k);
! 1102:
! 1103: if (DEBUGLEVEL >= 4) msgtimer("computation of the p-adic factorization");
! 1104:
! 1105: den = det(h); /* |den| = NP^k */
! 1106: hinv= adj(h);
! 1107: lt=(GEN)leading_term(polbase)[1];
! 1108:
! 1109: if (fl)
! 1110: {
! 1111: long x_a[] = { evaltyp(t_POL)|m_evallg(4), 0,0,0 };
! 1112: GEN mlt = gneg_i(lt), rem;
! 1113: x_a[1] = polbase[1]; setlgef(x_a, 4);
! 1114: x_a[3] = un;
! 1115: p1=cgetg(lg(rep)+1,t_VEC);
! 1116: polbase = unifpol(nf,polbase,1);
! 1117: for (m=1,i=1; i<lg(rep); i++)
! 1118: {
! 1119: p2=(GEN)rep[i]; if (lgef(p2)!=4) break;
! 1120:
! 1121: p3 = algtobasis(nf, gmul(mlt,(GEN)p2[2]));
! 1122: p3 = nf_bestlift(h,hinv,den,p3);
! 1123: p3 = gdiv(p3,lt);
! 1124: x_a[2] = lneg(p3); /* check P(p3) == 0 */
! 1125: p2 = poldivres(polbase, unifpol(nf,x_a,1), &rem);
! 1126: if (!signe(rem)) { polbase = p2; p1[m++] = (long)p3; }
! 1127: }
! 1128: tetpil=avma; rep=cgetg(m,t_VEC);
! 1129: for (i=1; i<m; i++) rep[i]=(long)basistoalg(nf,(GEN)p1[i]);
! 1130: return gerepile(av,tetpil,rep);
! 1131: }
! 1132:
! 1133: for (i=1; i<lg(rep); i++)
! 1134: rep[i] = (long)unifpol(nf,(GEN)rep[i],0);
! 1135:
! 1136: nfcmbf.pol = polmod;
! 1137: nfcmbf.lt = lt;
! 1138: nfcmbf.h = h;
! 1139: nfcmbf.hinv = hinv;
! 1140: nfcmbf.den = den;
! 1141: nfcmbf.fact = rep;
! 1142: nfcmbf.res = cgetg(lg(rep)+1,t_VEC);
! 1143: nfcmbf.nfact = 0;
! 1144: nfcmbf.nfactmod = lg(rep)-1;
! 1145: nf_combine_factors(nf,1,NULL,d-3,1);
! 1146:
! 1147: if (DEBUGLEVEL >= 4) msgtimer("computation of the factors");
! 1148:
! 1149: i = nfcmbf.nfact;
! 1150: if (lgef(nfcmbf.pol)>3)
! 1151: {
! 1152: nfcmbf.res[++i] = (long) nf_pol_divres(nf,nfcmbf.pol,nfcmbf.lt,NULL);
! 1153: nfcmbf.nfact = i;
! 1154: }
! 1155:
! 1156: tetpil=avma; rep=cgetg(i+1,t_VEC);
! 1157: for ( ; i>=1; i--)
! 1158: rep[i]=(long)unifpol(nf,(GEN)nfcmbf.res[i],1);
! 1159: return gerepile(av,tetpil,rep);
! 1160: }
! 1161:
! 1162: static int
! 1163: nf_combine_factors(GEN nf,long fxn,GEN psf,long dlim,long hint)
! 1164: {
! 1165: int val = 0; /* assume failure */
! 1166: GEN newf, newpsf = NULL;
! 1167: long av,newd,ltop,i;
! 1168:
! 1169: /* Assertion: fxn <= nfcmbf.nfactmod && dlim > 0 */
! 1170:
! 1171: /* first, try deeper factors without considering the current one */
! 1172: if (fxn != nfcmbf.nfactmod)
! 1173: {
! 1174: val = nf_combine_factors(nf,fxn+1,psf,dlim,hint);
! 1175: if (val && psf) return 1;
! 1176: }
! 1177:
! 1178: /* second, try including the current modular factor in the product */
! 1179: newf = (GEN)nfcmbf.fact[fxn];
! 1180: if (!newf) return val; /* modular factor already used */
! 1181: newd = degpol(newf);
! 1182: if (newd > dlim) return val; /* degree of new factor is too large */
! 1183:
! 1184: av = avma;
! 1185: if (newd % hint == 0)
! 1186: {
! 1187: GEN p, quot,rem;
! 1188:
! 1189: newpsf = nf_pol_mul(nf, psf? psf: nfcmbf.lt, newf);
! 1190: newpsf = nf_pol_lift(nfcmbf.h,nfcmbf.hinv,nfcmbf.den,newpsf);
! 1191: /* try out the new combination */
! 1192: ltop=avma;
! 1193: quot=nf_pol_divres(nf,nfcmbf.pol,newpsf,&rem);
! 1194: if (gcmp0(rem)) /* found a factor */
! 1195: {
! 1196: p = nf_pol_mul(nf,element_inv(nf,leading_term(newpsf)),newpsf);
! 1197: nfcmbf.res[++nfcmbf.nfact] = (long) p; /* store factor */
! 1198: nfcmbf.fact[fxn]=0; /* remove used modular factor */
! 1199:
! 1200: /* fix up target */
! 1201: p=gun; quot=unifpol(nf,quot,0);
! 1202: for (i=2; i<lgef(quot); i++)
! 1203: p = mpppcm(p, denom((GEN)quot[i]));
! 1204:
! 1205: nfcmbf.pol = nf_pol_mul(nf,p,quot);
! 1206: nfcmbf.lt = leading_term(nfcmbf.pol);
! 1207: return 1;
! 1208: }
! 1209: avma=ltop;
! 1210: }
! 1211:
! 1212: /* If room in degree limit + more modular factors to try, add more factors to
! 1213: * newpsf */
! 1214: if (newd < dlim && fxn < nfcmbf.nfactmod
! 1215: && nf_combine_factors(nf,fxn+1,newpsf,dlim-newd,hint))
! 1216: {
! 1217: nfcmbf.fact[fxn]=0; /* remove used modular factor */
! 1218: return 1;
! 1219: }
! 1220: avma = av; return val;
! 1221: }
! 1222:
! 1223: /* return the characteristic polynomial of alpha over nf, where alpha
! 1224: is an element of the algebra nf[X]/(T) given as a polynomial in X */
! 1225: GEN
! 1226: rnfcharpoly(GEN nf,GEN T,GEN alpha,int v)
! 1227: {
! 1228: long av = avma, vnf, vT, lT;
! 1229: GEN p1;
! 1230:
! 1231: nf=checknf(nf); vnf = varn(nf[1]);
! 1232: if (v<0) v = 0;
! 1233: T = fix_relative_pol(nf,T,1);
! 1234: if (typ(alpha) == t_POLMOD) alpha = lift_to_pol(alpha);
! 1235: lT = lgef(T);
! 1236: if (typ(alpha) != t_POL || varn(alpha) == vnf)
! 1237: return gerepileupto(av, gpowgs(gsub(polx[v], alpha), lT - 3));
! 1238: vT = varn(T);
! 1239: if (varn(alpha) != vT || v >= vnf)
! 1240: err(talker,"incorrect variables in rnfcharpoly");
! 1241: if (lgef(alpha) >= lT) alpha = gmod(alpha,T);
! 1242: if (lT <= 4)
! 1243: return gerepileupto(av, gsub(polx[v], alpha));
! 1244: p1 = caract2(unifpol(nf,T,1), unifpol(nf,alpha,1), v);
! 1245: return gerepileupto(av, unifpol(nf,p1,1));
! 1246: }
! 1247:
! 1248: #if 0
! 1249: /* return the minimal polynomial of alpha over nf, where alpha is
! 1250: an element of the algebra nf[X]/(T) given as a polynomial in X */
! 1251: GEN
! 1252: rnfminpoly(GEN nf,GEN T,GEN alpha,int n)
! 1253: {
! 1254: long av=avma,tetpil;
! 1255: GEN p1,p2;
! 1256:
! 1257: nf=checknf(nf); p1=rnfcharpoly(nf,T,alpha,n);
! 1258: tetpil=avma; p2=nf_pol_subres(nf,p1,deriv(p1,varn(T)));
! 1259: if (lgef(p2)==3) { avma=tetpil; return p1; }
! 1260:
! 1261: p1 = nf_pol_divres(nf,p1,p2,NULL);
! 1262: p2 = element_inv(nf,leading_term(p1));
! 1263: tetpil=avma; return gerepile(av,tetpil,unifpol(nf,nf_pol_mul(nf,p2,p1),1));
! 1264: }
! 1265: #endif
! 1266:
! 1267: /* relative Dedekind criterion over nf, applied to the order defined by a
! 1268: * root of irreducible polynomial T, modulo the prime ideal pr. Returns
! 1269: * [flag,basis,val], where basis is a pseudo-basis of the enlarged order,
! 1270: * flag is 1 iff this order is pr-maximal, and val is the valuation in pr of
! 1271: * the order discriminant
! 1272: */
! 1273: GEN
! 1274: rnfdedekind(GEN nf,GEN T,GEN pr)
! 1275: {
! 1276: long av=avma,vt,r,d,da,n,m,i,j;
! 1277: GEN p1,p2,p,tau,g,vecun,veczero,matid;
! 1278: GEN prhall,res,h,k,base,Ca;
! 1279:
! 1280: nf=checknf(nf); Ca=unifpol(nf,T,0);
! 1281: res=cgetg(4,t_VEC);
! 1282: if (typ(pr)==t_VEC && lg(pr)==3)
! 1283: { prhall = (GEN)pr[2]; pr = (GEN)pr[1]; }
! 1284: else
! 1285: prhall=nfmodprinit(nf,pr);
! 1286: p=(GEN)pr[1]; tau=gdiv((GEN)pr[5],p);
! 1287: n=degpol(nf[1]); m=degpol(T);
! 1288:
! 1289: vecun = gscalcol_i(gun,n);
! 1290: veczero = zerocol(n);
! 1291:
! 1292: p1=(GEN)nffactormod(nf,Ca,pr)[1];
! 1293: r=lg(p1); if (r < 2) err(constpoler,"rnfdedekind");
! 1294: g=lift((GEN)p1[1]);
! 1295: for (i=2; i<r; i++)
! 1296: g = nf_pol_mul(nf,g, lift((GEN)p1[i]));
! 1297: h=nfmod_pol_divres(nf,prhall,Ca,g,NULL);
! 1298: k=nf_pol_mul(nf,tau,gsub(Ca, nf_pol_mul(nf,lift(g),lift(h))));
! 1299: p2=nfmod_pol_gcd(nf,prhall,g,h);
! 1300: k= nfmod_pol_gcd(nf,prhall,p2,k);
! 1301:
! 1302: d=degpol(k); /* <= m */
! 1303: vt = idealval(nf,discsr(T),pr) - 2*d;
! 1304: res[3]=lstoi(vt);
! 1305: if (!d || vt<=1) res[1]=un; else res[1]=zero;
! 1306:
! 1307: base=cgetg(3,t_VEC);
! 1308: p1=cgetg(m+d+1,t_MAT); base[1]=(long)p1;
! 1309: p2=cgetg(m+d+1,t_VEC); base[2]=(long)p2;
! 1310: /* if d > 0, base[2] temporarily multiplied by p, for the final nfhermitemod
! 1311: * [ which requires integral ideals ] */
! 1312: matid = gscalmat(d? p: gun, n);
! 1313: for (j=1; j<=m; j++)
! 1314: {
! 1315: p2[j]=(long)matid;
! 1316: p1[j]=lgetg(m+1,t_COL);
! 1317: for (i=1; i<=m; i++)
! 1318: coeff(p1,i,j) = (i==j)?(long)vecun:(long)veczero;
! 1319: }
! 1320: if (d)
! 1321: {
! 1322: GEN pal = lift(nfmod_pol_divres(nf,prhall,Ca,k,NULL));
! 1323: GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */
! 1324: GEN nfx = unifpol(nf,polx[varn(T)],0);
! 1325:
! 1326: for ( ; j<=m+d; j++)
! 1327: {
! 1328: p1[j]=lgetg(m+1,t_COL);
! 1329: da=degpol(pal);
! 1330: for (i=1; i<=da+1; i++) coeff(p1,i,j)=pal[i+1];
! 1331: for ( ; i<=m; i++) coeff(p1,i,j)=(long)veczero;
! 1332: p2[j]=(long)prinvp;
! 1333: nf_pol_divres(nf,nf_pol_mul(nf,pal,nfx),T,&pal);
! 1334: }
! 1335: /* the modulus is integral */
! 1336: base = nfhermitemod(nf,base, gmul(gpowgs(p, m-d),
! 1337: idealpows(nf, prinvp, d)));
! 1338: base[2] = ldiv((GEN)base[2], p); /* cancel the factor p */
! 1339: }
! 1340: res[2]=(long)base; return gerepilecopy(av, res);
! 1341: }
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