Annotation of OpenXM_contrib/pari-2.2/src/basemath/polarit1.c, Revision 1.1
1.1 ! noro 1: /* $Id: polarit1.c,v 1.66 2001/10/01 17:19:06 bill 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: /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
! 19: /** (first part) **/
! 20: /** **/
! 21: /***********************************************************************/
! 22: #include "pari.h"
! 23: extern GEN get_bas_den(GEN bas);
! 24: extern GEN get_mul_table(GEN x,GEN bas,GEN invbas,GEN *T);
! 25: extern GEN pol_to_monic(GEN pol, GEN *lead);
! 26:
! 27: /* see splitgen() for how to use these two */
! 28: GEN
! 29: setloop(GEN a)
! 30: {
! 31: a=icopy(a); new_chunk(2); /* dummy to get two cells of extra space */
! 32: return a;
! 33: }
! 34:
! 35: /* assume a > 0 */
! 36: GEN
! 37: incpos(GEN a)
! 38: {
! 39: long i,l=lgefint(a);
! 40:
! 41: for (i=l-1; i>1; i--)
! 42: if (++a[i]) return a;
! 43: i=l+1; a--; /* use extra cell */
! 44: a[0]=evaltyp(1) | evallg(i);
! 45: a[1]=evalsigne(1) | evallgefint(i);
! 46: return a;
! 47: }
! 48:
! 49: GEN
! 50: incloop(GEN a)
! 51: {
! 52: long i,l;
! 53:
! 54: switch(signe(a))
! 55: {
! 56: case 0:
! 57: a--; /* use extra cell */
! 58: a[0]=evaltyp(t_INT) | evallg(3);
! 59: a[1]=evalsigne(1) | evallgefint(3);
! 60: a[2]=1; return a;
! 61:
! 62: case -1:
! 63: l=lgefint(a);
! 64: for (i=l-1; i>1; i--)
! 65: if (a[i]--) break;
! 66: if (a[2] == 0)
! 67: {
! 68: a++; /* save one cell */
! 69: a[0] = evaltyp(t_INT) | evallg(2);
! 70: a[1] = evalsigne(0) | evallgefint(2);
! 71: }
! 72: return a;
! 73:
! 74: default:
! 75: return incpos(a);
! 76: }
! 77: }
! 78:
! 79: /*******************************************************************/
! 80: /* */
! 81: /* DIVISIBILITE */
! 82: /* Return 1 if y | x, 0 otherwise */
! 83: /* */
! 84: /*******************************************************************/
! 85:
! 86: int
! 87: gdivise(GEN x, GEN y)
! 88: {
! 89: long av=avma;
! 90: x=gmod(x,y); avma=av; return gcmp0(x);
! 91: }
! 92:
! 93: int
! 94: poldivis(GEN x, GEN y, GEN *z)
! 95: {
! 96: long av = avma;
! 97: GEN p1 = poldivres(x,y,ONLY_DIVIDES);
! 98: if (p1) { *z = p1; return 1; }
! 99: avma=av; return 0;
! 100: }
! 101:
! 102: /*******************************************************************/
! 103: /* */
! 104: /* POLYNOMIAL EUCLIDEAN DIVISION */
! 105: /* */
! 106: /*******************************************************************/
! 107: /* Polynomial division x / y:
! 108: * if z = ONLY_REM return remainder, otherwise return quotient
! 109: * if z != NULL set *z to remainder
! 110: * *z is the last object on stack (and thus can be disposed of with cgiv
! 111: * instead of gerepile)
! 112: */
! 113: GEN
! 114: poldivres(GEN x, GEN y, GEN *pr)
! 115: {
! 116: ulong avy,av,av1;
! 117: long ty=typ(y),tx,vx,vy,dx,dy,dz,i,j,sx,lrem;
! 118: int remainder;
! 119: GEN z,p1,rem,y_lead,mod;
! 120: GEN (*f)(GEN,GEN);
! 121:
! 122: if (pr == ONLY_DIVIDES_EXACT)
! 123: { f = gdivexact; pr = ONLY_DIVIDES; }
! 124: else
! 125: f = gdiv;
! 126: if (is_scalar_t(ty))
! 127: {
! 128: if (pr == ONLY_REM) return gzero;
! 129: if (pr && pr != ONLY_DIVIDES) *pr=gzero;
! 130: return f(x,y);
! 131: }
! 132: tx=typ(x); vy=gvar9(y);
! 133: if (is_scalar_t(tx) || gvar9(x)>vy)
! 134: {
! 135: if (pr == ONLY_REM) return gcopy(x);
! 136: if (pr == ONLY_DIVIDES) return gcmp0(x)? gzero: NULL;
! 137: if (pr) *pr=gcopy(x);
! 138: return gzero;
! 139: }
! 140: if (tx!=t_POL || ty!=t_POL) err(typeer,"euclidean division (poldivres)");
! 141:
! 142: vx=varn(x);
! 143: if (vx<vy)
! 144: {
! 145: if (pr && pr != ONLY_DIVIDES)
! 146: {
! 147: p1 = zeropol(vx); if (pr == ONLY_REM) return p1;
! 148: *pr = p1;
! 149: }
! 150: return f(x,y);
! 151: }
! 152: if (!signe(y)) err(talker,"euclidean division by zero (poldivres)");
! 153:
! 154: dy=degpol(y); y_lead = (GEN)y[dy+2];
! 155: if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */
! 156: {
! 157: err(warner,"normalizing a polynomial with 0 leading term");
! 158: for (dy--; dy>=0; dy--)
! 159: {
! 160: y_lead = (GEN)y[dy+2];
! 161: if (!gcmp0(y_lead)) break;
! 162: }
! 163: }
! 164: if (!dy) /* y is constant */
! 165: {
! 166: if (pr && pr != ONLY_DIVIDES)
! 167: {
! 168: if (pr == ONLY_REM) return zeropol(vx);
! 169: *pr = zeropol(vx);
! 170: }
! 171: return f(x, constant_term(y));
! 172: }
! 173: dx=degpol(x);
! 174: if (vx>vy || dx<dy)
! 175: {
! 176: if (pr)
! 177: {
! 178: if (pr == ONLY_DIVIDES) return gcmp0(x)? gzero: NULL;
! 179: if (pr == ONLY_REM) return gcopy(x);
! 180: *pr = gcopy(x);
! 181: }
! 182: return zeropol(vy);
! 183: }
! 184: dz=dx-dy; av=avma; /* to avoid gsub's later on */
! 185: p1 = new_chunk(dy+3);
! 186: for (i=2; i<dy+3; i++)
! 187: {
! 188: GEN p2 = (GEN)y[i];
! 189: p1[i] = isexactzero(p2)? 0: (long)gneg_i(p2);
! 190: }
! 191: y = p1;
! 192: switch(typ(y_lead))
! 193: {
! 194: case t_INTMOD:
! 195: case t_POLMOD: y_lead = ginv(y_lead);
! 196: f = gmul; mod = gmodulcp(gun, (GEN)y_lead[1]);
! 197: break;
! 198: default: if (gcmp1(y_lead)) y_lead = NULL;
! 199: mod = NULL;
! 200: }
! 201: avy=avma; z=cgetg(dz+3,t_POL);
! 202: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 203: x += 2; y += 2; z += 2;
! 204:
! 205: p1 = (GEN)x[dx]; remainder = (pr == ONLY_REM);
! 206: z[dz]=y_lead? (long)f(p1,y_lead): lcopy(p1);
! 207: for (i=dx-1; i>=dy; i--)
! 208: {
! 209: av1=avma; p1=(GEN)x[i];
! 210: for (j=i-dy+1; j<=i && j<=dz; j++)
! 211: if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 212: if (y_lead) p1 = f(p1,y_lead);
! 213: if (!remainder) p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
! 214: z[i-dy] = (long)p1;
! 215: }
! 216: if (!pr) return gerepileupto(av,z-2);
! 217:
! 218: rem = (GEN)avma; av1 = (long)new_chunk(dx+3);
! 219: for (sx=0; ; i--)
! 220: {
! 221: p1 = (GEN)x[i];
! 222: /* we always enter this loop at least once */
! 223: for (j=0; j<=i && j<=dz; j++)
! 224: if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 225: if (mod && avma==av1) p1 = gmul(p1,mod);
! 226: if (!gcmp0(p1)) { sx = 1; break; } /* remainder is non-zero */
! 227: if (!isinexactreal(p1) && !isexactzero(p1)) break;
! 228: if (!i) break;
! 229: avma=av1;
! 230: }
! 231: if (pr == ONLY_DIVIDES)
! 232: {
! 233: if (sx) { avma=av; return NULL; }
! 234: avma = (long)rem;
! 235: return gerepileupto(av,z-2);
! 236: }
! 237: lrem=i+3; rem -= lrem;
! 238: if (avma==av1) { avma = (long)rem; p1 = gcopy(p1); }
! 239: else p1 = gerepileupto((long)rem,p1);
! 240: rem[0]=evaltyp(t_POL) | evallg(lrem);
! 241: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
! 242: rem += 2;
! 243: rem[i]=(long)p1;
! 244: for (i--; i>=0; i--)
! 245: {
! 246: av1=avma; p1 = (GEN)x[i];
! 247: for (j=0; j<=i && j<=dz; j++)
! 248: if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 249: if (mod && avma==av1) p1 = gmul(p1,mod);
! 250: rem[i]=avma==av1? lcopy(p1):lpileupto(av1,p1);
! 251: }
! 252: rem -= 2;
! 253: if (!sx) normalizepol_i(rem, lrem);
! 254: if (remainder) return gerepileupto(av,rem);
! 255: z -= 2;
! 256: {
! 257: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
! 258: gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
! 259: }
! 260: }
! 261:
! 262: /*******************************************************************/
! 263: /* */
! 264: /* ROOTS MODULO a prime p (no multiplicities) */
! 265: /* */
! 266: /*******************************************************************/
! 267: static GEN
! 268: mod(GEN x, GEN y)
! 269: {
! 270: GEN z = cgetg(3,t_INTMOD);
! 271: z[1]=(long)y; z[2]=(long)x; return z;
! 272: }
! 273:
! 274: static long
! 275: factmod_init(GEN *F, GEN pp, long *p)
! 276: {
! 277: GEN f = *F;
! 278: long i,d;
! 279: if (typ(f)!=t_POL || typ(pp)!=t_INT) err(typeer,"factmod");
! 280: if (expi(pp) > BITS_IN_LONG - 3) *p = 0;
! 281: else
! 282: {
! 283: *p = itos(pp);
! 284: if (*p < 2) err(talker,"not a prime in factmod");
! 285: }
! 286: f = gmul(f, mod(gun,pp));
! 287: if (!signe(f)) err(zeropoler,"factmod");
! 288: f = lift_intern(f); d = lgef(f);
! 289: for (i=2; i <d; i++)
! 290: if (typ(f[i])!=t_INT) err(impl,"factormod for general polynomials");
! 291: *F = f; return d-3;
! 292: }
! 293:
! 294: #define mods(x,y) mod(stoi(x),y)
! 295: static GEN
! 296: root_mod_2(GEN f)
! 297: {
! 298: int z1, z0 = !signe(constant_term(f));
! 299: long i,n;
! 300: GEN y;
! 301:
! 302: for (i=2, n=1; i < lgef(f); i++)
! 303: if (signe(f[i])) n++;
! 304: z1 = n & 1;
! 305: y = cgetg(z0+z1+1, t_COL); i = 1;
! 306: if (z0) y[i++] = (long)mods(0,gdeux);
! 307: if (z1) y[i] = (long)mods(1,gdeux);
! 308: return y;
! 309: }
! 310:
! 311: #define i_mod4(x) (signe(x)? mod4((GEN)(x)): 0)
! 312: static GEN
! 313: root_mod_4(GEN f)
! 314: {
! 315: long no,ne;
! 316: int z0 = !signe(constant_term(f));
! 317: int z2 = ((i_mod4(constant_term(f)) + 2*i_mod4(f[3])) & 3) == 0;
! 318: int i,z1,z3;
! 319: GEN y,p;
! 320:
! 321: for (ne=0,i=2; i<lgef(f); i+=2)
! 322: if (signe(f[i])) ne += mael(f,i,2);
! 323: for (no=0,i=3; i<lgef(f); i+=2)
! 324: if (signe(f[i])) no += mael(f,i,2);
! 325: no &= 3; ne &= 3;
! 326: z3 = (no == ne);
! 327: z1 = (no == ((4-ne)&3));
! 328: y=cgetg(1+z0+z1+z2+z3,t_COL); i = 1; p = stoi(4);
! 329: if (z0) y[i++] = (long)mods(0,p);
! 330: if (z1) y[i++] = (long)mods(1,p);
! 331: if (z2) y[i++] = (long)mods(2,p);
! 332: if (z3) y[i] = (long)mods(3,p);
! 333: return y;
! 334: }
! 335: #undef i_mod4
! 336:
! 337: /* p even, accept p = 4 for p-adic stuff */
! 338: static GEN
! 339: root_mod_even(GEN f, long p)
! 340: {
! 341: switch(p)
! 342: {
! 343: case 2: return root_mod_2(f);
! 344: case 4: return root_mod_4(f);
! 345: }
! 346: err(talker,"not a prime in rootmod");
! 347: return NULL; /* not reached */
! 348: }
! 349:
! 350: /* by checking f(0..p-1) */
! 351: GEN
! 352: rootmod2(GEN f, GEN pp)
! 353: {
! 354: GEN g,y,ss,q,r, x_minus_s;
! 355: long p,av = avma,av1,d,i,nbrac;
! 356:
! 357: if (!(d = factmod_init(&f, pp, &p))) { avma=av; return cgetg(1,t_COL); }
! 358: if (!p) err(talker,"prime too big in rootmod2");
! 359: if ((p & 1) == 0) { avma = av; return root_mod_even(f,p); }
! 360: x_minus_s = gadd(polx[varn(f)], stoi(-1));
! 361:
! 362: nbrac=1;
! 363: y=(GEN)gpmalloc((d+1)*sizeof(long));
! 364: if (gcmp0(constant_term(f))) y[nbrac++] = 0;
! 365: ss = icopy(gun); av1 = avma;
! 366: do
! 367: {
! 368: mael(x_minus_s,2,2) = ss[2];
! 369: /* one might do a FFT-type evaluation */
! 370: q = FpX_divres(f, x_minus_s, pp, &r);
! 371: if (signe(r)) avma = av1;
! 372: else
! 373: {
! 374: y[nbrac++] = ss[2]; f = q; av1 = avma;
! 375: }
! 376: ss[2]++;
! 377: }
! 378: while (nbrac<d && p>ss[2]);
! 379: if (nbrac == 1) { avma=av; return cgetg(1,t_COL); }
! 380: if (nbrac == d && p != ss[2])
! 381: {
! 382: g = mpinvmod((GEN)f[3], pp); setsigne(g,-1);
! 383: ss = modis(mulii(g, (GEN)f[2]), p);
! 384: y[nbrac++]=ss[2];
! 385: }
! 386: avma=av; g=cgetg(nbrac,t_COL);
! 387: if (isonstack(pp)) pp = icopy(pp);
! 388: for (i=1; i<nbrac; i++) g[i]=(long)mods(y[i],pp);
! 389: free(y); return g;
! 390: }
! 391:
! 392: /* by splitting */
! 393: GEN
! 394: rootmod(GEN f, GEN p)
! 395: {
! 396: long av = avma,tetpil,n,i,j,la,lb;
! 397: GEN y,pol,a,b,q,pol0;
! 398:
! 399: if (!factmod_init(&f, p, &i)) { avma=av; return cgetg(1,t_COL); }
! 400: i = p[lgefint(p)-1];
! 401: if ((i & 1) == 0) { avma = av; return root_mod_even(f,i); }
! 402: i=2; while (!signe(f[i])) i++;
! 403: if (i == 2) j = 1;
! 404: else
! 405: {
! 406: j = lgef(f) - (i-2);
! 407: if (j==3) /* f = x^n */
! 408: {
! 409: avma = av; y = cgetg(2,t_COL);
! 410: y[1] = (long)gmodulsg(0,p);
! 411: return y;
! 412: }
! 413: a = cgetg(j, t_POL); /* a = f / x^{v_x(f)} */
! 414: a[1] = evalsigne(1) | evalvarn(varn(f)) | evallgef(j);
! 415: f += i-2; for (i=2; i<j; i++) a[i]=f[i];
! 416: j = 2; f = a;
! 417: }
! 418: q = shifti(p,-1);
! 419: /* take gcd(x^(p-1) - 1, f) by splitting (x^q-1) * (x^q+1) */
! 420: b = FpXQ_pow(polx[varn(f)],q, f,p);
! 421: if (lgef(b)<3) err(talker,"not a prime in rootmod");
! 422: b = ZX_s_add(b,-1); /* b = x^((p-1)/2) - 1 mod f */
! 423: a = FpX_gcd(f,b, p);
! 424: b = ZX_s_add(b, 2); /* b = x^((p-1)/2) + 1 mod f */
! 425: b = FpX_gcd(f,b, p);
! 426: la = degpol(a);
! 427: lb = degpol(b); n = la + lb;
! 428: if (!n)
! 429: {
! 430: avma = av; y = cgetg(n+j,t_COL);
! 431: if (j>1) y[1] = (long)gmodulsg(0,p);
! 432: return y;
! 433: }
! 434: y = cgetg(n+j,t_COL);
! 435: if (j>1) { y[1] = zero; n++; }
! 436: y[j] = (long)FpX_normalize(b,p);
! 437: if (la) y[j+lb] = (long)FpX_normalize(a,p);
! 438: pol = gadd(polx[varn(f)], gun); pol0 = constant_term(pol);
! 439: while (j<=n)
! 440: {
! 441: a=(GEN)y[j]; la=degpol(a);
! 442: if (la==1)
! 443: y[j++] = lsubii(p, constant_term(a));
! 444: else if (la==2)
! 445: {
! 446: GEN d = subii(sqri((GEN)a[3]), shifti((GEN)a[2],2));
! 447: GEN e = mpsqrtmod(d,p), u = addis(q, 1); /* u = 1/2 */
! 448: y[j++] = lmodii(mulii(u,subii(e,(GEN)a[3])), p);
! 449: y[j++] = lmodii(mulii(u,negi(addii(e,(GEN)a[3]))), p);
! 450: }
! 451: else for (pol0[2]=1; ; pol0[2]++)
! 452: {
! 453: b = ZX_s_add(FpXQ_pow(pol,q, a,p), -1); /* pol^(p-1)/2 - 1 */
! 454: b = FpX_gcd(a,b, p); lb = degpol(b);
! 455: if (lb && lb<la)
! 456: {
! 457: b = FpX_normalize(b, p);
! 458: y[j+lb] = (long)FpX_div(a,b, p);
! 459: y[j] = (long)b; break;
! 460: }
! 461: }
! 462: }
! 463: tetpil = avma; y = gerepile(av,tetpil,sort(y));
! 464: if (isonstack(p)) p = icopy(p);
! 465: for (i=1; i<=n; i++) y[i] = (long)mod((GEN)y[i], p);
! 466: return y;
! 467: }
! 468:
! 469: GEN
! 470: rootmod0(GEN f, GEN p, long flag)
! 471: {
! 472: switch(flag)
! 473: {
! 474: case 0: return rootmod(f,p);
! 475: case 1: return rootmod2(f,p);
! 476: default: err(flagerr,"polrootsmod");
! 477: }
! 478: return NULL; /* not reached */
! 479: }
! 480:
! 481: /*******************************************************************/
! 482: /* */
! 483: /* FACTORISATION MODULO p */
! 484: /* */
! 485: /*******************************************************************/
! 486: static GEN spec_FpXQ_pow(GEN x, GEN p, GEN S);
! 487: /* Functions giving information on the factorisation. */
! 488:
! 489: /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
! 490: static GEN
! 491: Berlekamp_ker(GEN u, GEN p)
! 492: {
! 493: long i,j,d,N = degpol(u);
! 494: GEN vker,v,w,Q,p1,p2;
! 495: if (DEBUGLEVEL > 7) timer2();
! 496: Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);
! 497: w = v = FpXQ_pow(polx[varn(u)],p,u,p);
! 498: for (j=2; j<=N; j++)
! 499: {
! 500: Q[j] = lgetg(N+1,t_COL); p1 = (GEN)Q[j];
! 501: d = lgef(w)-1; p2 = w+1;
! 502: for (i=1; i<d ; i++) p1[i] = p2[i];
! 503: for ( ; i<=N; i++) p1[i] = zero;
! 504: p1[j] = laddis((GEN)p1[j], -1);
! 505: if (j < N)
! 506: {
! 507: ulong av = avma;
! 508: w = gerepileupto(av, FpX_res(gmul(w,v), u, p));
! 509: }
! 510: }
! 511: if (DEBUGLEVEL > 7) msgtimer("frobenius");
! 512: vker = FpM_ker(Q,p);
! 513: if (DEBUGLEVEL > 7) msgtimer("kernel");
! 514: return vker;
! 515: }
! 516:
! 517: /* f in ZZ[X] and p a prime number. */
! 518: long
! 519: FpX_is_squarefree(GEN f, GEN p)
! 520: {
! 521: long av = avma;
! 522: GEN z;
! 523: z = FpX_gcd(f,derivpol(f),p);
! 524: avma = av;
! 525: return lgef(z)==3;
! 526: }
! 527: /* idem
! 528: * leading term of f must be prime to p.
! 529: */
! 530: /* Compute the number of roots in Fp without counting multiplicity
! 531: * return -1 for 0 polynomial.
! 532: */
! 533: long
! 534: FpX_nbroots(GEN f, GEN p)
! 535: {
! 536: long av = avma, n=lgef(f);
! 537: GEN z;
! 538: if (n <= 4) return n-3;
! 539: f = FpX_red(f, p);
! 540: z = FpXQ_pow(polx[varn(f)], p, f, p);
! 541: z = FpX_sub(z,polx[varn(f)],NULL);
! 542: z = FpX_gcd(z,f,p),
! 543: avma = av; return degpol(z);
! 544: }
! 545: long
! 546: FpX_is_totally_split(GEN f, GEN p)
! 547: {
! 548: long av = avma, n=lgef(f);
! 549: GEN z;
! 550: if (n <= 4) return 1;
! 551: if (!is_bigint(p) && n-3 > p[2]) return 0;
! 552: f = FpX_red(f, p);
! 553: z = FpXQ_pow(polx[varn(f)], p, f, p);
! 554: avma = av; return lgef(z)==4 && gcmp1((GEN)z[3]) && !signe(z[2]);
! 555: }
! 556: /* u in ZZ[X] and p a prime number.
! 557: * u must be squarefree mod p.
! 558: * leading term of u must be prime to p. */
! 559: long
! 560: FpX_nbfact(GEN u, GEN p)
! 561: {
! 562: ulong av = avma;
! 563: GEN vker = Berlekamp_ker(u,p);
! 564: avma = av; return lg(vker)-1;
! 565: }
! 566:
! 567: /* Please use only use this function when you it is false, or that there is a
! 568: * lot of factors. If you believe f is irreducible or that it has few factors,
! 569: * then use `FpX_nbfact(f,p)==1' instead (faster).
! 570: */
! 571: static GEN factcantor0(GEN f, GEN pp, long flag);
! 572: long FpX_is_irred(GEN f, GEN p) { return !!factcantor0(f,p,2); }
! 573:
! 574: static GEN modulo;
! 575: static GEN gsmul(GEN a,GEN b){return FpX_mul(a,b,modulo);}
! 576: GEN
! 577: FpV_roots_to_pol(GEN V, GEN p, long v)
! 578: {
! 579: ulong ltop=avma;
! 580: long i;
! 581: GEN g=cgetg(lg(V),t_VEC);
! 582: for(i=1;i<lg(V);i++)
! 583: g[i]=(long)deg1pol(gun,negi((GEN)V[i]),v);
! 584: modulo=p;
! 585: g=divide_conquer_prod(g,&gsmul);
! 586: return gerepileupto(ltop,g);
! 587: }
! 588:
! 589: /************************************************************/
! 590: GEN
! 591: trivfact(void)
! 592: {
! 593: GEN y=cgetg(3,t_MAT);
! 594: y[1]=lgetg(1,t_COL);
! 595: y[2]=lgetg(1,t_COL); return y;
! 596: }
! 597:
! 598: static void
! 599: fqunclone(GEN x, GEN a, GEN p)
! 600: {
! 601: long i,j,lx = lgef(x);
! 602: for (i=2; i<lx; i++)
! 603: {
! 604: GEN p1 = (GEN)x[i];
! 605: if (typ(p1) == t_POLMOD) { p1[1] = (long)a; p1 = (GEN)p1[2]; }
! 606: if (typ(p1) == t_INTMOD) p1[1] = (long)p;
! 607: else /* t_POL */
! 608: for (j = lgef(p1)-1; j > 1; j--)
! 609: {
! 610: GEN p2 = (GEN)p1[j];
! 611: if (typ(p2) == t_INTMOD) p2[1] = (long)p;
! 612: }
! 613: }
! 614: }
! 615:
! 616: static GEN
! 617: try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)
! 618: {
! 619: GEN w2, w = FpXQ_pow(w0,q, pol,p);
! 620: long s;
! 621: if (gcmp1(w)) return w0;
! 622: for (s=1; s<r; s++,w=w2)
! 623: {
! 624: w2 = gsqr(w);
! 625: w2 = FpX_res(w2, pol, p);
! 626: if (gcmp1(w2)) break;
! 627: }
! 628: return gcmp_1(w)? NULL: w;
! 629: }
! 630:
! 631: /* INPUT:
! 632: * m integer (converted to polynomial w in Z[X] by stopoly)
! 633: * p prime; q = (p^d-1) / 2^r
! 634: * t[0] polynomial of degree k*d product of k irreducible factors of degree d
! 635: * t[0] is expected to be normalized (leading coeff = 1)
! 636: * OUTPUT:
! 637: * t[0],t[1]...t[k-1] the k factors, normalized
! 638: */
! 639: static void
! 640: split(long m, GEN *t, long d, GEN p, GEN q, long r, GEN S)
! 641: {
! 642: long ps,l,v,dv,av0,av;
! 643: GEN w,w0;
! 644:
! 645: dv=degpol(*t); if (dv==d) return;
! 646: v=varn(*t); av0=avma; ps = p[2];
! 647: for(av=avma;;avma=av)
! 648: {
! 649: if (ps==2)
! 650: {
! 651: w0=w=gpuigs(polx[v],m-1); m+=2;
! 652: for (l=1; l<d; l++)
! 653: w = gadd(w0, spec_FpXQ_pow(w, p, S));
! 654: }
! 655: else
! 656: {
! 657: w = FpX_res(stopoly(m,ps,v),*t, p);
! 658: m++; w = try_pow(w,*t,p,q,r);
! 659: if (!w) continue;
! 660: w = ZX_s_add(w, -1);
! 661: }
! 662: w = FpX_gcd(*t,w, p);
! 663: l = degpol(w); if (l && l!=dv) break;
! 664: }
! 665: w = FpX_normalize(w, p);
! 666: w = gerepileupto(av0, w);
! 667: l /= d; t[l]=FpX_div(*t,w,p); *t=w;
! 668: split(m,t+l,d,p,q,r,S);
! 669: split(m,t, d,p,q,r,S);
! 670: }
! 671:
! 672: static void
! 673: splitgen(GEN m, GEN *t, long d, GEN p, GEN q, long r)
! 674: {
! 675: long l,v,dv,av;
! 676: GEN w;
! 677:
! 678: dv=degpol(*t); if (dv==d) return;
! 679: v=varn(*t); m=setloop(m); m=incpos(m);
! 680: av=avma;
! 681: for(;; avma=av, m=incpos(m))
! 682: {
! 683: w = FpX_res(stopoly_gen(m,p,v),*t, p);
! 684: w = try_pow(w,*t,p,q,r);
! 685: if (!w) continue;
! 686: w = ZX_s_add(w,-1);
! 687: w = FpX_gcd(*t,w, p); l=degpol(w);
! 688: if (l && l!=dv) break;
! 689:
! 690: }
! 691: w = FpX_normalize(w, p);
! 692: w = gerepileupto(av, w);
! 693: l /= d; t[l]=FpX_div(*t,w,p); *t=w;
! 694: splitgen(m,t+l,d,p,q,r);
! 695: splitgen(m,t, d,p,q,r);
! 696: }
! 697:
! 698: /* return S = [ x^p, x^2p, ... x^(n-1)p ] mod (p, T), n = degree(T) > 0 */
! 699: static GEN
! 700: init_pow_p_mod_pT(GEN p, GEN T)
! 701: {
! 702: long i, n = degpol(T), v = varn(T);
! 703: GEN p1, S = cgetg(n, t_VEC);
! 704: if (n == 1) return S;
! 705: S[1] = (long)FpXQ_pow(polx[v], p, T, p);
! 706: /* use as many squarings as possible */
! 707: for (i=2; i < n; i+=2)
! 708: {
! 709: p1 = gsqr((GEN)S[i>>1]);
! 710: S[i] = (long)FpX_res(p1, T, p);
! 711: if (i == n-1) break;
! 712: p1 = gmul((GEN)S[i], (GEN)S[1]);
! 713: S[i+1] = (long)FpX_res(p1, T, p);
! 714: }
! 715: return S;
! 716: }
! 717:
! 718: /* compute x^p, S is as above */
! 719: static GEN
! 720: spec_FpXQ_pow(GEN x, GEN p, GEN S)
! 721: {
! 722: long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);
! 723: GEN x0 = x+2, z;
! 724: z = (GEN)x0[0];
! 725: if (dx < 0) err(talker, "zero polynomial in FpXQ_pow. %Z not prime", p);
! 726: for (i = 1; i <= dx; i++)
! 727: {
! 728: GEN d, c = (GEN)x0[i]; /* assume coeffs in [0, p-1] */
! 729: if (!signe(c)) continue;
! 730: d = (GEN)S[i]; if (!gcmp1(c)) d = gmul(c,d);
! 731: z = gadd(z, d);
! 732: if (low_stack(lim, stack_lim(av,1)))
! 733: {
! 734: if(DEBUGMEM>1) err(warnmem,"spec_FpXQ_pow");
! 735: z = gerepileupto(av, z);
! 736: }
! 737: }
! 738: z = FpX_red(z, p);
! 739: return gerepileupto(av, z);
! 740: }
! 741:
! 742: /* factor f mod pp.
! 743: * If (flag = 1) return the degrees, not the factors
! 744: * If (flag = 2) return NULL if f is not irreducible
! 745: */
! 746: static GEN
! 747: factcantor0(GEN f, GEN pp, long flag)
! 748: {
! 749: long i,j,k,d,e,vf,p,nbfact,tetpil,av = avma;
! 750: GEN ex,y,f2,g,g1,u,v,pd,q;
! 751: GEN *t;
! 752:
! 753: if (!(d = factmod_init(&f, pp, &p))) { avma=av; return trivfact(); }
! 754: /* to hold factors and exponents */
! 755: t = (GEN*)cgetg(d+1,t_VEC); ex = new_chunk(d+1);
! 756: vf=varn(f); e = nbfact = 1;
! 757: for(;;)
! 758: {
! 759: f2 = FpX_gcd(f,derivpol(f), pp);
! 760: if (flag > 1 && lgef(f2) > 3) return NULL;
! 761: g1 = FpX_div(f,f2,pp);
! 762: k = 0;
! 763: while (lgef(g1)>3)
! 764: {
! 765: long du,dg;
! 766: GEN S;
! 767: k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }
! 768: u = g1; g1 = FpX_gcd(f2,g1, pp);
! 769: if (lgef(g1)>3)
! 770: {
! 771: u = FpX_div( u,g1,pp);
! 772: f2= FpX_div(f2,g1,pp);
! 773: }
! 774: du = degpol(u);
! 775: if (du <= 0) continue;
! 776:
! 777: /* here u is square-free (product of irred. of multiplicity e * k) */
! 778: S = init_pow_p_mod_pT(pp, u);
! 779: pd=gun; v=polx[vf];
! 780: for (d=1; d <= du>>1; d++)
! 781: {
! 782: if (!flag) pd = mulii(pd,pp);
! 783: v = spec_FpXQ_pow(v, pp, S);
! 784: g = FpX_gcd(gadd(v, gneg(polx[vf])), u, pp);
! 785: dg = degpol(g);
! 786: if (dg <= 0) continue;
! 787:
! 788: /* Ici g est produit de pol irred ayant tous le meme degre d; */
! 789: j=nbfact+dg/d;
! 790:
! 791: if (flag)
! 792: {
! 793: if (flag > 1) return NULL;
! 794: for ( ; nbfact<j; nbfact++) { t[nbfact]=(GEN)d; ex[nbfact]=e*k; }
! 795: }
! 796: else
! 797: {
! 798: long r;
! 799: g = FpX_normalize(g, pp);
! 800: t[nbfact]=g; q = subis(pd,1); /* also ok for p=2: unused */
! 801: r = vali(q); q = shifti(q,-r);
! 802: /* le premier parametre est un entier variable m qui sera
! 803: * converti en un polynome w dont les coeff sont ses digits en
! 804: * base p (initialement m = p --> X) pour faire pgcd de g avec
! 805: * w^(p^d-1)/2 jusqu'a casser. p = 2 is treated separately.
! 806: */
! 807: if (p)
! 808: split(p,t+nbfact,d,pp,q,r,S);
! 809: else
! 810: splitgen(pp,t+nbfact,d,pp,q,r);
! 811: for (; nbfact<j; nbfact++) ex[nbfact]=e*k;
! 812: }
! 813: du -= dg;
! 814: u = FpX_div(u,g,pp);
! 815: v = FpX_res(v,u,pp);
! 816: }
! 817: if (du)
! 818: {
! 819: t[nbfact] = flag? (GEN)du: FpX_normalize(u, pp);
! 820: ex[nbfact++]=e*k;
! 821: }
! 822: }
! 823: j = lgef(f2); if (j==3) break;
! 824:
! 825: e*=p; j=(j-3)/p+3; setlg(f,j); setlgef(f,j);
! 826: for (i=2; i<j; i++) f[i]=f2[p*(i-2)+2];
! 827: }
! 828: if (flag > 1) { avma = av; return gun; } /* irreducible */
! 829: tetpil=avma; y=cgetg(3,t_MAT);
! 830: if (!flag)
! 831: {
! 832: y[1]=(long)t; setlg(t, nbfact);
! 833: y[2]=(long)ex; (void)sort_factor(y,cmpii);
! 834: }
! 835: u=cgetg(nbfact,t_COL); y[1]=(long)u;
! 836: v=cgetg(nbfact,t_COL); y[2]=(long)v;
! 837: if (flag)
! 838: for (j=1; j<nbfact; j++)
! 839: {
! 840: u[j] = lstoi((long)t[j]);
! 841: v[j] = lstoi(ex[j]);
! 842: }
! 843: else
! 844: for (j=1; j<nbfact; j++)
! 845: {
! 846: u[j] = (long)FpX(t[j], pp);
! 847: v[j] = lstoi(ex[j]);
! 848: }
! 849: return gerepile(av,tetpil,y);
! 850: }
! 851:
! 852: GEN
! 853: factcantor(GEN f, GEN p)
! 854: {
! 855: return factcantor0(f,p,0);
! 856: }
! 857:
! 858: GEN
! 859: simplefactmod(GEN f, GEN p)
! 860: {
! 861: return factcantor0(f,p,1);
! 862: }
! 863:
! 864: /* vector of polynomials (in v) whose coeffs are given by the columns of x */
! 865: GEN
! 866: mat_to_vecpol(GEN x, long v)
! 867: {
! 868: long i,j, lx = lg(x), lcol = lg(x[1]);
! 869: GEN y = cgetg(lx, t_VEC);
! 870:
! 871: for (j=1; j<lx; j++)
! 872: {
! 873: GEN p1, col = (GEN)x[j];
! 874: long k = lcol;
! 875:
! 876: while (k-- && gcmp0((GEN)col[k]));
! 877: i=k+2; p1=cgetg(i,t_POL);
! 878: p1[1] = evalsigne(1) | evallgef(i) | evalvarn(v);
! 879: col--; for (k=2; k<i; k++) p1[k] = col[k];
! 880: y[j] = (long)p1;
! 881: }
! 882: return y;
! 883: }
! 884:
! 885: /* matrix whose entries are given by the coeffs of the polynomials in
! 886: * vector v (considered as degree n-1 polynomials) */
! 887: GEN
! 888: vecpol_to_mat(GEN v, long n)
! 889: {
! 890: long i,j,d,N = lg(v);
! 891: GEN p1,w, y = cgetg(N, t_MAT);
! 892: if (typ(v) != t_VEC) err(typeer,"vecpol_to_mat");
! 893: n++;
! 894: for (j=1; j<N; j++)
! 895: {
! 896: p1 = cgetg(n,t_COL); y[j] = (long)p1;
! 897: w = (GEN)v[j];
! 898: if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }
! 899: else
! 900: {
! 901: d=lgef(w)-1; w++;
! 902: for (i=1; i<d; i++) p1[i] = w[i];
! 903: }
! 904: for ( ; i<n; i++) p1[i] = zero;
! 905: }
! 906: return y;
! 907: }
! 908: /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
! 909: GEN
! 910: mat_to_polpol(GEN x, long v,long w)
! 911: {
! 912: long i,j, lx = lg(x), lcol = lg(x[1]);
! 913: GEN y = cgetg(lx+1, t_POL);
! 914: y[1]=evalsigne(1) | evallgef(lx+1) | evalvarn(v);
! 915: y++;
! 916: for (j=1; j<lx; j++)
! 917: {
! 918: GEN p1, col = (GEN)x[j];
! 919: long k;
! 920: i=lcol+1; p1=cgetg(i,t_POL);
! 921: p1[1] = evalsigne(1) | evallgef(i) | evalvarn(w);
! 922: col--; for (k=2; k<i; k++) p1[k] = col[k];
! 923: y[j] = (long)normalizepol_i(p1,i);
! 924: }
! 925: return normalizepol_i(--y,lx+1);
! 926: }
! 927:
! 928: /* matrix whose entries are given by the coeffs of the polynomial v in
! 929: * two variables (considered as degree n polynomials) */
! 930: GEN
! 931: polpol_to_mat(GEN v, long n)
! 932: {
! 933: long i,j,d,N = lgef(v)-1;
! 934: GEN p1,w, y = cgetg(N, t_MAT);
! 935: if (typ(v) != t_POL) err(typeer,"polpol_to_mat");
! 936: n++;v++;
! 937: for (j=1; j<N; j++)
! 938: {
! 939: p1 = cgetg(n,t_COL); y[j] = (long)p1;
! 940: w = (GEN)v[j];
! 941: if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }
! 942: else
! 943: {
! 944: d=lgef(w)-1; w++;
! 945: for (i=1; i<d; i++) p1[i] = w[i];
! 946: }
! 947: for ( ; i<n; i++) p1[i] = zero;
! 948: }
! 949: return y;
! 950: }
! 951:
! 952:
! 953: /* set x <-- x + c*y mod p */
! 954: static void
! 955: split_berlekamp_addmul(GEN x, GEN y, long c, long p)
! 956: {
! 957: long i,lx,ly,l;
! 958: if (!c) return;
! 959: lx = lgef(x); ly = lgef(y); l = min(lx,ly);
! 960: if (p & ~MAXHALFULONG)
! 961: {
! 962: for (i=2; i<l; i++) x[i] = ((ulong)x[i]+ (ulong)mulssmod(c,y[i],p)) % p;
! 963: for ( ; i<ly; i++) x[i] = mulssmod(c,y[i],p);
! 964: }
! 965: else
! 966: {
! 967: for (i=2; i<l; i++) x[i] = ((ulong)x[i] + (ulong)(c*y[i])) % p;
! 968: for ( ; i<ly; i++) x[i] = (c*y[i]) % p;
! 969: }
! 970: do i--; while (i>1 && !x[i]);
! 971: if (i==1) setsigne(x,0); else { setsigne(x,1); setlgef(x,i+1); }
! 972: }
! 973:
! 974: long
! 975: split_berlekamp(GEN *t, GEN pp, GEN pps2)
! 976: {
! 977: GEN u = *t, p1, p2, vker,pol;
! 978: long av,N = degpol(u), d,i,kk,l1,l2,p, vu = varn(u);
! 979: ulong av0 = avma;
! 980:
! 981: vker = Berlekamp_ker(u,pp);
! 982: vker = mat_to_vecpol(vker,vu);
! 983: d = lg(vker)-1;
! 984: p = is_bigint(pp)? 0: pp[2];
! 985: if (p)
! 986: {
! 987: avma = av0; p1 = cgetg(d+1, t_VEC); /* hack: hidden gerepile */
! 988: for (i=1; i<=d; i++) p1[i] = (long)pol_to_small((GEN)vker[i]);
! 989: vker = p1;
! 990: }
! 991: pol = cgetg(N+3,t_POL);
! 992:
! 993: for (kk=1; kk<d; )
! 994: {
! 995: GEN polt;
! 996: if (p)
! 997: {
! 998: if (p==2)
! 999: {
! 1000: pol[2] = ((mymyrand() & 0x1000) == 0);
! 1001: pol[1] = evallgef(pol[2]? 3: 2);
! 1002: for (i=2; i<=d; i++)
! 1003: split_berlekamp_addmul(pol,(GEN)vker[i],(mymyrand()&0x1000)?0:1, p);
! 1004: }
! 1005: else
! 1006: {
! 1007: pol[2] = mymyrand()%p; /* vker[1] = 1 */
! 1008: pol[1] = evallgef(pol[2]? 3: 2);
! 1009: for (i=2; i<=d; i++)
! 1010: split_berlekamp_addmul(pol,(GEN)vker[i],mymyrand()%p, p);
! 1011: }
! 1012: polt = small_to_pol(pol,vu);
! 1013: }
! 1014: else
! 1015: {
! 1016: pol[2] = (long)genrand(pp);
! 1017: pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu);
! 1018: for (i=2; i<=d; i++)
! 1019: pol = gadd(pol, gmul((GEN)vker[i], genrand(pp)));
! 1020: polt = FpX_red(pol,pp);
! 1021: }
! 1022: for (i=1; i<=kk && kk<d; i++)
! 1023: {
! 1024: p1=t[i-1]; l1=degpol(p1);
! 1025: if (l1>1)
! 1026: {
! 1027: av = avma; p2 = FpX_res(polt, p1, pp);
! 1028: if (lgef(p2) <= 3) { avma=av; continue; }
! 1029: p2 = FpXQ_pow(p2,pps2, p1,pp);
! 1030: if (!signe(p2)) err(talker,"%Z not a prime in split_berlekamp",pp);
! 1031: p2 = ZX_s_add(p2, -1);
! 1032: p2 = FpX_gcd(p1,p2, pp); l2=degpol(p2);
! 1033: if (l2>0 && l2<l1)
! 1034: {
! 1035: p2 = FpX_normalize(p2, pp);
! 1036: t[i-1] = p2; kk++;
! 1037: t[kk-1] = FpX_div(p1,p2,pp);
! 1038: if (DEBUGLEVEL > 7) msgtimer("new factor");
! 1039: }
! 1040: else avma = av;
! 1041: }
! 1042: }
! 1043: }
! 1044: return d;
! 1045: }
! 1046:
! 1047: GEN
! 1048: factmod0(GEN f, GEN pp)
! 1049: {
! 1050: long i,j,k,e,p,N,nbfact,av = avma,tetpil,d;
! 1051: GEN pps2,ex,y,f2,p1,g1,u, *t;
! 1052:
! 1053: if (!(d = factmod_init(&f, pp, &p))) { avma=av; return trivfact(); }
! 1054: /* to hold factors and exponents */
! 1055: t = (GEN*)cgetg(d+1,t_VEC); ex = cgetg(d+1,t_VECSMALL);
! 1056: e = nbfact = 1;
! 1057: pps2 = shifti(pp,-1);
! 1058:
! 1059: for(;;)
! 1060: {
! 1061: f2 = FpX_gcd(f,derivpol(f), pp);
! 1062: g1 = lgef(f2)==3? f: FpX_div(f,f2,pp);
! 1063: k = 0;
! 1064: while (lgef(g1)>3)
! 1065: {
! 1066: k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }
! 1067: p1 = FpX_gcd(f2,g1, pp); u = g1; g1 = p1;
! 1068: if (lgef(p1)!=3)
! 1069: {
! 1070: u = FpX_div( u,p1,pp);
! 1071: f2= FpX_div(f2,p1,pp);
! 1072: }
! 1073: N = degpol(u);
! 1074: if (N)
! 1075: {
! 1076: /* here u is square-free (product of irred. of multiplicity e * k) */
! 1077: t[nbfact] = FpX_normalize(u,pp);
! 1078: d = (N==1)? 1: split_berlekamp(t+nbfact, pp, pps2);
! 1079: for (j=0; j<d; j++) ex[nbfact+j] = e*k;
! 1080: nbfact += d;
! 1081: }
! 1082: }
! 1083: if (!p) break;
! 1084: j=(degpol(f2))/p+3; if (j==3) break;
! 1085:
! 1086: e*=p; setlg(f,j); setlgef(f,j);
! 1087: for (i=2; i<j; i++) f[i] = f2[p*(i-2)+2];
! 1088: }
! 1089: tetpil=avma; y=cgetg(3,t_VEC);
! 1090: setlg((GEN)t, nbfact);
! 1091: setlg(ex, nbfact);
! 1092: y[1]=lcopy((GEN)t);
! 1093: y[2]=lcopy(ex);
! 1094: (void)sort_factor(y,cmpii);
! 1095: return gerepile(av,tetpil,y);
! 1096: }
! 1097: GEN
! 1098: factmod(GEN f, GEN pp)
! 1099: {
! 1100: long tetpil,av=avma;
! 1101: long nbfact;
! 1102: long j;
! 1103: GEN y,u,v;
! 1104: GEN z=factmod0(f,pp),t=(GEN)z[1],ex=(GEN)z[2];
! 1105: nbfact=lg(t);
! 1106: tetpil=avma; y=cgetg(3,t_MAT);
! 1107: u=cgetg(nbfact,t_COL); y[1]=(long)u;
! 1108: v=cgetg(nbfact,t_COL); y[2]=(long)v;
! 1109: for (j=1; j<nbfact; j++)
! 1110: {
! 1111: u[j] = (long)FpX((GEN)t[j], pp);
! 1112: v[j] = lstoi(ex[j]);
! 1113: }
! 1114: return gerepile(av,tetpil,y);
! 1115: }
! 1116: GEN
! 1117: factormod0(GEN f, GEN p, long flag)
! 1118: {
! 1119: switch(flag)
! 1120: {
! 1121: case 0: return factmod(f,p);
! 1122: case 1: return simplefactmod(f,p);
! 1123: default: err(flagerr,"factormod");
! 1124: }
! 1125: return NULL; /* not reached */
! 1126: }
! 1127:
! 1128: /*******************************************************************/
! 1129: /* */
! 1130: /* Recherche de racines p-adiques */
! 1131: /* */
! 1132: /*******************************************************************/
! 1133: /* make f suitable for [root|factor]padic */
! 1134: static GEN
! 1135: padic_pol_to_int(GEN f)
! 1136: {
! 1137: long i, l = lgef(f);
! 1138: f = gdiv(f,content(f));
! 1139: for (i=2; i<l; i++)
! 1140: switch(typ(f[i]))
! 1141: {
! 1142: case t_INT: break;
! 1143: case t_PADIC: f[i] = ltrunc((GEN)f[i]); break;
! 1144: default: err(talker,"incorrect coeffs in padic_pol_to_int");
! 1145: }
! 1146: return f;
! 1147: }
! 1148:
! 1149: /* return invlead * (x + O(pr)), x in Z or Z_p, pr = p^r */
! 1150: static GEN
! 1151: int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead)
! 1152: {
! 1153: GEN p1,y;
! 1154: long v,sx, av = avma;
! 1155:
! 1156: if (typ(x) == t_PADIC)
! 1157: {
! 1158: v = valp(x);
! 1159: if (r >= precp(x) + v) return invlead? gmul(x, invlead): gcopy(x);
! 1160: sx = !gcmp0(x);
! 1161: p1 = (GEN)x[4];
! 1162: }
! 1163: else
! 1164: {
! 1165: sx = signe(x);
! 1166: if (!sx) return gzero;
! 1167: v = pvaluation(x,p,&p1);
! 1168: }
! 1169: y = cgetg(5,t_PADIC);
! 1170: if (sx && v < r)
! 1171: {
! 1172: y[4] = lmodii(p1,pr); r -= v;
! 1173: }
! 1174: else
! 1175: {
! 1176: y[4] = zero; v = r; r = 0;
! 1177: }
! 1178: y[3] = (long)pr;
! 1179: y[2] = (long)p;
! 1180: y[1] = evalprecp(r)|evalvalp(v);
! 1181: return invlead? gerepileupto(av, gmul(invlead,y)): y;
! 1182: }
! 1183:
! 1184: /* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff
! 1185: * is a power of p), x in Z[X] (or Z_p[X]) */
! 1186: static GEN
! 1187: pol_to_padic(GEN x, GEN pr, GEN p, long r)
! 1188: {
! 1189: long v = 0,i,lx = lgef(x);
! 1190: GEN z = cgetg(lx,t_POL), lead = leading_term(x);
! 1191:
! 1192: if (gcmp1(lead)) lead = NULL;
! 1193: else
! 1194: {
! 1195: long av = avma;
! 1196: v = ggval(lead,p);
! 1197: if (v) lead = gdiv(lead, gpowgs(p,v));
! 1198: lead = int_to_padic(lead,p,pr,r,NULL);
! 1199: lead = gerepileupto(av, ginv(lead));
! 1200: }
! 1201: for (i=lx-1; i>1; i--)
! 1202: z[i] = (long)int_to_padic((GEN)x[i],p,pr,r,lead);
! 1203: z[1] = x[1]; return z;
! 1204: }
! 1205:
! 1206: static GEN
! 1207: padic_trivfact(GEN x, GEN p, long r)
! 1208: {
! 1209: GEN p1, y = cgetg(3,t_MAT);
! 1210: p1=cgetg(2,t_COL); y[1]=(long)p1; p = icopy(p);
! 1211: p1[1]=(long)pol_to_padic(x,gpowgs(p,r),p,r);
! 1212: p1=cgetg(2,t_COL); y[2]=(long)p1;
! 1213: p1[1]=un; return y;
! 1214: }
! 1215:
! 1216: /* a etant un p-adique, retourne le vecteur des racines p-adiques de f
! 1217: * congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
! 1218: * (ou a 4 si p=2).
! 1219: */
! 1220: GEN
! 1221: apprgen(GEN f, GEN a)
! 1222: {
! 1223: GEN fp,p1,p,P,pro,x,x2,u,ip;
! 1224: long av=avma,tetpil,v,Ps,i,j,k,lu,n,fl2;
! 1225:
! 1226: if (typ(f)!=t_POL) err(notpoler,"apprgen");
! 1227: if (gcmp0(f)) err(zeropoler,"apprgen");
! 1228: if (typ(a) != t_PADIC) err(rootper1);
! 1229: f = padic_pol_to_int(f);
! 1230: fp=derivpol(f); p1=ggcd(f,fp);
! 1231: if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
! 1232: p=(GEN)a[2]; p1=poleval(f,a);
! 1233: v=ggval(p1,p); if (v <= 0) err(rootper2);
! 1234: fl2=egalii(p,gdeux);
! 1235: if (fl2 && v==1) err(rootper2);
! 1236: v=ggval(poleval(fp,a),p);
! 1237: if (!v) /* simple zero */
! 1238: {
! 1239: while (!gcmp0(p1))
! 1240: {
! 1241: a = gsub(a,gdiv(p1,poleval(fp,a)));
! 1242: p1 = poleval(f,a);
! 1243: }
! 1244: tetpil=avma; pro=cgetg(2,t_VEC); pro[1]=lcopy(a);
! 1245: return gerepile(av,tetpil,pro);
! 1246: }
! 1247: n=degpol(f); pro=cgetg(n+1,t_VEC);
! 1248:
! 1249: if (is_bigint(p)) err(impl,"apprgen for p>=2^31");
! 1250: x = ggrandocp(p, valp(a) | precp(a));
! 1251: if (fl2)
! 1252: {
! 1253: x2=ggrandocp(p,2); P = stoi(4);
! 1254: }
! 1255: else
! 1256: {
! 1257: x2=ggrandocp(p,1); P = p;
! 1258: }
! 1259: f = poleval(f, gadd(a,gmul(P,polx[varn(f)])));
! 1260: if (!gcmp0(f)) f = gdiv(f,gpuigs(p,ggval(f, p)));
! 1261: Ps = itos(P);
! 1262: for (j=0,i=0; i<Ps; i++)
! 1263: {
! 1264: ip=stoi(i);
! 1265: if (gcmp0(poleval(f,gadd(ip,x2))))
! 1266: {
! 1267: u=apprgen(f,gadd(x,ip)); lu=lg(u);
! 1268: for (k=1; k<lu; k++)
! 1269: {
! 1270: j++; pro[j]=ladd(a,gmul(P,(GEN)u[k]));
! 1271: }
! 1272: }
! 1273: }
! 1274: setlg(pro,j+1); return gerepilecopy(av,pro);
! 1275: }
! 1276:
! 1277: /* Retourne le vecteur des racines p-adiques de f en precision r */
! 1278: GEN
! 1279: rootpadic(GEN f, GEN p, long r)
! 1280: {
! 1281: GEN fp,y,z,p1,pr,rac;
! 1282: long lx,i,j,k,n,av=avma,tetpil,fl2;
! 1283:
! 1284: if (typ(f)!=t_POL) err(notpoler,"rootpadic");
! 1285: if (gcmp0(f)) err(zeropoler,"rootpadic");
! 1286: if (r<=0) err(rootper4);
! 1287: f = padic_pol_to_int(f);
! 1288: fp=derivpol(f); p1=ggcd(f,fp);
! 1289: if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
! 1290: fl2=egalii(p,gdeux); rac=(fl2 && r>=2)? rootmod(f,stoi(4)): rootmod(f,p);
! 1291: lx=lg(rac); p=gclone(p);
! 1292: if (r==1)
! 1293: {
! 1294: tetpil=avma; y=cgetg(lx,t_COL);
! 1295: for (i=1; i<lx; i++)
! 1296: {
! 1297: z=cgetg(5,t_PADIC); y[i]=(long)z;
! 1298: z[1] = evalprecp(1)|evalvalp(0);
! 1299: z[2] = z[3] = (long)p;
! 1300: z[4] = lcopy(gmael(rac,i,2));
! 1301: }
! 1302: return gerepile(av,tetpil,y);
! 1303: }
! 1304: n=degpol(f); y=cgetg(n+1,t_COL);
! 1305: j=0; pr = NULL;
! 1306: z = cgetg(5,t_PADIC);
! 1307: z[2] = (long)p;
! 1308: for (i=1; i<lx; i++)
! 1309: {
! 1310: p1 = gmael(rac,i,2);
! 1311: if (signe(p1))
! 1312: {
! 1313: if (!fl2 || mod2(p1))
! 1314: {
! 1315: z[1] = evalvalp(0)|evalprecp(r);
! 1316: z[4] = (long)p1;
! 1317: }
! 1318: else
! 1319: {
! 1320: z[1] = evalvalp(1)|evalprecp(r);
! 1321: z[4] = un;
! 1322: }
! 1323: if (!pr) pr=gpuigs(p,r);
! 1324: z[3] = (long)pr;
! 1325: }
! 1326: else
! 1327: {
! 1328: z[1] = evalvalp(r);
! 1329: z[3] = un;
! 1330: z[4] = (long)p1;
! 1331: }
! 1332: p1 = apprgen(f,z);
! 1333: for (k=1; k<lg(p1); k++) y[++j]=p1[k];
! 1334: }
! 1335: setlg(y,j+1); return gerepilecopy(av,y);
! 1336: }
! 1337: /*************************************************************************/
! 1338: /* rootpadicfast */
! 1339: /*************************************************************************/
! 1340:
! 1341: /*lift accelerator. The author of the idea is unknown.*/
! 1342: long hensel_lift_accel(long n, long *pmask)
! 1343: {
! 1344: long a,j;
! 1345: long mask;
! 1346: mask=0;
! 1347: for(j=BITS_IN_LONG-1, a=n ;; j--)
! 1348: {
! 1349: mask|=(a&1)<<j;
! 1350: a=(a>>1)+(a&1);
! 1351: if (a==1) break;
! 1352: }
! 1353: *pmask=mask>>j;
! 1354: return BITS_IN_LONG-j;
! 1355: }
! 1356: /*
! 1357: SPEC:
! 1358: q is an integer > 1
! 1359: e>=0
! 1360: f in ZZ[X], with leading term prime to q.
! 1361: S must be a simple root mod p for all p|q.
! 1362:
! 1363: return roots of f mod q^e, as integers (implicitly mod Q)
! 1364: */
! 1365:
! 1366: /* STANDARD USE
! 1367: There exists p a prime number and a>0 such that
! 1368: q=p^a
! 1369: f in ZZ[X], with leading term prime to p.
! 1370: S must be a simple root mod p.
! 1371:
! 1372: return p-adics roots of f with prec b, as integers (implicitly mod q^e)
! 1373: */
! 1374:
! 1375: GEN
! 1376: rootpadiclift(GEN T, GEN S, GEN p, long e)
! 1377: {
! 1378: ulong ltop=avma;
! 1379: long x;
! 1380: GEN qold, q, qm1;
! 1381: GEN W, Tr, Sr, Wr = gzero;
! 1382: long i, nb, mask;
! 1383: x = varn(T);
! 1384: qold = p ; q = p; qm1 = gun;
! 1385: nb=hensel_lift_accel(e, &mask);
! 1386: Tr = FpX_red(T,q);
! 1387: W=FpX_eval(deriv(Tr, x),S,q);
! 1388: W=mpinvmod(W,q);
! 1389: for(i=0;i<nb;i++)
! 1390: {
! 1391: qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
! 1392: q = mulii(qm1, p);
! 1393: Tr = FpX_red(T,q);
! 1394: Sr = S;
! 1395: if (i)
! 1396: {
! 1397: W = modii(mulii(Wr,FpX_eval(deriv(Tr,x),Sr,qold)),qold);
! 1398: W = subii(gdeux,W);
! 1399: W = modii(mulii(Wr, W),qold);
! 1400: }
! 1401: Wr = W;
! 1402: S = subii(Sr, mulii(Wr, FpX_eval(Tr, Sr,q)));
! 1403: S = modii(S,q);
! 1404: qold = q;
! 1405: }
! 1406: return gerepileupto(ltop,S);
! 1407: }
! 1408: /*
! 1409: * Apply rootpadiclift to all roots in S and trace trick.
! 1410: * Elements of S must be distinct simple roots mod p for all p|q.
! 1411: */
! 1412:
! 1413: GEN
! 1414: rootpadicliftroots(GEN f, GEN S, GEN q, long e)
! 1415: {
! 1416: GEN y;
! 1417: long i,n=lg(S);
! 1418: if (n==1)
! 1419: return gcopy(S);
! 1420: y=cgetg(n,typ(S));
! 1421: for (i=1; i<n-1; i++)
! 1422: y[i]=(long) rootpadiclift(f, (GEN) S[i], q, e);
! 1423: if (n!=lgef(f)-2)/* non totally split*/
! 1424: y[n-1]=(long) rootpadiclift(f, (GEN) S[n-1], q, e);
! 1425: else/* distinct-->totally split-->use trace trick */
! 1426: {
! 1427: ulong av=avma;
! 1428: GEN z;
! 1429: z=(GEN)f[lgef(f)-2];/*-trace(roots)*/
! 1430: for(i=1; i<n-1;i++)
! 1431: z=addii(z,(GEN) y[i]);
! 1432: z=modii(negi(z),gpowgs(q,e));
! 1433: y[n-1]=lpileupto(av,z);
! 1434: }
! 1435: return y;
! 1436: }
! 1437: /*
! 1438: p is a prime number, pr a power of p,
! 1439:
! 1440: f in ZZ[X], with leading term prime to p.
! 1441: f must have no multiple roots mod p.
! 1442:
! 1443: return p-adics roots of f with prec pr, as integers (implicitly mod pr)
! 1444:
! 1445: */
! 1446: GEN
! 1447: rootpadicfast(GEN f, GEN p, long e)
! 1448: {
! 1449: ulong ltop=avma;
! 1450: GEN S,y;
! 1451: S=lift(rootmod(f,p));/*no multiplicity*/
! 1452: if (lg(S)==1)/*no roots*/
! 1453: {
! 1454: avma=ltop;
! 1455: return cgetg(1,t_COL);
! 1456: }
! 1457: S=gclone(S);
! 1458: avma=ltop;
! 1459: y=rootpadicliftroots(f,S,p,e);
! 1460: gunclone(S);
! 1461: return y;
! 1462: }
! 1463: /* Same as rootpadiclift for the polynomial X^n-a,
! 1464: * but here, n can be much larger.
! 1465: * TODO: generalize to sparse polynomials.
! 1466: */
! 1467: GEN
! 1468: padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)
! 1469: {
! 1470: ulong ltop=avma;
! 1471: GEN qold, q, qm1;
! 1472: GEN W, Sr, Wr = gzero;
! 1473: long i, nb, mask;
! 1474: qold = p ; q = p; qm1 = gun;
! 1475: nb = hensel_lift_accel(e, &mask);
! 1476: W = modii(mulii(n,powmodulo(S,subii(n,gun),q)),q);
! 1477: W = mpinvmod(W,q);
! 1478: for(i=0;i<nb;i++)
! 1479: {
! 1480: qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
! 1481: q = mulii(qm1, p);
! 1482: Sr = S;
! 1483: if (i)
! 1484: {
! 1485: W = modii(mulii(Wr,mulii(n,powmodulo(Sr,subii(n,gun),qold))),qold);
! 1486: W = subii(gdeux,W);
! 1487: W = modii(mulii(Wr, W),qold);
! 1488: }
! 1489: Wr = W;
! 1490: S = subii(Sr, mulii(Wr, subii(powmodulo(Sr,n,q),a)));
! 1491: S = modii(S,q);
! 1492: qold = q;
! 1493: }
! 1494: return gerepileupto(ltop,S);
! 1495: }
! 1496: /**************************************************************************/
! 1497: static long
! 1498: getprec(GEN x, long prec, GEN *p)
! 1499: {
! 1500: long i,e;
! 1501: GEN p1;
! 1502:
! 1503: for (i = lgef(x)-1; i>1; i--)
! 1504: {
! 1505: p1=(GEN)x[i];
! 1506: if (typ(p1)==t_PADIC)
! 1507: {
! 1508: e=valp(p1); if (signe(p1[4])) e += precp(p1);
! 1509: if (e<prec) prec = e; *p = (GEN)p1[2];
! 1510: }
! 1511: }
! 1512: return prec;
! 1513: }
! 1514:
! 1515: /* a appartenant a une extension finie de Q_p, retourne le vecteur des
! 1516: * racines de f congrues a a modulo p dans le cas ou on suppose f(a) congru a
! 1517: * 0 modulo p (ou a 4 si p=2).
! 1518: */
! 1519: GEN
! 1520: apprgen9(GEN f, GEN a)
! 1521: {
! 1522: GEN fp,p1,p,pro,x,x2,u,ip,t,vecg;
! 1523: long av=avma,tetpil,v,ps_1,i,j,k,lu,n,prec,d,va,fl2;
! 1524:
! 1525: if (typ(f)!=t_POL) err(notpoler,"apprgen9");
! 1526: if (gcmp0(f)) err(zeropoler,"apprgen9");
! 1527: if (typ(a)==t_PADIC) return apprgen(f,a);
! 1528: if (typ(a)!=t_POLMOD || typ(a[2])!=t_POL) err(rootper1);
! 1529: fp=derivpol(f); p1=ggcd(f,fp);
! 1530: if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }
! 1531: t=(GEN)a[1];
! 1532: prec = getprec((GEN)a[2], BIGINT, &p);
! 1533: prec = getprec(t, prec, &p);
! 1534: if (prec==BIGINT) err(rootper1);
! 1535:
! 1536: p1=poleval(f,a); v=ggval(lift_intern(p1),p); if (v<=0) err(rootper2);
! 1537: fl2=egalii(p,gdeux);
! 1538: if (fl2 && v==1) err(rootper2);
! 1539: v=ggval(lift_intern(poleval(fp,a)), p);
! 1540: if (!v)
! 1541: {
! 1542: while (!gcmp0(p1))
! 1543: {
! 1544: a = gsub(a, gdiv(p1,poleval(fp,a)));
! 1545: p1 = poleval(f,a);
! 1546: }
! 1547: tetpil=avma; pro=cgetg(2,t_COL); pro[1]=lcopy(a);
! 1548: return gerepile(av,tetpil,pro);
! 1549: }
! 1550: n=degpol(f); pro=cgetg(n+1,t_COL); j=0;
! 1551:
! 1552: if (is_bigint(p)) err(impl,"apprgen9 for p>=2^31");
! 1553: x=gmodulcp(ggrandocp(p,prec), t);
! 1554: if (fl2)
! 1555: {
! 1556: ps_1=3; x2=ggrandocp(p,2); p=stoi(4);
! 1557: }
! 1558: else
! 1559: {
! 1560: ps_1=itos(p)-1; x2=ggrandocp(p,1);
! 1561: }
! 1562: f = poleval(f,gadd(a,gmul(p,polx[varn(f)])));
! 1563: if (!gcmp0(f)) f=gdiv(f,gpuigs(p,ggval(f,p)));
! 1564: d=degpol(t); vecg=cgetg(d+1,t_COL);
! 1565: for (i=1; i<=d; i++)
! 1566: vecg[i] = (long)setloop(gzero);
! 1567: va=varn(t);
! 1568: for(;;) /* loop through F_q */
! 1569: {
! 1570: ip=gmodulcp(gtopoly(vecg,va),t);
! 1571: if (gcmp0(poleval(f,gadd(ip,x2))))
! 1572: {
! 1573: u=apprgen9(f,gadd(ip,x)); lu=lg(u);
! 1574: for (k=1; k<lu; k++)
! 1575: {
! 1576: j++; pro[j]=ladd(a,gmul(p,(GEN)u[k]));
! 1577: }
! 1578: }
! 1579: for (i=d; i; i--)
! 1580: {
! 1581: p1 = (GEN)vecg[i];
! 1582: if (p1[2] != ps_1) { (void)incloop(p1); break; }
! 1583: affsi(0, p1);
! 1584: }
! 1585: if (!i) break;
! 1586: }
! 1587: setlg(pro,j+1); return gerepilecopy(av,pro);
! 1588: }
! 1589:
! 1590: /*****************************************/
! 1591: /* Factorisation p-adique d'un polynome */
! 1592: /*****************************************/
! 1593: int
! 1594: cmp_padic(GEN x, GEN y)
! 1595: {
! 1596: long vx, vy;
! 1597: if (x == gzero) return -1;
! 1598: if (y == gzero) return 1;
! 1599: vx = valp(x);
! 1600: vy = valp(y);
! 1601: if (vx < vy) return 1;
! 1602: if (vx > vy) return -1;
! 1603: return cmpii((GEN)x[4], (GEN)y[4]);
! 1604: }
! 1605:
! 1606: /* factorise le polynome T=nf[1] dans Zp avec la precision pr */
! 1607: static GEN
! 1608: padicff2(GEN nf,GEN p,long pr)
! 1609: {
! 1610: long N=degpol(nf[1]),i,j,d,l;
! 1611: GEN mat,V,D,fa,p1,pk,dec_p,pke,a,theta;
! 1612:
! 1613: pk=gpuigs(p,pr); dec_p=primedec(nf,p);
! 1614: l=lg(dec_p); fa=cgetg(l,t_COL);
! 1615: for (i=1; i<l; i++)
! 1616: {
! 1617: p1 = (GEN)dec_p[i];
! 1618: pke = idealpows(nf,p1, pr * itos((GEN)p1[3]));
! 1619: p1=smith2(pke); V=(GEN)p1[3]; D=(GEN)p1[1];
! 1620: for (d=1; d<=N; d++)
! 1621: if (! egalii(gcoeff(V,d,d),pk)) break;
! 1622: a=ginv(D); theta=gmael(nf,8,2); mat=cgetg(d,t_MAT);
! 1623: for (j=1; j<d; j++)
! 1624: {
! 1625: p1 = gmul(D, element_mul(nf,theta,(GEN)a[j]));
! 1626: setlg(p1,d); mat[j]=(long)p1;
! 1627: }
! 1628: fa[i]=(long)caradj(mat,0,NULL);
! 1629: }
! 1630: a = cgetg(l,t_COL); pk = icopy(pk);
! 1631: for (i=1; i<l; i++)
! 1632: a[i] = (long)pol_to_padic((GEN)fa[i],pk,p,pr);
! 1633: return a;
! 1634: }
! 1635:
! 1636: static GEN
! 1637: padicff(GEN x,GEN p,long pr)
! 1638: {
! 1639: GEN q,basden,bas,invbas,mul,dx,nf,mat;
! 1640: long n=degpol(x),av=avma;
! 1641:
! 1642: nf=cgetg(10,t_VEC); nf[1]=(long)x; dx=discsr(x);
! 1643: mat=cgetg(3,t_MAT); mat[1]=lgetg(3,t_COL); mat[2]=lgetg(3,t_COL);
! 1644: coeff(mat,1,1)=(long)p; coeff(mat,1,2)=lstoi(pvaluation(dx,p,&q));
! 1645: coeff(mat,2,1)=(long)q; coeff(mat,2,2)=un;
! 1646: bas=allbase4(x,(long)mat,(GEN*)(nf+3),NULL);
! 1647: if (!carrecomplet(divii(dx,(GEN)nf[3]),(GEN*)(nf+4)))
! 1648: err(bugparier,"factorpadic2 (incorrect discriminant)");
! 1649: basden = get_bas_den(bas);
! 1650: invbas = QM_inv(vecpol_to_mat(bas,n), gun);
! 1651: mul = get_mul_table(x,basden,invbas,NULL);
! 1652: nf[7]=(long)bas;
! 1653: nf[8]=(long)invbas;
! 1654: nf[9]=(long)mul; nf[2]=nf[5]=nf[6]=zero;
! 1655: return gerepileupto(av,padicff2(nf,p,pr));
! 1656: }
! 1657:
! 1658: GEN
! 1659: factorpadic2(GEN x, GEN p, long r)
! 1660: {
! 1661: long av=avma,av2,k,i,j,i1,f,nbfac;
! 1662: GEN res,p1,p2,y,d,a,ap,t,v,w;
! 1663: GEN *fa;
! 1664:
! 1665: if (typ(x)!=t_POL) err(notpoler,"factorpadic");
! 1666: if (gcmp0(x)) err(zeropoler,"factorpadic");
! 1667: if (r<=0) err(rootper4);
! 1668:
! 1669: if (lgef(x)==3) return trivfact();
! 1670: if (!gcmp1(leading_term(x)))
! 1671: err(impl,"factorpadic2 for non-monic polynomial");
! 1672: if (lgef(x)==4) return padic_trivfact(x,p,r);
! 1673: y=cgetg(3,t_MAT);
! 1674: fa = (GEN*)new_chunk(lgef(x)-2);
! 1675: d=content(x); a=gdiv(x,d);
! 1676: ap=derivpol(a); t=ggcd(a,ap); v=gdeuc(a,t);
! 1677: w=gdeuc(ap,t); j=0; f=1; nbfac=0;
! 1678: while (f)
! 1679: {
! 1680: j++; w=gsub(w,derivpol(v)); f=signe(w);
! 1681: if (f) { res=ggcd(v,w); v=gdeuc(v,res); w=gdeuc(w,res); }
! 1682: else res=v;
! 1683: fa[j]=(lgef(res)>3) ? padicff(res,p,r) : cgetg(1,t_COL);
! 1684: nbfac += (lg(fa[j])-1);
! 1685: }
! 1686: av2=avma; y=cgetg(3,t_MAT);
! 1687: p1=cgetg(nbfac+1,t_COL); y[1]=(long)p1;
! 1688: p2=cgetg(nbfac+1,t_COL); y[2]=(long)p2;
! 1689: for (i=1,k=0; i<=j; i++)
! 1690: for (i1=1; i1<lg(fa[i]); i1++)
! 1691: {
! 1692: p1[++k]=lcopy((GEN)fa[i][i1]); p2[k]=lstoi(i);
! 1693: }
! 1694: y = gerepile(av,av2,y);
! 1695: sort_factor(y, cmp_padic); return y;
! 1696: }
! 1697:
! 1698: /*******************************************************************/
! 1699: /* */
! 1700: /* FACTORISATION P-adique avec ROUND 4 */
! 1701: /* */
! 1702: /*******************************************************************/
! 1703: extern GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
! 1704: extern GEN nilord(GEN p, GEN fx, long mf, GEN gx, long flag);
! 1705: extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pev, long e);
! 1706:
! 1707: static GEN
! 1708: squarefree(GEN f, GEN *ex)
! 1709: {
! 1710: GEN T,V,W,A,B;
! 1711: long n,i,k;
! 1712:
! 1713: T=ggcd(derivpol(f),f); V=gdeuc(f,T);
! 1714: n=lgef(f)-2; A=cgetg(n,t_COL); B=cgetg(n,t_COL);
! 1715: k=1; i=1;
! 1716: do
! 1717: {
! 1718: W=ggcd(T,V); T=gdeuc(T,W);
! 1719: if (lgef(V) != lgef(W))
! 1720: {
! 1721: A[i]=ldeuc(V,W); B[i]=k; i++;
! 1722: }
! 1723: k++; V=W;
! 1724: }
! 1725: while (lgef(V)>3);
! 1726: setlg(A,i); *ex=B; return A;
! 1727: }
! 1728:
! 1729: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 1730:
! 1731: /* reverse x in place */
! 1732: static void
! 1733: polreverse(GEN x)
! 1734: {
! 1735: long i, j;
! 1736: if (typ(x) != t_POL) err(typeer,"polreverse");
! 1737: for (i=2, j=lgef(x)-1; i<j; i++, j--) swap(x[i], x[j]);
! 1738: (void)normalizepol(x);
! 1739: }
! 1740:
! 1741: GEN
! 1742: factorpadic4(GEN f,GEN p,long prec)
! 1743: {
! 1744: GEN w,g,poly,fx,y,p1,p2,ex,pols,exps,ppow,lead;
! 1745: long v=varn(f),n=degpol(f),av,tetpil,mfx,i,k,j,r,pr;
! 1746: int reverse = 0;
! 1747:
! 1748: if (typ(f)!=t_POL) err(notpoler,"factorpadic");
! 1749: if (typ(p)!=t_INT) err(arither1);
! 1750: if (gcmp0(f)) err(zeropoler,"factorpadic");
! 1751: if (prec<=0) err(rootper4);
! 1752:
! 1753: if (n==0) return trivfact();
! 1754: av=avma; f = padic_pol_to_int(f);
! 1755: if (n==1) return gerepileupto(av, padic_trivfact(f,p,prec));
! 1756: lead = leading_term(f); pr = prec;
! 1757: if (!gcmp1(lead))
! 1758: {
! 1759: long val = ggval(lead,p), val1 = ggval(constant_term(f),p);
! 1760: if (val1 < val)
! 1761: {
! 1762: reverse = 1; polreverse(f);
! 1763: /* take care of loss of precision from leading coeff of factor
! 1764: * (whose valuation is <= val) */
! 1765: pr += val;
! 1766: val = val1;
! 1767: }
! 1768: pr += val * (n-1);
! 1769: }
! 1770: f = pol_to_monic(f, &lead);
! 1771:
! 1772: poly=squarefree(f,&ex);
! 1773: pols=cgetg(n+1,t_COL);
! 1774: exps=cgetg(n+1,t_COL); n = lg(poly);
! 1775: for (j=1,i=1; i<n; i++)
! 1776: {
! 1777: long av1 = avma;
! 1778: fx=(GEN)poly[i]; mfx=ggval(discsr(fx),p);
! 1779: w = (GEN)factmod(fx,p)[1];
! 1780: if (!mfx)
! 1781: { /* no repeated factors: Hensel lift */
! 1782: p1 = hensel_lift_fact(fx, lift_intern(w), NULL, p, gpowgs(p,pr), pr);
! 1783: p2 = stoi(ex[i]);
! 1784: for (k=1; k<lg(p1); k++,j++)
! 1785: {
! 1786: pols[j] = p1[k];
! 1787: exps[j] = (long)p2;
! 1788: }
! 1789: continue;
! 1790: }
! 1791: /* use Round 4 */
! 1792: r = lg(w)-1;
! 1793: g = lift_intern((GEN)w[r]);
! 1794: p2 = (r == 1)? nilord(p,fx,mfx,g,pr)
! 1795: : Decomp(p,fx,mfx,polx[v],fx,g, (pr<=mfx)? mfx+1: pr);
! 1796: if (p2)
! 1797: {
! 1798: p2 = gerepileupto(av1,p2);
! 1799: p1 = (GEN)p2[1];
! 1800: p2 = (GEN)p2[2];
! 1801: for (k=1; k<lg(p1); k++,j++)
! 1802: {
! 1803: pols[j]=p1[k];
! 1804: exps[j]=lmulis((GEN)p2[k],ex[i]);
! 1805: }
! 1806: }
! 1807: else
! 1808: {
! 1809: avma=av1;
! 1810: pols[j]=(long)fx;
! 1811: exps[j]=lstoi(ex[i]); j++;
! 1812: }
! 1813: }
! 1814: if (lead)
! 1815: {
! 1816: p1 = gmul(polx[v],lead);
! 1817: for (i=1; i<j; i++)
! 1818: {
! 1819: p2 = poleval((GEN)pols[i], p1);
! 1820: pols[i] = ldiv(p2, content(p2));
! 1821: }
! 1822: }
! 1823:
! 1824: tetpil=avma; y=cgetg(3,t_MAT);
! 1825: p1 = cgetg(j,t_COL); ppow = gpowgs(p,prec); p = icopy(p);
! 1826: for (i=1; i<j; i++)
! 1827: {
! 1828: if (reverse) polreverse((GEN)pols[i]);
! 1829: p1[i] = (long)pol_to_padic((GEN)pols[i],ppow,p,prec);
! 1830: }
! 1831: y[1]=(long)p1; setlg(exps,j);
! 1832: y[2]=lcopy(exps); y = gerepile(av,tetpil,y);
! 1833: sort_factor(y, cmp_padic); return y;
! 1834: }
! 1835:
! 1836: GEN
! 1837: factorpadic0(GEN f,GEN p,long r,long flag)
! 1838: {
! 1839: switch(flag)
! 1840: {
! 1841: case 0: return factorpadic4(f,p,r);
! 1842: case 1: return factorpadic2(f,p,r);
! 1843: default: err(flagerr,"factorpadic");
! 1844: }
! 1845: return NULL; /* not reached */
! 1846: }
! 1847:
! 1848: /*******************************************************************/
! 1849: /* */
! 1850: /* FACTORISATION DANS F_q */
! 1851: /* */
! 1852: /*******************************************************************/
! 1853: extern GEN to_Kronecker(GEN P, GEN Q);
! 1854: extern GEN from_Kronecker(GEN z, GEN pol);
! 1855: static GEN spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S);
! 1856:
! 1857: static GEN
! 1858: to_fq(GEN x, GEN a, GEN p)
! 1859: {
! 1860: long i,lx = lgef(x);
! 1861: GEN z = cgetg(3,t_POLMOD), pol = cgetg(lx,t_POL);
! 1862: pol[1] = x[1];
! 1863: if (lx == 2) setsigne(pol, 0);
! 1864: else
! 1865: for (i=2; i<lx; i++) pol[i] = (long)mod((GEN)x[i], p);
! 1866: /* assume deg(pol) < deg(a) */
! 1867: z[1] = (long)a;
! 1868: z[2] = (long)pol; return z;
! 1869: }
! 1870:
! 1871: /* x POLMOD over Fq, return lift(x^n) */
! 1872: static GEN
! 1873: Kronecker_powmod(GEN x, GEN mod, GEN n)
! 1874: {
! 1875: long lim,av,av0 = avma, i,j,m,v = varn(x);
! 1876: GEN y, p1, p = NULL, pol = NULL;
! 1877:
! 1878: for (i=lgef(mod)-1; i>1; i--)
! 1879: {
! 1880: p1 = (GEN)mod[i];
! 1881: if (typ(p1) == t_POLMOD) { pol = (GEN)p1[1] ; break; }
! 1882: }
! 1883: if (!pol) err(talker,"need POLMOD coeffs in Kronecker_powmod");
! 1884: for (i=lgef(pol)-1; i>1; i--)
! 1885: {
! 1886: p1 = (GEN)pol[i];
! 1887: if (typ(p1) == t_INTMOD) { p = (GEN)p1[1] ; break; }
! 1888: }
! 1889: if (!p) err(talker,"need Fq coeffs in Kronecker_powmod");
! 1890: x = lift_intern(to_Kronecker(x,pol));
! 1891:
! 1892: /* adapted from powgi */
! 1893: av=avma; lim=stack_lim(av,1);
! 1894: p1 = n+2; m = *p1;
! 1895:
! 1896: y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
! 1897: for (i=lgefint(n)-2;;)
! 1898: {
! 1899: for (; j; m<<=1,j--)
! 1900: {
! 1901: y = gsqr(y);
! 1902: y = from_Kronecker(FpX(y,p), pol);
! 1903: setvarn(y, v);
! 1904: y = gres(y, mod);
! 1905: y = lift_intern(to_Kronecker(y,pol));
! 1906:
! 1907: if (m<0)
! 1908: {
! 1909: y = gmul(y,x);
! 1910: y = from_Kronecker(FpX(y,p), pol);
! 1911: setvarn(y, v);
! 1912: y = gres(y, mod);
! 1913: y = lift_intern(to_Kronecker(y,pol));
! 1914: }
! 1915: if (low_stack(lim, stack_lim(av,1)))
! 1916: {
! 1917: if(DEBUGMEM>1) err(warnmem,"Kronecker_powmod");
! 1918: y = gerepilecopy(av, y);
! 1919: }
! 1920: }
! 1921: if (--i == 0) break;
! 1922: m = *++p1, j = BITS_IN_LONG;
! 1923: }
! 1924: y = from_Kronecker(FpX(y,p),pol);
! 1925: setvarn(y, v); return gerepileupto(av0, y);
! 1926: }
! 1927:
! 1928: GEN
! 1929: FpX_rand(long d1, long v, GEN p)
! 1930: {
! 1931: long i, d = d1+2;
! 1932: GEN y;
! 1933: y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
! 1934: for (i=2; i<d; i++) y[i] = (long)genrand(p);
! 1935: normalizepol_i(y,d); return y;
! 1936: }
! 1937:
! 1938: /* return a random polynomial in F_q[v], degree < d1 */
! 1939: static GEN
! 1940: FqX_rand(long d1, long v, GEN p, GEN a)
! 1941: {
! 1942: long i,j, d = d1+2, k = lgef(a)-1;
! 1943: GEN y,t;
! 1944:
! 1945: y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
! 1946: t = cgetg(k,t_POL); t[1] = a[1];
! 1947: for (i=2; i<d; i++)
! 1948: {
! 1949: for (j=2; j<k; j++) t[j] = (long)genrand(p);
! 1950: normalizepol_i(t,k); y[i]=(long)to_fq(t,a,p);
! 1951: }
! 1952: normalizepol_i(y,d); return y;
! 1953: }
! 1954:
! 1955: /* split into r factors of degree d */
! 1956: static void
! 1957: split9(GEN *t, long d, GEN p, GEN q, GEN a, GEN S)
! 1958: {
! 1959: long l,v,av,is2,cnt, dt = degpol(*t), da = degpol(a);
! 1960: GEN w,w0;
! 1961:
! 1962: if (dt == d) return;
! 1963: v = varn(*t);
! 1964: if (DEBUGLEVEL > 6) timer2();
! 1965: av = avma; is2 = egalii(p, gdeux);
! 1966: for(cnt = 1;;cnt++)
! 1967: { /* splits *t with probability ~ 1 - 2^(1-r) */
! 1968: w = w0 = FqX_rand(dt,v, p,a);
! 1969: for (l=1; l<d; l++) /* sum_{0<i<d} w^(q^i), result in (F_q)^r */
! 1970: w = gadd(w0, spec_Fq_pow_mod_pol(w, p, a, S));
! 1971: if (is2)
! 1972: {
! 1973: w0 = w;
! 1974: for (l=1; l<da; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */
! 1975: w = gadd(w0, gres(gsqr(w), *t));
! 1976: }
! 1977: else
! 1978: {
! 1979: w = Kronecker_powmod(w, *t, shifti(q,-1));
! 1980: /* w in {-1,0,1}^r */
! 1981: if (lgef(w) == 3) continue;
! 1982: w[2] = ladd((GEN)w[2], gun);
! 1983: }
! 1984: w = ggcd(*t,w); l = degpol(w);
! 1985: if (l && l != dt) break;
! 1986: avma = av;
! 1987: }
! 1988: w = gerepileupto(av,w);
! 1989: if (DEBUGLEVEL > 6)
! 1990: fprintferr("[split9] time for splitting: %ld (%ld trials)\n",timer2(),cnt);
! 1991: l /= d; t[l]=gdeuc(*t,w); *t=w;
! 1992: split9(t+l,d,p,q,a,S);
! 1993: split9(t ,d,p,q,a,S);
! 1994: }
! 1995:
! 1996: /* to "compare" (real) scalars and t_INTMODs */
! 1997: static int
! 1998: cmp_coeff(GEN x, GEN y)
! 1999: {
! 2000: if (typ(x) == t_INTMOD) x = (GEN)x[2];
! 2001: if (typ(y) == t_INTMOD) y = (GEN)y[2];
! 2002: return gcmp(x,y);
! 2003: }
! 2004:
! 2005: int
! 2006: cmp_pol(GEN x, GEN y)
! 2007: {
! 2008: long fx[3], fy[3];
! 2009: long i,lx,ly;
! 2010: int fl;
! 2011: if (typ(x) == t_POLMOD) x = (GEN)x[2];
! 2012: if (typ(y) == t_POLMOD) y = (GEN)y[2];
! 2013: if (typ(x) == t_POL) lx = lgef(x); else { lx = 3; fx[2] = (long)x; x = fx; }
! 2014: if (typ(y) == t_POL) ly = lgef(y); else { ly = 3; fy[2] = (long)y; y = fy; }
! 2015: if (lx > ly) return 1;
! 2016: if (lx < ly) return -1;
! 2017: for (i=lx-1; i>1; i--)
! 2018: if ((fl = cmp_coeff((GEN)x[i], (GEN)y[i]))) return fl;
! 2019: return 0;
! 2020: }
! 2021:
! 2022: /* assume n > 1, X a POLMOD over Fq */
! 2023: /* return S = [ X^q, X^2q, ... X^(n-1)q ] mod T (in Fq[X]) in Kronecker form */
! 2024: static GEN
! 2025: init_pow_q_mod_pT(GEN Xmod, GEN q, GEN a, GEN T)
! 2026: {
! 2027: long i, n = degpol(T);
! 2028: GEN p1, S = cgetg(n, t_VEC);
! 2029:
! 2030: S[1] = (long)Kronecker_powmod((GEN)Xmod[2], (GEN)Xmod[1], q);
! 2031: #if 1 /* use as many squarings as possible */
! 2032: for (i=2; i < n; i+=2)
! 2033: {
! 2034: p1 = gsqr((GEN)S[i>>1]);
! 2035: S[i] = lres(p1, T);
! 2036: if (i == n-1) break;
! 2037: p1 = gmul((GEN)S[i], (GEN)S[1]);
! 2038: S[i+1] = lres(p1, T);
! 2039: }
! 2040: #else
! 2041: for (i=2; i < n; i++)
! 2042: {
! 2043: p1 = gmul((GEN)S[i-1], (GEN)S[1]);
! 2044: S[i] = lres(p1, T);
! 2045: }
! 2046: #endif
! 2047: for (i=1; i < n; i++)
! 2048: S[i] = (long)lift_intern(to_Kronecker((GEN)S[i], a));
! 2049: return S;
! 2050: }
! 2051:
! 2052: /* compute x^q, S is as above */
! 2053: static GEN
! 2054: spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)
! 2055: {
! 2056: long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);
! 2057: GEN x0 = x+2, z,c;
! 2058:
! 2059: c = (GEN)x0[0];
! 2060: z = lift_intern(lift(c));
! 2061: for (i = 1; i <= dx; i++)
! 2062: {
! 2063: GEN d;
! 2064: c = (GEN)x0[i];
! 2065: if (gcmp0(c)) continue;
! 2066: d = (GEN)S[i];
! 2067: if (!gcmp1(c)) d = gmul(lift_intern(lift(c)),d);
! 2068: z = gadd(z, d);
! 2069: if (low_stack(lim, stack_lim(av,1)))
! 2070: {
! 2071: if(DEBUGMEM>1) err(warnmem,"spec_Fq_pow_mod_pol");
! 2072: z = gerepileupto(av, z);
! 2073: }
! 2074: }
! 2075: z = FpX(z, p);
! 2076: z = from_Kronecker(z, a);
! 2077: setvarn(z, varn(x)); return gerepileupto(av, z);
! 2078: }
! 2079:
! 2080: static long isabsolutepol(GEN f, GEN pp, GEN a)
! 2081: {
! 2082: int i,res=1;
! 2083: GEN c;
! 2084: for(i=2; i<lg(f); i++)
! 2085: {
! 2086: c = (GEN) f[i];
! 2087: switch(typ(c))
! 2088: {
! 2089: case t_INT: /* OK*/
! 2090: break;
! 2091: case t_INTMOD:
! 2092: if (gcmp((GEN)c[1],pp))
! 2093: err(typeer,"factmod9");
! 2094: break;
! 2095: case t_POLMOD:
! 2096: if (gcmp((GEN)c[1],a))
! 2097: err(typeer,"factmod9");
! 2098: isabsolutepol((GEN)c[1],pp,gzero);
! 2099: isabsolutepol((GEN)c[2],pp,gzero);
! 2100: if (degpol(c[1])>0)
! 2101: res = 0;
! 2102: break;
! 2103: case t_POL:
! 2104: isabsolutepol(c,pp,gzero);
! 2105: if (degpol(c)>0)
! 2106: res = 0;
! 2107: break;
! 2108: default:
! 2109: err(typeer,"factmod9");
! 2110: }
! 2111: }
! 2112: return res;
! 2113: }
! 2114:
! 2115: GEN
! 2116: factmod9(GEN f, GEN pp, GEN a)
! 2117: {
! 2118: long av = avma, tetpil,p,i,j,k,d,e,vf,va,nbfact,nbf,pk;
! 2119: GEN S,ex,y,f2,f3,df1,df2,g,g1,xmod,u,v,qqd,qq,unfp,unfq, *t;
! 2120: GEN frobinv,X;
! 2121:
! 2122: if (typ(a)!=t_POL || typ(f)!=t_POL || gcmp0(a)) err(typeer,"factmod9");
! 2123: vf=varn(f); va=varn(a);
! 2124: if (va<=vf) err(talker,"polynomial variable must be of higher priority than finite field\nvariable in factorff");
! 2125: if (isabsolutepol(f, pp, a))
! 2126: {
! 2127: GEN z= Fp_factor_rel0(simplify(lift(lift(f))), pp, lift(a));
! 2128: GEN t=(GEN)z[1],ex=(GEN)z[2];
! 2129: unfp = gmodulsg(1,pp);
! 2130: unfq = gmodulcp(gmul(unfp, polun[va]), gmul(unfp,a));
! 2131: nbfact=lg(t);
! 2132: tetpil=avma;
! 2133: y=cgetg(3,t_MAT);
! 2134: u=cgetg(nbfact,t_COL); y[1]=(long)u;
! 2135: v=cgetg(nbfact,t_COL); y[2]=(long)v;
! 2136: for (j=1; j<nbfact; j++)
! 2137: {
! 2138: u[j] = lmul((GEN)t[j],unfq);
! 2139: v[j] = lstoi(ex[j]);
! 2140: }
! 2141: return gerepile(av,tetpil,y);
! 2142: }
! 2143: p = is_bigint(pp)? 0: itos(pp);
! 2144: unfp=gmodulsg(1,pp); a=gmul(unfp,a);
! 2145: unfq=gmodulo(gmul(unfp,polun[va]), a); a = (GEN)unfq[1];
! 2146: f = gmul(unfq,f); if (!signe(f)) err(zeropoler,"factmod9");
! 2147: d = degpol(f); if (!d) { avma=av; gunclone(a); return trivfact(); }
! 2148:
! 2149: S = df2 = NULL; /* gcc -Wall */
! 2150: pp = gmael(a,2,1); /* out of the stack */
! 2151: t = (GEN*)cgetg(d+1,t_VEC); ex = new_chunk(d+1);
! 2152:
! 2153: frobinv = gpowgs(pp, lgef(a)-4);
! 2154: xmod = cgetg(3,t_POLMOD);
! 2155: X = gmul(polx[vf],unfq);
! 2156: xmod[2] = (long)X;
! 2157: qq=gpuigs(pp,degpol(a));
! 2158: e = nbfact = 1;
! 2159: pk=1; df1=derivpol(f); f3=NULL;
! 2160: for(;;)
! 2161: {
! 2162: long du,dg;
! 2163: while (gcmp0(df1))
! 2164: { /* needs d >= pp: p = 0 can't happen */
! 2165: pk *= p; e=pk;
! 2166: j=(degpol(f))/p+3; setlg(f,j); setlgef(f,j);
! 2167: for (i=2; i<j; i++) f[i] = (long)powgi((GEN)f[p*(i-2)+2], frobinv);
! 2168: df1=derivpol(f); f3=NULL;
! 2169: }
! 2170: f2 = f3? f3: ggcd(f,df1);
! 2171: if (lgef(f2)==3) u = f;
! 2172: else
! 2173: {
! 2174: g1=gdeuc(f,f2); df2=derivpol(f2);
! 2175: if (gcmp0(df2)) { u=g1; f3=f2; }
! 2176: else
! 2177: {
! 2178: f3=ggcd(f2,df2);
! 2179: if (lgef(f3)==3) u=gdeuc(g1,f2);
! 2180: else
! 2181: u=gdeuc(g1,gdeuc(f2,f3));
! 2182: }
! 2183: }
! 2184: /* u is square-free (product of irreducibles of multiplicity e) */
! 2185: qqd=gun; xmod[1]=(long)u;
! 2186:
! 2187: du = degpol(u); v = X;
! 2188: if (du > 1) S = init_pow_q_mod_pT(xmod, qq, a, u);
! 2189: for (d=1; d <= du>>1; d++)
! 2190: {
! 2191: qqd=mulii(qqd,qq);
! 2192: v = spec_Fq_pow_mod_pol(v, pp, a, S);
! 2193: g = ggcd(gsub(v,X),u);
! 2194: dg = degpol(g);
! 2195: if (dg <= 0) continue;
! 2196:
! 2197: /* all factors of g have degree d */
! 2198: j = nbfact+dg/d;
! 2199:
! 2200: t[nbfact] = g;
! 2201: split9(t+nbfact,d,pp,qq,a,S);
! 2202: for (; nbfact<j; nbfact++) ex[nbfact]=e;
! 2203: du -= dg;
! 2204: u = gdeuc(u,g);
! 2205: v = gres(v,u); xmod[1] = (long)u;
! 2206: }
! 2207: if (du) { t[nbfact]=u; ex[nbfact++]=e; }
! 2208: if (lgef(f2) == 3) break;
! 2209:
! 2210: f=f2; df1=df2; e += pk;
! 2211: }
! 2212:
! 2213: nbf=nbfact; tetpil=avma; y=cgetg(3,t_MAT);
! 2214: for (j=1; j<nbfact; j++)
! 2215: {
! 2216: t[j]=gdiv((GEN)t[j],leading_term(t[j]));
! 2217: for (k=1; k<j; k++)
! 2218: if (ex[k] && gegal(t[j],t[k]))
! 2219: {
! 2220: ex[k] += ex[j]; ex[j]=0;
! 2221: nbf--; break;
! 2222: }
! 2223: }
! 2224: u=cgetg(nbf,t_COL); y[1]=(long)u;
! 2225: v=cgetg(nbf,t_COL); y[2]=(long)v;
! 2226: for (j=1,k=0; j<nbfact; j++)
! 2227: if (ex[j])
! 2228: {
! 2229: k++;
! 2230: u[k]=(long)t[j];
! 2231: v[k]=lstoi(ex[j]);
! 2232: }
! 2233: y = gerepile(av,tetpil,y);
! 2234: u=(GEN)y[1];
! 2235: { /* put a back on the stack */
! 2236: GEN tokill = a;
! 2237: a = forcecopy(a);
! 2238: gunclone(tokill);
! 2239: }
! 2240: pp = (GEN)leading_term(a)[1];
! 2241: for (j=1; j<nbf; j++) fqunclone((GEN)u[j], a, pp);
! 2242: (void)sort_factor(y, cmp_pol); return y;
! 2243: }
! 2244: /* See also: Isomorphisms between finite field and relative
! 2245: * factorization in polarit3.c */
! 2246:
! 2247: /*******************************************************************/
! 2248: /* */
! 2249: /* RACINES COMPLEXES */
! 2250: /* l represente la longueur voulue pour les parties */
! 2251: /* reelles et imaginaires des racines de x */
! 2252: /* */
! 2253: /*******************************************************************/
! 2254: GEN square_free_factorization(GEN pol);
! 2255: static GEN laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC);
! 2256: GEN zrhqr(GEN a,long PREC);
! 2257:
! 2258: GEN
! 2259: rootsold(GEN x, long l)
! 2260: {
! 2261: long av1=avma,i,j,f,g,gg,fr,deg,l0,l1,l2,l3,l4,ln;
! 2262: long exc,expmin,m,deg0,k,ti,h,ii,e,e1,emax,v;
! 2263: GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7,p8;
! 2264: GEN p9,p10,p11,p12,p14,p15,pa,pax,pb,pp,pq,ps;
! 2265:
! 2266: if (typ(x)!=t_POL) err(typeer,"rootsold");
! 2267: v=varn(x); deg0=degpol(x); expmin=12 - bit_accuracy(l);
! 2268: if (!signe(x)) err(zeropoler,"rootsold");
! 2269: y=cgetg(deg0+1,t_COL); if (!deg0) return y;
! 2270: for (i=1; i<=deg0; i++)
! 2271: {
! 2272: p1=cgetg(3,t_COMPLEX); p1[1]=lgetr(l); p1[2]=lgetr(l); y[i]=(long)p1;
! 2273: for (j=3; j<l; j++) ((GEN)p1[2])[j]=((GEN)p1[1])[j]=0;
! 2274: }
! 2275: g=1; gg=1;
! 2276: for (i=2; i<=deg0+2; i++)
! 2277: {
! 2278: ti=typ(x[i]);
! 2279: if (ti==t_REAL) gg=0;
! 2280: else if (ti==t_QUAD)
! 2281: {
! 2282: p2=gmael3(x,i,1,2);
! 2283: if (gsigne(p2)>0) g=0;
! 2284: } else if (ti != t_INT && ti != t_INTMOD && !is_frac_t(ti)) g=0;
! 2285: }
! 2286: l1=avma; p2=cgetg(3,t_COMPLEX);
! 2287: p2[1]=lmppi(DEFAULTPREC); p2[2]=ldivrs((GEN)p2[1],10);
! 2288: p11=cgetg(4,t_POL); p11[1]=evalsigne(1)+evallgef(4);
! 2289: setvarn(p11,v); p11[3]=un;
! 2290:
! 2291: p12=cgetg(5,t_POL); p12[1]=evalsigne(1)+evallgef(5);
! 2292: setvarn(p12,v); p12[4]=un;
! 2293: for (i=2; i<=deg0+2 && gcmp0((GEN)x[i]); i++) gaffsg(0,(GEN)y[i-1]);
! 2294: k=i-2;
! 2295: if (k!=deg0)
! 2296: {
! 2297: if (k)
! 2298: {
! 2299: j=deg0+3-k; pax=cgetg(j,t_POL);
! 2300: pax[1] = evalsigne(1) | evalvarn(v) | evallgef(j);
! 2301: for (i=2; i<j; i++) pax[i]=x[i+k];
! 2302: }
! 2303: else pax=x;
! 2304: xd0=deriv(pax,v); m=1; pa=pax;
! 2305: pq = NULL; /* for lint */
! 2306: if (gg) { pp=ggcd(pax,xd0); h=isnonscalar(pp); if (h) pq=gdeuc(pax,pp); }
! 2307: else{ pp=gun; h=0; }
! 2308: do
! 2309: {
! 2310: if (h)
! 2311: {
! 2312: pa=pp; pb=pq; pp=ggcd(pa,deriv(pa,v)); h=isnonscalar(pp);
! 2313: if (h) pq=gdeuc(pa,pp); else pq=pa; ps=gdeuc(pb,pq);
! 2314: }
! 2315: else ps=pa;
! 2316: /* calcul des racines d'ordre exactement m */
! 2317: deg=degpol(ps);
! 2318: if (deg)
! 2319: {
! 2320: l3=avma; e=gexpo((GEN)ps[deg+2]); emax=e;
! 2321: for (i=2; i<deg+2; i++)
! 2322: {
! 2323: p3=(GEN)(ps[i]);
! 2324: e1=gexpo(p3); if (e1>emax) emax=e1;
! 2325: }
! 2326: e=emax-e; if (e<0) e=0; avma=l3; if (ps!=pax) xd0=deriv(ps,v);
! 2327: xdabs=cgetg(deg+2,t_POL); xdabs[1]=xd0[1];
! 2328: for (i=2; i<deg+2; i++)
! 2329: {
! 2330: l3=avma; p3=(GEN)xd0[i];
! 2331: p4=gabs(greal(p3),l);
! 2332: p5=gabs(gimag(p3),l); l4=avma;
! 2333: xdabs[i]=lpile(l3,l4,gadd(p4,p5));
! 2334: }
! 2335: l0=avma; xc=gcopy(ps); xd=gcopy(xd0); l2=avma;
! 2336: for (i=1; i<=deg; i++)
! 2337: {
! 2338: if (i==deg)
! 2339: {
! 2340: p1=(GEN)y[k+m*i]; gdivz(gneg_i((GEN)xc[2]),(GEN)xc[3],p1);
! 2341: p14=(GEN)p1[1]; p15=(GEN)p1[2];
! 2342: }
! 2343: else
! 2344: {
! 2345: p3=gshift(p2,e); p4=poleval(xc,p3); p5=gnorm(p4); exc=0;
! 2346: while (exc >= -20)
! 2347: {
! 2348: p7 = gneg_i(gdiv(p4, poleval(xd,p3)));
! 2349: l3 = avma;
! 2350: if (gcmp0(p5)) exc = -32;
! 2351: else exc = expo(gnorm(p7))-expo(gnorm(p3));
! 2352: avma = l3;
! 2353: for (j=1; j<=10; j++)
! 2354: {
! 2355: p8=gadd(p3,p7); p9=poleval(xc,p8); p10=gnorm(p9);
! 2356: if (exc < -20 || cmprr(p10,p5) < 0)
! 2357: {
! 2358: GEN *gptr[3];
! 2359: p3=p8; p4=p9; p5=p10;
! 2360: gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5;
! 2361: gerepilemanysp(l2,l3,gptr,3);
! 2362: break;
! 2363: }
! 2364: gshiftz(p7,-2,p7); avma=l3;
! 2365: }
! 2366: if (j > 10)
! 2367: {
! 2368: avma=av1;
! 2369: if (DEBUGLEVEL)
! 2370: {
! 2371: fprintferr("too many iterations in rootsold(): ");
! 2372: fprintferr("using roots2()\n"); flusherr();
! 2373: }
! 2374: return roots2(x,l);
! 2375: }
! 2376: }
! 2377: p1=(GEN)y[k+m*i]; setlg(p1[1],3); setlg(p1[2],3); gaffect(p3,p1);
! 2378: avma=l2; p14=(GEN)p1[1]; p15=(GEN)p1[2];
! 2379: for (ln=4; ln<=l; ln=(ln<<1)-2)
! 2380: {
! 2381: setlg(p14,ln); setlg(p15,ln);
! 2382: if (gcmp0(p14)) { settyp(p14,t_INT); p14[1]=2; }
! 2383: if (gcmp0(p15)) { settyp(p15,t_INT); p15[1]=2; }
! 2384: p4=poleval(xc,p1);
! 2385: p5=poleval(xd,p1); p6=gneg_i(gdiv(p4,p5));
! 2386: settyp(p14,t_REAL); settyp(p15,t_REAL);
! 2387: gaffect(gadd(p1,p6),p1); avma=l2;
! 2388: }
! 2389: }
! 2390: setlg(p14,l); setlg(p15,l);
! 2391: p7=gcopy(p1); p14=(GEN)(p7[1]); p15=(GEN)(p7[2]);
! 2392: setlg(p14,l+1); setlg(p15,l+1);
! 2393: if (gcmp0(p14)) { settyp(p14,t_INT); p14[1]=2; }
! 2394: if (gcmp0(p15)) { settyp(p15,t_INT); p15[1]=2; }
! 2395: for (ii=1; ii<=5; ii++)
! 2396: {
! 2397: p4=poleval(ps,p7); p5=poleval(xd0,p7);
! 2398: p6=gneg_i(gdiv(p4,p5)); p7=gadd(p7,p6);
! 2399: p14=(GEN)(p7[1]); p15=(GEN)(p7[2]);
! 2400: if (gcmp0(p14)) { settyp(p14,t_INT); p14[1]=2; }
! 2401: if (gcmp0(p15)) { settyp(p15,t_INT); p15[1]=2; }
! 2402: }
! 2403: gaffect(p7,p1); p4=poleval(ps,p7);
! 2404: p6=gdiv(p4,poleval(xdabs,gabs(p7,l)));
! 2405: if (gexpo(p6)>=expmin)
! 2406: {
! 2407: avma=av1;
! 2408: if (DEBUGLEVEL)
! 2409: {
! 2410: fprintferr("internal error in rootsold(): using roots2()\n");
! 2411: flusherr();
! 2412: }
! 2413: return roots2(x,l);
! 2414: }
! 2415: avma=l2;
! 2416: if (expo(p1[2])<expmin && g)
! 2417: {
! 2418: gaffect(gzero,(GEN)p1[2]);
! 2419: for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]);
! 2420: p11[2]=lneg((GEN)p1[1]);
! 2421: l4=avma; xc=gerepile(l0,l4,gdeuc(xc,p11));
! 2422: }
! 2423: else
! 2424: {
! 2425: for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]);
! 2426: if (g)
! 2427: {
! 2428: p1=gconj(p1);
! 2429: for (j=1; j<=m; j++) gaffect(p1,(GEN)y[k+i*m+j]);
! 2430: i++;
! 2431: p12[2]=lnorm(p1); p12[3]=lmulsg(-2,(GEN)p1[1]); l4=avma;
! 2432: xc=gerepile(l0,l4,gdeuc(xc,p12));
! 2433: }
! 2434: else
! 2435: {
! 2436: p11[2]=lneg(p1); l4=avma;
! 2437: xc=gerepile(l0,l4,gdeuc(xc,p11));
! 2438: }
! 2439: }
! 2440: xd=deriv(xc,v); l2=avma;
! 2441: }
! 2442: k += deg*m;
! 2443: }
! 2444: m++;
! 2445: }
! 2446: while (k!=deg0);
! 2447: }
! 2448: avma=l1;
! 2449: for (j=2; j<=deg0; j++)
! 2450: {
! 2451: p1 = (GEN)y[j];
! 2452: if (gcmp0((GEN)p1[2])) fr=0; else fr=1;
! 2453: for (k=j-1; k>=1; k--)
! 2454: {
! 2455: p2 = (GEN)y[k];
! 2456: if (gcmp0((GEN)p2[2])) f=0; else f=1;
! 2457: if (f<fr) break;
! 2458: if (f==fr && gcmp((GEN)p2[1],(GEN)p1[1]) <= 0) break;
! 2459: y[k+1]=y[k];
! 2460: }
! 2461: y[k+1]=(long)p1;
! 2462: }
! 2463: return y;
! 2464: }
! 2465:
! 2466: GEN
! 2467: roots2(GEN pol,long PREC)
! 2468: {
! 2469: ulong av = avma;
! 2470: long N,flagexactpol,flagrealpol,flagrealrac,ti,i,j;
! 2471: long nbpol,k,av1,multiqol,deg,nbroot,fr,f;
! 2472: GEN p1,p2,rr,EPS,qol,qolbis,x,b,c,*ad,v,tabqol;
! 2473:
! 2474: if (typ(pol)!=t_POL) err(typeer,"roots2");
! 2475: if (!signe(pol)) err(zeropoler,"roots2");
! 2476: N=degpol(pol);
! 2477: if (!N) return cgetg(1,t_COL);
! 2478: if (N==1)
! 2479: {
! 2480: p1=gmul(realun(PREC),(GEN)pol[3]);
! 2481: p2=gneg_i(gdiv((GEN)pol[2],p1));
! 2482: return gerepilecopy(av,p2);
! 2483: }
! 2484: EPS=realun(3); setexpo(EPS, 12 - bit_accuracy(PREC));
! 2485: flagrealpol=1; flagexactpol=1;
! 2486: for (i=2; i<=N+2; i++)
! 2487: {
! 2488: ti=typ(pol[i]);
! 2489: if (ti!=t_INT && ti!=t_INTMOD && !is_frac_t(ti))
! 2490: {
! 2491: flagexactpol=0;
! 2492: if (ti!=t_REAL) flagrealpol=0;
! 2493: }
! 2494: if (ti==t_QUAD)
! 2495: {
! 2496: p1=gmael3(pol,i,1,2);
! 2497: flagrealpol = (gsigne(p1)>0)? 0 : 1;
! 2498: }
! 2499: }
! 2500: rr=cgetg(N+1,t_COL);
! 2501: for (i=1; i<=N; i++)
! 2502: {
! 2503: p1 = cgetc(PREC); rr[i] = (long)p1;
! 2504: for (j=3; j<PREC; j++) mael(p1,2,j)=mael(p1,1,j)=0;
! 2505: }
! 2506: if (flagexactpol) tabqol=square_free_factorization(pol);
! 2507: else
! 2508: {
! 2509: tabqol=cgetg(3,t_MAT);
! 2510: tabqol[1]=lgetg(2,t_COL); mael(tabqol,1,1)=un;
! 2511: tabqol[2]=lgetg(2,t_COL); mael(tabqol,2,1)=lcopy(pol);
! 2512: }
! 2513: nbpol=lg(tabqol[1])-1; nbroot=0;
! 2514: for (k=1; k<=nbpol; k++)
! 2515: {
! 2516: av1=avma; qol=gmael(tabqol,2,k); qolbis=gcopy(qol);
! 2517: multiqol=itos(gmael(tabqol,1,k)); deg=degpol(qol);
! 2518: for (j=deg; j>=1; j--)
! 2519: {
! 2520: x=gzero; flagrealrac=0;
! 2521: if (j==1) x=gneg_i(gdiv((GEN)qolbis[2],(GEN)qolbis[3]));
! 2522: else
! 2523: {
! 2524: x=laguer(qolbis,j,x,EPS,PREC);
! 2525: if (x == NULL) goto RLAB;
! 2526: }
! 2527: if (flagexactpol)
! 2528: {
! 2529: x=gprec(x,(long)((PREC-1)*pariK));
! 2530: x=laguer(qol,deg,x,gmul2n(EPS,-32),PREC+1);
! 2531: }
! 2532: else x=laguer(qol,deg,x,EPS,PREC);
! 2533: if (x == NULL) goto RLAB;
! 2534:
! 2535: if (typ(x)==t_COMPLEX &&
! 2536: gcmp(gabs(gimag(x),PREC),gmul2n(gmul(EPS,gabs(greal(x),PREC)),1))<=0)
! 2537: { x[2]=zero; flagrealrac=1; }
! 2538: else if (j==1 && flagrealpol)
! 2539: { x[2]=zero; flagrealrac=1; }
! 2540: else if (typ(x)!=t_COMPLEX) flagrealrac=1;
! 2541:
! 2542: for (i=1; i<=multiqol; i++) gaffect(x,(GEN)rr[nbroot+i]);
! 2543: nbroot+=multiqol;
! 2544: if (!flagrealpol || flagrealrac)
! 2545: {
! 2546: ad = (GEN*) new_chunk(j+1);
! 2547: for (i=0; i<=j; i++) ad[i]=(GEN)qolbis[i+2];
! 2548: b=(GEN)ad[j];
! 2549: for (i=j-1; i>=0; i--)
! 2550: {
! 2551: c=(GEN)ad[i]; ad[i]=b;
! 2552: b=gadd(gmul((GEN)rr[nbroot],b),c);
! 2553: }
! 2554: v=cgetg(j+1,t_VEC); for (i=1; i<=j; i++) v[i]=(long)ad[j-i];
! 2555: qolbis=gtopoly(v,varn(qolbis));
! 2556: if (flagrealpol)
! 2557: for (i=2; i<=j+1; i++)
! 2558: if (typ(qolbis[i])==t_COMPLEX) mael(qolbis,i,2)=zero;
! 2559: }
! 2560: else
! 2561: {
! 2562: ad = (GEN*) new_chunk(j-1); ad[j-2]=(GEN)qolbis[j+2];
! 2563: p1=gmulsg(2,greal((GEN)rr[nbroot])); p2=gnorm((GEN)rr[nbroot]);
! 2564: ad[j-3]=gadd((GEN)qolbis[j+1],gmul(p1,ad[j-2]));
! 2565: for (i=j-2; i>=2; i--)
! 2566: ad[i-2] = gadd((GEN)qolbis[i+2],gsub(gmul(p1,ad[i-1]),gmul(p2,ad[i])));
! 2567: v=cgetg(j,t_VEC); for (i=1; i<=j-1; i++) v[i]=(long)ad[j-1-i];
! 2568: qolbis=gtopoly(v,varn(qolbis));
! 2569: for (i=2; i<=j; i++)
! 2570: if (typ(qolbis[i])==t_COMPLEX) mael(qolbis,i,2)=zero;
! 2571: for (i=1; i<=multiqol; i++)
! 2572: gaffect(gconj((GEN)rr[nbroot]), (GEN)rr[nbroot+i]);
! 2573: nbroot+=multiqol; j--;
! 2574: }
! 2575: }
! 2576: avma=av1;
! 2577: }
! 2578: for (j=2; j<=N; j++)
! 2579: {
! 2580: x=(GEN)rr[j]; if (gcmp0((GEN)x[2])) fr=0; else fr=1;
! 2581: for (i=j-1; i>=1; i--)
! 2582: {
! 2583: if (gcmp0(gmael(rr,i,2))) f=0; else f=1;
! 2584: if (f<fr) break;
! 2585: if (f==fr && gcmp(greal((GEN)rr[i]),greal(x)) <= 0) break;
! 2586: rr[i+1]=rr[i];
! 2587: }
! 2588: rr[i+1]=(long)x;
! 2589: }
! 2590: return gerepilecopy(av,rr);
! 2591:
! 2592: RLAB:
! 2593: avma = av;
! 2594: for(i=2;i<=N+2;i++)
! 2595: {
! 2596: ti = typ(pol[i]);
! 2597: if (!is_intreal_t(ti)) err(talker,"too many iterations in roots");
! 2598: }
! 2599: if (DEBUGLEVEL)
! 2600: {
! 2601: fprintferr("too many iterations in roots2() ( laguer() ): \n");
! 2602: fprintferr(" real coefficients polynomial, using zrhqr()\n");
! 2603: flusherr();
! 2604: }
! 2605: return zrhqr(pol,PREC);
! 2606: }
! 2607:
! 2608: #define MR 8
! 2609: #define MT 10
! 2610:
! 2611: static GEN
! 2612: laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)
! 2613: {
! 2614: long av = avma, av1,MAXIT,iter,i,j;
! 2615: GEN rac,erre,I,x,abx,abp,abm,dx,x1,b,d,f,g,h,sq,gp,gm,g2,*ffrac;
! 2616:
! 2617: MAXIT=MR*MT; rac=cgetg(3,t_COMPLEX);
! 2618: rac[1]=lgetr(PREC); rac[2]=lgetr(PREC);
! 2619: av1 = avma;
! 2620: I=cgetg(3,t_COMPLEX); I[1]=un; I[2]=un;
! 2621: ffrac=(GEN*)new_chunk(MR+1); for (i=0; i<=MR; i++) ffrac[i]=cgetr(PREC);
! 2622: affrr(dbltor(0.0), ffrac[0]); affrr(dbltor(0.5), ffrac[1]);
! 2623: affrr(dbltor(0.25),ffrac[2]); affrr(dbltor(0.75),ffrac[3]);
! 2624: affrr(dbltor(0.13),ffrac[4]); affrr(dbltor(0.38),ffrac[5]);
! 2625: affrr(dbltor(0.62),ffrac[6]); affrr(dbltor(0.88),ffrac[7]);
! 2626: affrr(dbltor(1.0),ffrac[8]);
! 2627: x=y0;
! 2628: for (iter=1; iter<=MAXIT; iter++)
! 2629: {
! 2630: b=(GEN)pol[N+2]; erre=QuickNormL1(b,PREC);
! 2631: d=gzero; f=gzero; abx=QuickNormL1(x,PREC);
! 2632: for (j=N-1; j>=0; j--)
! 2633: {
! 2634: f=gadd(gmul(x,f),d); d=gadd(gmul(x,d),b);
! 2635: b=gadd(gmul(x,b),(GEN)pol[j+2]);
! 2636: erre=gadd(QuickNormL1(b,PREC),gmul(abx,erre));
! 2637: }
! 2638: erre=gmul(erre,EPS);
! 2639: if (gcmp(QuickNormL1(b,PREC),erre)<=0)
! 2640: {
! 2641: gaffect(x,rac); avma = av1; return rac;
! 2642: }
! 2643: g=gdiv(d,b); g2=gsqr(g); h=gsub(g2, gmul2n(gdiv(f,b),1));
! 2644: sq=gsqrt(gmulsg(N-1,gsub(gmulsg(N,h),g2)),PREC);
! 2645: gp=gadd(g,sq); gm=gsub(g,sq); abp=gnorm(gp); abm=gnorm(gm);
! 2646: if (gcmp(abp,abm)<0) gp=gcopy(gm);
! 2647: if (gsigne(gmax(abp,abm))==1)
! 2648: dx = gdivsg(N,gp);
! 2649: else
! 2650: dx = gmul(gadd(gun,abx),gexp(gmulgs(I,iter),PREC));
! 2651: x1=gsub(x,dx);
! 2652: if (gcmp(QuickNormL1(gsub(x,x1),PREC),EPS)<0)
! 2653: {
! 2654: gaffect(x,rac); avma = av1; return rac;
! 2655: }
! 2656: if (iter%MT) x=gcopy(x1); else x=gsub(x,gmul(ffrac[iter/MT],dx));
! 2657: }
! 2658: avma=av; return NULL;
! 2659: }
! 2660:
! 2661: #undef MR
! 2662: #undef MT
! 2663:
! 2664: /***********************************************************************/
! 2665: /** **/
! 2666: /** ROOTS of a polynomial with REAL coeffs **/
! 2667: /** **/
! 2668: /***********************************************************************/
! 2669: #define RADIX 1L
! 2670: #define COF 0.95
! 2671:
! 2672: /* ONLY FOR REAL COEFFICIENTS MATRIX : replace the matrix x with
! 2673: a symmetric matrix a with the same eigenvalues */
! 2674: static GEN
! 2675: balanc(GEN x)
! 2676: {
! 2677: ulong av = avma;
! 2678: long last,i,j, sqrdx = (RADIX<<1), n = lg(x);
! 2679: GEN r,c,cofgen,a;
! 2680:
! 2681: a = dummycopy(x);
! 2682: last = 0; cofgen = dbltor(COF);
! 2683: while (!last)
! 2684: {
! 2685: last = 1;
! 2686: for (i=1; i<n; i++)
! 2687: {
! 2688: r = c = gzero;
! 2689: for (j=1; j<n; j++)
! 2690: if (j!=i)
! 2691: {
! 2692: c = gadd(c, gabs(gcoeff(a,j,i),0));
! 2693: r = gadd(r, gabs(gcoeff(a,i,j),0));
! 2694: }
! 2695: if (!gcmp0(r) && !gcmp0(c))
! 2696: {
! 2697: GEN g, s = gmul(cofgen, gadd(c,r));
! 2698: long ex = 0;
! 2699: g = gmul2n(r,-RADIX); while (gcmp(c,g) < 0) {ex++; c=gmul2n(c, sqrdx);}
! 2700: g = gmul2n(r, RADIX); while (gcmp(c,g) > 0) {ex--; c=gmul2n(c,-sqrdx);}
! 2701: if (gcmp(gadd(c,r), gmul2n(s,ex)) < 0)
! 2702: {
! 2703: last = 0;
! 2704: for (j=1; j<n; j++) coeff(a,i,j)=lmul2n(gcoeff(a,i,j),-ex);
! 2705: for (j=1; j<n; j++) coeff(a,j,i)=lmul2n(gcoeff(a,j,i), ex);
! 2706: }
! 2707: }
! 2708: }
! 2709: }
! 2710: return gerepilecopy(av, a);
! 2711: }
! 2712:
! 2713: #define SIGN(a,b) ((b)>=0.0 ? fabs(a) : -fabs(a))
! 2714: static GEN
! 2715: hqr(GEN mat) /* find all the eigenvalues of the matrix mat */
! 2716: {
! 2717: long nn,n,m,l,k,j,its,i,mmin,flj,flk;
! 2718: double **a,p,q,r,s,t,u,v,w,x,y,z,anorm,*wr,*wi;
! 2719: const double eps = 0.000001;
! 2720: GEN eig;
! 2721:
! 2722: n=lg(mat)-1; a=(double**)gpmalloc(sizeof(double*)*(n+1));
! 2723: for (i=1; i<=n; i++) a[i]=(double*)gpmalloc(sizeof(double)*(n+1));
! 2724: for (j=1; j<=n; j++)
! 2725: for (i=1; i<=n; i++) a[i][j]=gtodouble((GEN)((GEN)mat[j])[i]);
! 2726: wr=(double*)gpmalloc(sizeof(double)*(n+1));
! 2727: wi=(double*)gpmalloc(sizeof(double)*(n+1));
! 2728:
! 2729: anorm=fabs(a[1][1]);
! 2730: for (i=2; i<=n; i++) for (j=(i-1); j<=n; j++) anorm+=fabs(a[i][j]);
! 2731: nn=n; t=0.0;
! 2732: if (DEBUGLEVEL>3) { fprintferr("* Finding eigenvalues\n"); flusherr(); }
! 2733: while (nn>=1)
! 2734: {
! 2735: its=0;
! 2736: do
! 2737: {
! 2738: for (l=nn; l>=2; l--)
! 2739: {
! 2740: s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.0) s=anorm;
! 2741: if ((double)(fabs(a[l][l-1])+s)==s) break;
! 2742: }
! 2743: x=a[nn][nn];
! 2744: if (l==nn){ wr[nn]=x+t; wi[nn--]=0.0; }
! 2745: else
! 2746: {
! 2747: y=a[nn-1][nn-1];
! 2748: w=a[nn][nn-1]*a[nn-1][nn];
! 2749: if (l == nn-1)
! 2750: {
! 2751: p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x+=t;
! 2752: if (q>=0.0 || fabs(q)<=eps)
! 2753: {
! 2754: z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z;
! 2755: if (fabs(z)>eps) wr[nn]=x-w/z;
! 2756: wi[nn-1]=wi[nn]=0.0;
! 2757: }
! 2758: else{ wr[nn-1]=wr[nn]=x+p; wi[nn-1]=-(wi[nn]=z); }
! 2759: nn-=2;
! 2760: }
! 2761: else
! 2762: {
! 2763: p = q = r = 0.0; /* for lint */
! 2764: if (its==30) err(talker,"too many iterations in hqr");
! 2765: if (its==10 || its==20)
! 2766: {
! 2767: t+=x; for (i=1; i<=nn; i++) a[i][i]-=x;
! 2768: s = fabs(a[nn][nn-1]) + fabs(a[nn-1][nn-2]);
! 2769: y=x=0.75*s; w=-0.4375*s*s;
! 2770: }
! 2771: its++;
! 2772: for (m=nn-2; m>=l; m--)
! 2773: {
! 2774: z=a[m][m]; r=x-z; s=y-z;
! 2775: p=(r*s-w)/a[m+1][m]+a[m][m+1];
! 2776: q=a[m+1][m+1]-z-r-s;
! 2777: r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); p/=s; q/=s; r/=s;
! 2778: if (m==l) break;
! 2779: u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
! 2780: v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
! 2781: if ((double)(u+v)==v) break;
! 2782: }
! 2783: for (i=m+2; i<=nn; i++){ a[i][i-2]=0.0; if (i!=(m+2)) a[i][i-3]=0.0; }
! 2784: for (k=m; k<=nn-1; k++)
! 2785: {
! 2786: if (k!=m)
! 2787: {
! 2788: p=a[k][k-1]; q=a[k+1][k-1];
! 2789: r = (k != nn-1)? a[k+2][k-1]: 0.0;
! 2790: x = fabs(p)+fabs(q)+fabs(r);
! 2791: if (x != 0.0) { p/=x; q/=x; r/=x; }
! 2792: }
! 2793: s = SIGN(sqrt(p*p+q*q+r*r),p);
! 2794: if (s == 0.0) continue;
! 2795:
! 2796: if (k==m)
! 2797: { if (l!=m) a[k][k-1] = -a[k][k-1]; }
! 2798: else
! 2799: a[k][k-1] = -s*x;
! 2800: p+=s; x=p/s; y=q/s; z=r/s; q/=p; r/=p;
! 2801: for (j=k; j<=nn; j++)
! 2802: {
! 2803: p = a[k][j]+q*a[k+1][j];
! 2804: if (k != nn-1) { p+=r*a[k+2][j]; a[k+2][j]-=p*z; }
! 2805: a[k+1][j] -= p*y; a[k][j] -= p*x;
! 2806: }
! 2807: mmin = (nn < k+3)? nn: k+3;
! 2808: for (i=l; i<=mmin; i++)
! 2809: {
! 2810: p = x*a[i][k]+y*a[i][k+1];
! 2811: if (k != nn-1) { p+=z*a[i][k+2]; a[i][k+2]-=p*r; }
! 2812: a[i][k+1] -= p*q; a[i][k] -= p;
! 2813: }
! 2814: }
! 2815: }
! 2816: }
! 2817: }
! 2818: while (l<nn-1);
! 2819: }
! 2820: for (j=2; j<=n; j++) /* ordering the roots */
! 2821: {
! 2822: x=wr[j]; y=wi[j]; if (y) flj=1; else flj=0;
! 2823: for (k=j-1; k>=1; k--)
! 2824: {
! 2825: if (wi[k]) flk=1; else flk=0;
! 2826: if (flk<flj) break;
! 2827: if (!flk && !flj && wr[k]<=x) break;
! 2828: if (flk&&flj&& wr[k]<x) break;
! 2829: if (flk&&flj&& wr[k]==x && wi[k]>0) break;
! 2830: wr[k+1]=wr[k]; wi[k+1]=wi[k];
! 2831: }
! 2832: wr[k+1]=x; wi[k+1]=y;
! 2833: }
! 2834: if (DEBUGLEVEL>3) { fprintferr("* Eigenvalues computed\n"); flusherr(); }
! 2835: for (i=1; i<=n; i++) free(a[i]); free(a); eig=cgetg(n+1,t_COL);
! 2836: for (i=1; i<=n; i++)
! 2837: {
! 2838: if (wi[i])
! 2839: {
! 2840: GEN p1 = cgetg(3,t_COMPLEX);
! 2841: eig[i]=(long)p1;
! 2842: p1[1]=(long)dbltor(wr[i]);
! 2843: p1[2]=(long)dbltor(wi[i]);
! 2844: }
! 2845: else eig[i]=(long)dbltor(wr[i]);
! 2846: }
! 2847: free(wr); free(wi); return eig;
! 2848: }
! 2849:
! 2850: /* ONLY FOR POLYNOMIAL WITH REAL COEFFICIENTS : give the roots of the
! 2851: * polynomial a (first, the real roots, then the non real roots) in
! 2852: * increasing order of their real parts MULTIPLE ROOTS ARE FORBIDDEN.
! 2853: */
! 2854: GEN
! 2855: zrhqr(GEN a,long prec)
! 2856: {
! 2857: ulong av = avma;
! 2858: long i,j,prec2, n = degpol(a), ex = -bit_accuracy(prec);
! 2859: GEN aa,b,p1,rt,rr,hess,x,dx,y,newval,oldval;
! 2860:
! 2861: hess = cgetg(n+1,t_MAT);
! 2862: for (j=1; j<=n; j++)
! 2863: {
! 2864: p1 = cgetg(n+1,t_COL); hess[j] = (long)p1;
! 2865: p1[1] = lneg(gdiv((GEN)a[n-j+2],(GEN)a[n+2]));
! 2866: for (i=2; i<=n; i++) p1[i] = (i==(j+1))? un: zero;
! 2867: }
! 2868: rt = hqr(balanc(hess));
! 2869: prec2 = 2*prec; /* polishing the roots */
! 2870: aa = gprec_w(a, prec2);
! 2871: b = derivpol(aa); rr = cgetg(n+1,t_COL);
! 2872: for (i=1; i<=n; i++)
! 2873: {
! 2874: x = gprec_w((GEN)rt[i], prec2);
! 2875: for (oldval=NULL;; oldval=newval, x=y)
! 2876: { /* Newton iteration */
! 2877: dx = poleval(b,x);
! 2878: if (gexpo(dx) < ex)
! 2879: err(talker,"polynomial has probably multiple roots in zrhqr");
! 2880: y = gsub(x, gdiv(poleval(aa,x),dx));
! 2881: newval = gabs(poleval(aa,y),prec2);
! 2882: if (gexpo(newval) < ex || (oldval && gcmp(newval,oldval) > 0)) break;
! 2883: }
! 2884: if (DEBUGLEVEL>3) fprintferr("%ld ",i);
! 2885: rr[i] = (long)cgetc(prec); gaffect(y, (GEN)rr[i]);
! 2886: }
! 2887: if (DEBUGLEVEL>3) { fprintferr("\npolished roots = %Z",rr); flusherr(); }
! 2888: return gerepilecopy(av, rr);
! 2889: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to polarit1.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 / basemath