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