Annotation of OpenXM_contrib/pari/src/basemath/polarit2.c, Revision 1.1
1.1 ! maekawa 1: /***********************************************************************/
! 2: /** **/
! 3: /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
! 4: /** (second part) **/
! 5: /** **/
! 6: /***********************************************************************/
! 7: /* $Id: polarit2.c,v 1.2 1999/09/20 13:30:05 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: GEN
! 11: polsym(GEN x, long n)
! 12: {
! 13: long av1,av2,dx=lgef(x)-3,i,k;
! 14: GEN s,y,x_lead;
! 15:
! 16: if (n<0) err(impl,"polsym of a negative n");
! 17: if (typ(x) != t_POL) err(typeer,"polsym");
! 18: if (!signe(x)) err(zeropoler,"polsym");
! 19: y=cgetg(n+2,t_COL); y[1]=lstoi(dx);
! 20: x_lead=(GEN)x[dx+2]; if (gcmp1(x_lead)) x_lead=NULL;
! 21: for (k=1; k<=n; k++)
! 22: {
! 23: av1=avma; s = (dx>=k)? gmulsg(k,(GEN)x[dx+2-k]): gzero;
! 24: for (i=1; i<k && i<=dx; i++)
! 25: s = gadd(s,gmul((GEN)y[k-i+1],(GEN)x[dx+2-i]));
! 26: if (x_lead) s = gdiv(s,x_lead);
! 27: av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
! 28: }
! 29: return y;
! 30: }
! 31:
! 32: static int (*polcmp_coeff_cmp)(GEN,GEN);
! 33:
! 34: /* assume x and y are polynomials in the same variable whose coeffs can be
! 35: * compared (used to sort polynomial factorizations)
! 36: */
! 37: static int
! 38: polcmp(GEN x, GEN y)
! 39: {
! 40: long i, lx = lgef(x), ly = lgef(y);
! 41: int fl;
! 42: #if 0 /* for efficiency */
! 43: if (typ(x) != t_POL || typ(y) != t_POL)
! 44: err(talker,"not a polynomial in polcmp");
! 45: #endif
! 46: if (lx > ly) return 1;
! 47: if (lx < ly) return -1;
! 48: for (i=lx-1; i>1; i--)
! 49: if ((fl = polcmp_coeff_cmp((GEN)x[i], (GEN)y[i]))) return fl;
! 50: return 0;
! 51: }
! 52:
! 53: /* sort factorization so that factors appear by decreasing degree, in place */
! 54: GEN
! 55: sort_factor(GEN y, int (*cmp)(GEN,GEN))
! 56: {
! 57: GEN w,c1,c2, p1 = (GEN)y[1], p2 = (GEN)y[2];
! 58: long av=avma, n=lg(p1), i;
! 59: int (*old)(GEN,GEN) = polcmp_coeff_cmp;
! 60:
! 61: c1 = new_chunk(n); c2 = new_chunk(n);
! 62: polcmp_coeff_cmp = cmp;
! 63: w = gen_sort(p1, cmp_IND | cmp_C, polcmp);
! 64: for (i=1; i<n; i++) { c1[i]=p1[w[i]]; c2[i]=p2[w[i]]; }
! 65: for (i=1; i<n; i++) { p1[i]=c1[i]; p2[i]=c2[i]; }
! 66: polcmp_coeff_cmp = old; avma = av; return y;
! 67: }
! 68:
! 69: /* for internal use */
! 70: GEN
! 71: centermod(GEN x, GEN p)
! 72: {
! 73: long av,i,lx;
! 74: GEN y,p1,ps2;
! 75:
! 76: ps2=shifti(p,-1);
! 77: switch(typ(x))
! 78: {
! 79: case t_INT:
! 80: y=modii(x,p);
! 81: if (cmpii(y,ps2)>0) return subii(y,p);
! 82: return y;
! 83:
! 84: case t_POL: lx=lgef(x);
! 85: y=cgetg(lx,t_POL); y[1]=x[1];
! 86: for (i=2; i<lx; i++)
! 87: {
! 88: av=avma; p1=modii((GEN)x[i],p);
! 89: if (cmpii(p1,ps2)>0) p1=subii(p1,p);
! 90: y[i]=lpileupto(av,p1);
! 91: }
! 92: return y;
! 93:
! 94: case t_COL: lx=lg(x);
! 95: y=cgetg(lx,t_COL);
! 96: for (i=1; i<lx; i++)
! 97: {
! 98: p1=modii((GEN)x[i],p);
! 99: if (cmpii(p1,ps2)>0) p1=subii(p1,p);
! 100: y[i]=(long)p1;
! 101: }
! 102: return y;
! 103: }
! 104: return x;
! 105: }
! 106:
! 107: static GEN
! 108: decpol(GEN x, long klim)
! 109: {
! 110: short int pos[200];
! 111: long av=avma,av1,tete,k,kin,i,j,i1,i2,fl,d,nbfact;
! 112: GEN res,p1,p2;
! 113:
! 114: kin=1; res=cgetg(lgef(x)-2,t_VEC); nbfact=0;
! 115: p1=roots(x,DEFAULTPREC); d=lg(p1)-1; if (!klim) klim=d;
! 116: do
! 117: {
! 118: fl=1;
! 119: for (k=kin; k+k <= d && k<=klim; k++)
! 120: {
! 121: for (i=0; i<=k; i++) pos[i]=i;
! 122: do
! 123: {
! 124: av1=avma; p2=gzero; j=0;
! 125: for (i1=1; i1<=k; i1++) p2=gadd(p2,(GEN)p1[pos[i1]]);
! 126: if (gexpo(gimag(p2))<-20 && gexpo(gsub(p2,ground(p2)))<-20)
! 127: {
! 128: p2=gun;
! 129: for (i1=1; i1<=k; i1++)
! 130: p2=gmul(p2,gsub(polx[0],(GEN)p1[pos[i1]]));
! 131: p2 = ground(p2);
! 132: if (gcmp0(gimag(p2)) && gcmp0(gmod(x,p2)))
! 133: {
! 134: res[++nbfact]=(long)p2; x=gdiv(x,p2);
! 135: kin=k; p2=cgetg(d-k+1,t_COL);
! 136: for (i=i1=i2=1; i<=d; i++)
! 137: {
! 138: if (i1<=k && i==pos[i1]) i1++;
! 139: else p2[i2++]=p1[i];
! 140: }
! 141: p1=p2; d-=k; fl=0; break;
! 142: }
! 143: }
! 144: avma=av1; pos[k]++;
! 145: while (pos[k-j] > d-j) { j++; pos[k-j]++; }
! 146: for (i=k-j+1; i<=k; i++) pos[i]=i+pos[k-j]-k+j;
! 147: }
! 148: while (j<k);
! 149: if (!fl) break;
! 150: }
! 151: if (lgef(x)<=3) break;
! 152: }
! 153: while (!fl || (k+k <= d && k<=klim));
! 154: if (lgef(x)>3) res[++nbfact]=(long)x;
! 155: setlg(res,nbfact+1); tete=avma;
! 156: return gerepile(av,tete,greal(res));
! 157: }
! 158:
! 159: /* This code originally by Richard Schroeppel */
! 160:
! 161: /* Note that PARI's idea of the maximum possible coefficient involves the
! 162: * limit on the degree (klim). Consider revising this. If I don't respect
! 163: * the degree limit when testing potential factors, there's the possibility
! 164: * that I might identify a high degree factor that isn't irreducible, because
! 165: * it's lower degree divisors were missed because they had a coefficient
! 166: * outside the Borne limit for klim, but the higher degree factor had it's
! 167: * coefficients within Borne. This would still have the property that any
! 168: * factors of degree <= klim were guaranteed irr, but higher degrees
! 169: * (> 2*klim) might not be irr.
! 170: *
! 171: * The subroutine:
! 172: * fxn points at the first unconsidered factor for the current combination
! 173: * psf is the product-so-far, or 0 for a null product
! 174: * dlim is the degree limit remaining for unconsidered divisors
! 175: * other arguments are "global" and must already be setup
! 176: * as factors are found, they are put in cmbf_fax; the count is kept in
! 177: * cmbf_faxn; and they are divided out of cmbf_target; the degree and
! 178: * leading coefficient are updated; and the constituent modular factors
! 179: * are deleted from cmbf_modfax.
! 180: * exit value is 1 if any factors are found.
! 181: * If psf is 0, all factors made up from pieces at or after fxn will be
! 182: * found & removed. If psf is not 0, only the factor which is a
! 183: * continuation of psf will be found.
! 184: */
! 185: /* setup for calling cmbf = combine_factors */
! 186: static GEN cmbf_target; /* target poly. Assume content removed */
! 187: static GEN cmbf_lc; /* leading coefficient */
! 188: static GEN cmbf_abslc; /* |lc| */
! 189: static GEN cmbf_abslcxtarget;/* abslc * target */
! 190: static GEN cmbf_mod; /* modulus */
! 191: static GEN cmbf_fax; /* result array + one cell for leftover cst. */
! 192: static GEN cmbf_modfax; /* array of modular factors. Each has LC 1.1 based
! 193: indexing. Product should be congruent to a/lc(a). */
! 194: static long cmbf_degree; /* degree(target) */
! 195: static long cmbf_modfaxn; /* number of modular factors */
! 196: static long cmbf_faxn; /* last used cell in fax. # of factors found */
! 197: static GEN cmbf_modhalf;
! 198:
! 199: /* assume same variables, y normalized, non constant */
! 200: static GEN
! 201: polidivis(GEN x, GEN y, GEN bound)
! 202: {
! 203: long dx,dy,dz,i,j,av, vx = varn(x);
! 204: GEN z,p1,y_lead;
! 205:
! 206: dy=lgef(y)-3;
! 207: dx=lgef(x)-3;
! 208: dz=dx-dy; if (dz<0) return NULL;
! 209: z=cgetg(dz+3,t_POL);
! 210: x += 2; y += 2; z += 2;
! 211: y_lead = (GEN)y[dy];
! 212: if (gcmp1(y_lead)) y_lead = NULL;
! 213:
! 214: p1 = (GEN)x[dx];
! 215: z[dz]=y_lead? ldivii(p1,y_lead): licopy(p1);
! 216: for (i=dx-1; i>=dy; i--)
! 217: {
! 218: av=avma; p1=(GEN)x[i];
! 219: for (j=i-dy+1; j<=i && j<=dz; j++)
! 220: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 221: if (y_lead) { p1 = gdiv(p1,y_lead); if (typ(p1)!=t_INT) return NULL; }
! 222: if (absi_cmp(p1, bound) > 0) return NULL;
! 223: p1 = gerepileupto(av,p1);
! 224: z[i-dy] = (long)p1;
! 225: }
! 226: av = avma;
! 227: for (;; i--)
! 228: {
! 229: p1 = (GEN)x[i];
! 230: /* we always enter this loop at least once */
! 231: for (j=0; j<=i && j<=dz; j++)
! 232: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 233: if (!gcmp0(p1)) return NULL;
! 234: avma = av;
! 235: if (!i) break;
! 236: }
! 237: z -= 2;
! 238: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 239: return z;
! 240: }
! 241:
! 242: static int
! 243: try_div(GEN x, GEN y, GEN *q)
! 244: {
! 245: if (resii((GEN)x[2], (GEN)y[2]) != gzero)
! 246: {
! 247: if (DEBUGLEVEL>6) fprintferr(".");
! 248: return 0;
! 249: }
! 250: *q = polidivis(x,y,cmbf_modhalf);
! 251: if (*q) return 1;
! 252: if (DEBUGLEVEL>6) fprintferr("*");
! 253: return 0;
! 254: }
! 255:
! 256: static int
! 257: cmbf(long fxn, GEN psf, long dlim, long hint)
! 258: {
! 259: long newd,av,val2,val=0; /* assume failure */
! 260: GEN newf,newpsf,quo,cont;
! 261:
! 262: if (dlim <= 0 || fxn > cmbf_modfaxn) return 0;
! 263:
! 264: /* first, try deeper factors without considering the current one */
! 265: if (fxn < cmbf_modfaxn)
! 266: {
! 267: val = cmbf(fxn+1,psf,dlim,hint);
! 268: if (val && psf) return val;
! 269: }
! 270:
! 271: /* second, try including the current modular factor in the product */
! 272: newf = (GEN)cmbf_modfax[fxn];
! 273: if (!newf) return val; /* modular factor already used */
! 274:
! 275: newd = lgef(newf)-3; if (newd > dlim) return val;
! 276: av = avma;
! 277: if (!(newd%hint))
! 278: {
! 279: if (!psf) psf = cmbf_abslc;
! 280: newpsf = centermod(gmul(psf,newf),cmbf_mod);
! 281:
! 282: /* try out the new combination */
! 283: if (try_div(cmbf_abslcxtarget, newpsf, &quo))
! 284: {
! 285: /* found a factor */
! 286: cont = content(newpsf);
! 287: if (signe(leading_term(newpsf)) < 0) cont = negi(cont);
! 288: newpsf = gdiv(newpsf,cont);
! 289: cmbf_fax[++cmbf_faxn] = (long)newpsf; /* store factor */
! 290: cmbf_modfax[fxn] = 0; /* remove used modular factor */
! 291: if (DEBUGLEVEL > 2)
! 292: { fprintferr("\n"); msgtimer("to find factor %Z",newpsf); }
! 293: /* fix up target */
! 294: cmbf_target = gdiv(quo,leading_term(newpsf));
! 295: cmbf_degree = lgef(cmbf_target)-3;
! 296: cmbf_lc = (GEN)cmbf_target[cmbf_degree+2];
! 297: cmbf_abslc = absi(cmbf_lc);
! 298: cmbf_abslcxtarget = gmul(cmbf_abslc,cmbf_target);
! 299: return 1;
! 300: }
! 301: }
! 302: /* degree too large, or no more modular factors? */
! 303: if (newd == dlim || fxn == cmbf_modfaxn) return val;
! 304: /* include more factors */
! 305: val2 = cmbf(fxn+1,newpsf,dlim-newd,hint);
! 306: if (val2) cmbf_modfax[fxn] = 0; /* remove used modular factor */
! 307: else avma = av;
! 308: return val || val2;
! 309: }
! 310:
! 311: static long
! 312: min_deg(long jmax,ulong tabbit[])
! 313: {
! 314: long j,k=1,j1=2;
! 315:
! 316: for (j=jmax; j>=0; j--)
! 317: {
! 318: for ( ; k<=15; k++)
! 319: {
! 320: if (tabbit[j]&j1) return ((jmax-j)<<4)+k;
! 321: j1<<=1;
! 322: }
! 323: k=0; j1=1;
! 324: }
! 325: return (jmax<<4)+15;
! 326: }
! 327:
! 328: static void
! 329: dotab(long d,long jmax,ulong *tabkbit,ulong *tmp)
! 330: {
! 331: ulong rem,pro;
! 332: long j,d1,d2;
! 333:
! 334: d1=d>>4; d2=d&15; rem=0;
! 335: for (j=jmax-d1; j>=0; j--)
! 336: {
! 337: pro = tabkbit[j+d1]<<d2;
! 338: tmp[j] = (pro&0xffff)+rem; rem=pro>>16;
! 339: }
! 340: for (j=jmax-d1; j>=0; j--) tabkbit[j] |= tmp[j];
! 341: }
! 342:
! 343: /* lift factorisation mod p to mod p^e = pev */
! 344: GEN
! 345: hensel_lift_fact(GEN pol, GEN fact, GEN p, GEN pev, long e)
! 346: {
! 347: long ev,i, nf = lg(fact);
! 348: GEN aprov = pol, res = cgetg(nf, t_VEC);
! 349:
! 350: for (i=nf-1; i; i--)
! 351: {
! 352: GEN ae,be,u,v,ae2,be2,s,t,pe,pe2,z,g;
! 353: long ltop = avma;
! 354:
! 355: ae = (GEN)fact[i]; /* lead coeff(ae) = 1 */
! 356: be = Fp_poldivres(aprov,ae,p, NULL);
! 357: g = (GEN)Fp_pol_extgcd(ae,be,p,&u,&v)[2]; /* deg g = 0 */
! 358: if (!gcmp1(g))
! 359: {
! 360: g = mpinvmod(g, p);
! 361: u = gmul(u, g);
! 362: v = gmul(v, g);
! 363: }
! 364: for(pe=p,ev=1;;)
! 365: {
! 366: ev<<=1; pe2=(ev>=e)? pev: sqri(pe);
! 367: g = gadd(aprov, gneg_i(gmul(ae,be)));
! 368:
! 369: g = Fp_pol_red(g, pe2); g = gdivexact(g, pe);
! 370: z = Fp_pol_red(gmul(v,g), pe);
! 371: t = Fp_poldivres(z,ae,pe, &s);
! 372: t = gadd(gmul(u,g), gmul(t,be));
! 373: t = Fp_pol_red(t, pe);
! 374: be2 = gadd(be, gmul(t,pe));
! 375: ae2 = gadd(ae, gmul(s,pe)); /* already reduced mod pe2 */
! 376: if (ev>=e) break;
! 377:
! 378: g = gadd(gun, gneg_i(gadd(gmul(u,ae2),gmul(v,be2))));
! 379:
! 380: g = Fp_pol_red(g, pe2); g = gdivexact(g, pe);
! 381: z = Fp_pol_red(gmul(v,g), pe);
! 382: t = Fp_poldivres(z,ae,pe, &s);
! 383: t = gadd(gmul(u,g), gmul(t,be));
! 384: t = Fp_pol_red(t, pe);
! 385: u = gadd(u, gmul(t,pe));
! 386: v = gadd(v, gmul(s,pe));
! 387: pe=pe2; ae=ae2; be=be2;
! 388: }
! 389: ae = gerepileupto(ltop, ae2);
! 390: aprov = Fp_poldivres(aprov,ae,pev, NULL);
! 391: ltop = avma;
! 392: res[i] = (long)ae;
! 393: if (DEBUGLEVEL > 4)
! 394: fprintferr("...lifting factor of degree %3ld. Time = %ld\n",
! 395: lgef(ae)-3,timer2());
! 396: }
! 397: return res;
! 398: }
! 399:
! 400: long split_berlekamp(GEN Q, GEN *t, GEN pp, GEN pps2);
! 401:
! 402: /* assume degree(a) > 0, a(0) != 0, and a squarefree */
! 403: static GEN
! 404: squff(GEN a, long klim, long hint)
! 405: {
! 406: GEN Q,prime,primes2,factorsmodp,p1,p2,y,g,z,w,pe,*tabd,*tabdnew;
! 407: long av=avma,tetpil,va=varn(a),da=lgef(a)-3;
! 408: long chosenp,p,nfacp,lbit,i,j,d,e,np,nmax,lgg,nf,nft;
! 409: ulong *tabbit, *tabkbit, *tmp;
! 410: byteptr pt=diffptr;
! 411: const int NOMBDEP = 5;
! 412:
! 413: if (hint < 0) return decpol(a,klim);
! 414: if (DEBUGLEVEL > 2) timer2();
! 415: lbit=(da>>4)+1; nmax=da+1; i=da>>1;
! 416: if (!klim || klim>i) klim=i;
! 417: tabdnew = (GEN*)new_chunk(nmax);
! 418: tabbit = (ulong*)new_chunk(lbit);
! 419: tabkbit = (ulong*)new_chunk(lbit);
! 420: tmp = (ulong*)new_chunk(lbit);
! 421: p = np = 0; prime = icopy(gun);
! 422: while (np < NOMBDEP)
! 423: {
! 424: GEN minuspolx;
! 425: p += *pt++; if (!*pt) err(primer1);
! 426: if (!smodis((GEN)a[da+2],p)) continue;
! 427: prime[2] = p; z = Fp_pol(a, prime);
! 428: if (gcmp0(discsr(z))) continue;
! 429: z = lift_intern(z);
! 430:
! 431: for (j=0; j<lbit-1; j++) tabkbit[j]=0;
! 432: tabkbit[j]=1;
! 433: d=0; e=da; nfacp=0;
! 434: w = polx[va]; minuspolx = gneg(w);
! 435: while (d < (e>>1))
! 436: {
! 437: d++; w = Fp_pow_mod_pol(w, prime, z, prime);
! 438: g = Fp_pol_gcd(z, gadd(w, minuspolx), prime);
! 439: tabdnew[d]=g; lgg=lgef(g)-3;
! 440: if (lgg>0)
! 441: {
! 442: z = Fp_poldivres(z, g, prime, NULL);
! 443: e = lgef(z)-3;
! 444: w = Fp_poldivres(w, z, prime, ONLY_REM);
! 445: lgg /= d; nfacp += lgg;
! 446: if (DEBUGLEVEL>5)
! 447: fprintferr(" %3ld %s of degree %3ld\n",
! 448: lgg, lgg==1?"factor": "factors",d);
! 449: for (j=1; j<=lgg; j++) dotab(d,lbit-1,tabkbit,tmp);
! 450: }
! 451: }
! 452: if (e>0)
! 453: {
! 454: if (DEBUGLEVEL>5)
! 455: fprintferr(" %3ld factor of degree %3ld\n",1,e);
! 456: tabdnew[e]=z; nfacp++;
! 457: dotab(e,lbit-1,tabkbit,tmp);
! 458: }
! 459: if (np)
! 460: for (j=0; j<lbit; j++) tabbit[j] &= tabkbit[j];
! 461: else
! 462: for (j=0; j<lbit; j++) tabbit[j] = tabkbit[j];
! 463: if (DEBUGLEVEL > 4)
! 464: fprintferr("...tried prime %3ld (%-3ld %s). Time = %ld\n",
! 465: p,nfacp,nfacp==1?"factor": "factors",timer2());
! 466: if (min_deg(lbit-1,tabbit) > klim)
! 467: {
! 468: avma=av; y=cgetg(2,t_COL); y[1]=lcopy(a); return y;
! 469: }
! 470: if (nfacp<nmax)
! 471: {
! 472: nmax=nfacp; tabd=tabdnew; chosenp=p;
! 473: for (j=d+1; j<e; j++) tabd[j]=polun[va];
! 474: tabdnew = (GEN*)new_chunk(da);
! 475: }
! 476: np++;
! 477: }
! 478: prime[2]=chosenp; primes2 = shifti(prime, -1);
! 479: nf=nmax; nft=1;
! 480: y=cgetg(nf+1,t_COL); factorsmodp=cgetg(nf+1,t_COL);
! 481: Q = cgetg(da+1,t_MAT);
! 482: for (i=1; i<=da; i++) Q[i] = lgetg(da+1, t_COL);
! 483: p1 = (GEN)Q[1]; for (i=1; i<=da; i++) p1[i] = zero;
! 484: for (d=1; nft<=nf; d++)
! 485: {
! 486: g=tabd[d]; lgg=lgef(g)-3;
! 487: if (lgg)
! 488: {
! 489: long n = lgg/d;
! 490: factorsmodp[nft] = (long)normalize_mod_p(g, prime);
! 491: if (n==1) nft++;
! 492: else
! 493: {
! 494: (void)split_berlekamp(Q, (GEN*)(factorsmodp+nft),prime,primes2);
! 495: nft += n;
! 496: }
! 497: }
! 498: }
! 499: if (DEBUGLEVEL > 4) msgtimer("splitting mod p = %ld",chosenp);
! 500:
! 501: p1=gzero; for (i=2; i<=da+2; i++) p1=addii(p1,sqri((GEN)a[i]));
! 502: p2=absi((GEN)a[da+2]);
! 503: p1=addii(p2,addsi(1,racine(p1)));
! 504: p1=mulii(p1,binome(stoi(da-1),klim));
! 505: p1=shifti(mulii(p2,p1),1);
! 506: for (e=1,pe=gun; ; e++)
! 507: {
! 508: pe = mulis(pe,chosenp);
! 509: if (cmpii(pe,p1) >= 0) break;
! 510: }
! 511: if (DEBUGLEVEL > 4) fprintferr("coeff bound: %Z\n",p1);
! 512:
! 513: cmbf_target = a;
! 514: cmbf_degree = lgef(a)-3;
! 515: cmbf_lc = (GEN)a[lgef(a)-1];
! 516: cmbf_abslc = absi(cmbf_lc);
! 517: cmbf_abslcxtarget = gmul(cmbf_abslc,a);
! 518: cmbf_mod = pe;
! 519: cmbf_modhalf = shifti(pe,-1);
! 520: cmbf_modfax = hensel_lift_fact(a,factorsmodp,prime,pe,e);
! 521: cmbf_modfaxn = nf;
! 522: cmbf_fax = cgetg(nf+2,t_COL);
! 523: cmbf_faxn = 0;
! 524: /* sorting factors decreasing by degree helps if klim is used
! 525: * If not, can start with first arg of 2 instead of 1, saving some time.
! 526: * Should be sending tabbit through for more efficiency ???
! 527: */
! 528: cmbf(1,NULL,klim,hint); /* The Call */
! 529: if (cmbf_degree)
! 530: {
! 531: if (signe(cmbf_lc)<0) cmbf_target = gneg_i(cmbf_target);
! 532: cmbf_fax[++cmbf_faxn] = (long)cmbf_target; /* leftover factor */
! 533: }
! 534: if (DEBUGLEVEL>6) fprintferr("\n");
! 535: tetpil=avma; y=cgetg(cmbf_faxn+1,t_COL);
! 536: for (i=1; i<=cmbf_faxn; i++) y[i]=lcopy((GEN)cmbf_fax[i]);
! 537: return gerepile(av,tetpil,y);
! 538: }
! 539:
! 540: /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs
! 541: * de degre <= klim
! 542: */
! 543: GEN
! 544: factpol(GEN x, long klim, long hint)
! 545: {
! 546: GEN fa,p1,d,t,v,w, y = cgetg(3,t_MAT);
! 547: long av=avma,av2,lx,vx,i,j,k,ex,nbfac,zval;
! 548:
! 549: if (typ(x)!=t_POL) err(notpoler,"factpol");
! 550: if (!signe(x)) err(zeropoler,"factpol");
! 551:
! 552: ex = nbfac = 0;
! 553: p1 = x+2; while (gcmp0((GEN)*p1)) p1++;
! 554: zval = p1 - (x + 2);
! 555: lx = lgef(x) - zval;
! 556: vx = varn(x);
! 557: if (zval)
! 558: {
! 559: x = cgetg(lx, t_POL); p1 -= 2;
! 560: for (i=2; i<lx; i++) x[i] = p1[i];
! 561: x[1] = evalsigne(1)|evalvarn(vx)|evallgef(lx);
! 562: nbfac++;
! 563: }
! 564: /* now x(0) != 0 */
! 565: if (lx==3) goto END;
! 566: p1 = cgetg(1,t_VEC); fa=cgetg(lx,t_VEC);
! 567: for (i=1; i<lx; i++) fa[i] = (long)p1;
! 568: d=content(x); if (gsigne(leading_term(x)) < 0) d = gneg_i(d);
! 569: if (!gcmp1(d)) x=gdiv(x,d);
! 570: if (lx==4) { nbfac++; ex++; fa[1] = (long)concatsp(p1,x); goto END; }
! 571:
! 572: w=derivpol(x); t=modulargcd(x,w);
! 573: if (!gcmp1(t)) { x=gdeuc(x,t); w=gdeuc(w,t); }
! 574: k=1;
! 575: while (k)
! 576: {
! 577: ex++; w=gadd(w, gneg_i(derivpol(x))); k=signe(w);
! 578: if (k) { t=modulargcd(x,w); x=gdeuc(x,t); w=gdeuc(w,t); } else t=x;
! 579: if (lgef(t) > 3)
! 580: {
! 581: fa[ex] = (long)squff(t,klim,hint);
! 582: nbfac += lg(fa[ex])-1;
! 583: }
! 584: }
! 585: END: av2=avma;
! 586: v=cgetg(nbfac+1,t_COL); y[1]=(long)v;
! 587: w=cgetg(nbfac+1,t_COL); y[2]=(long)w;
! 588: if (zval) { v[1]=lpolx[vx]; w[1]=lstoi(zval); k=1; } else k=0;
! 589: for (i=1; i<=ex; i++)
! 590: for (j=1; j<lg(fa[i]); j++)
! 591: {
! 592: k++; v[k]=lcopy(gmael(fa,i,j)); w[k]=lstoi(i);
! 593: }
! 594: gerepilemanyvec(av,av2,y+1,2);
! 595: return sort_factor(y, cmpii);
! 596: }
! 597:
! 598: GEN
! 599: factpol2(GEN x, long klim)
! 600: {
! 601: return factpol(x,klim,-1);
! 602: }
! 603:
! 604: /***********************************************************************/
! 605: /** **/
! 606: /** FACTORISATION **/
! 607: /** **/
! 608: /***********************************************************************/
! 609: #define LT 17
! 610: #define assign_or_fail(x,y) {\
! 611: if (y==NULL) y=x; else if (!gegal(x,y)) return 0;\
! 612: }
! 613: #define tsh 6
! 614: #define typs(x,y) ((x << tsh) | y)
! 615: #define typ1(x) (x >> tsh)
! 616: #define typ2(x) (x & ((1<<tsh)-1))
! 617:
! 618: static long
! 619: poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)
! 620: {
! 621: long t[LT]; /* code pour 0,1,2,3,61,62,63,67,7,81,82,83,86,87,91,93,97 */
! 622: long tx = typ(x),lx,i,j,s,pa=BIGINT;
! 623: GEN pcx=NULL, p=NULL,pol=NULL,p1,p2;
! 624:
! 625: if (is_scalar_t(tx))
! 626: {
! 627: if (tx == t_POLMOD) return 0;
! 628: x = scalarpol(x,0);
! 629: }
! 630: for (i=2; i<LT; i++) t[i]=0; /* t[0..1] unused */
! 631: lx = lgef(x);
! 632: for (i=2; i<lx; i++)
! 633: {
! 634: p1=(GEN)x[i];
! 635: switch(typ(p1))
! 636: {
! 637: case t_INT: case t_FRAC: case t_FRACN:
! 638: break;
! 639: case t_REAL:
! 640: s = precision(p1); if (s < pa) pa = s;
! 641: t[2]=1; break;
! 642: case t_INTMOD:
! 643: assign_or_fail((GEN)p1[1],p);
! 644: t[3]=1; break;
! 645: case t_COMPLEX:
! 646: if (!pcx)
! 647: {
! 648: pcx = cgetg(5,t_POL); /* x^2 + 1 */
! 649: pcx[1] = evalsigne(1)|evalvarn(0)|m_evallgef(5),
! 650: pcx[2]=pcx[4]=un; pcx[3]=zero;
! 651: }
! 652: for (j=1; j<=2; j++)
! 653: {
! 654: p2 = (GEN)p1[j];
! 655: switch(typ(p2))
! 656: {
! 657: case t_INT: case t_FRAC: case t_FRACN:
! 658: assign_or_fail(pcx,pol);
! 659: t[4]=1; break;
! 660: case t_REAL:
! 661: s = precision(p2); if (s < pa) pa = s;
! 662: t[5]=1; break;
! 663: case t_INTMOD:
! 664: assign_or_fail((GEN)p2[1],p);
! 665: assign_or_fail(pcx,pol);
! 666: t[6]=1; break;
! 667: case t_PADIC:
! 668: s = precp(p2) + valp(p2); if (s < pa) pa = s;
! 669: assign_or_fail((GEN)p2[2],p);
! 670: assign_or_fail(pcx,pol);
! 671: t[7]=1; break;
! 672: default: return 0;
! 673: }
! 674: }
! 675: break;
! 676: case t_PADIC:
! 677: s = precp(p1) + valp(p1); if (s < pa) pa = s;
! 678: assign_or_fail((GEN)p1[2],p);
! 679: t[8]=1; break;
! 680: case t_QUAD:
! 681: for (j=2; j<=3; j++)
! 682: {
! 683: p2 = (GEN)p1[j];
! 684: switch(typ(p2))
! 685: {
! 686: case t_INT: case t_FRAC: case t_FRACN:
! 687: assign_or_fail((GEN)p1[1],pol);
! 688: t[9]=1; break;
! 689: case t_REAL:
! 690: s = precision(p2); if (s < pa) pa = s;
! 691: if (gsigne(discsr((GEN)p1[1]))>0) t[10]=1; else t[12]=1;
! 692: break;
! 693: case t_INTMOD:
! 694: assign_or_fail((GEN)p2[1],p);
! 695: assign_or_fail((GEN)p1[1],pol);
! 696: t[11]=1; break;
! 697: case t_PADIC:
! 698: s = precp(p2) + valp(p2); if (s < pa) pa = s;
! 699: assign_or_fail((GEN)p2[2],p);
! 700: assign_or_fail((GEN)p1[1],pol);
! 701: t[13]=1; break;
! 702: default: return 0;
! 703: }
! 704: }
! 705: break;
! 706: case t_POLMOD:
! 707: assign_or_fail((GEN)p1[1],pol);
! 708: for (j=1; j<=2; j++)
! 709: {
! 710: GEN pbis = NULL, polbis = NULL;
! 711: long pabis;
! 712: switch(poltype((GEN)p1[j],&pbis,&polbis,&pabis))
! 713: {
! 714: case t_INT: t[14]=1; break;
! 715: case t_INTMOD: t[15]=1; break;
! 716: case t_PADIC: t[16]=1; if (pabis<pa) pa=pabis; break;
! 717: default: return 0;
! 718: }
! 719: if (pbis) assign_or_fail(pbis,p);
! 720: if (polbis) assign_or_fail(polbis,pol);
! 721: }
! 722: break;
! 723: default: return 0;
! 724: }
! 725: }
! 726: if (t[5]||t[12])
! 727: {
! 728: if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
! 729: *ptpa=pa; return t_COMPLEX;
! 730: }
! 731: if (t[2]||t[10])
! 732: {
! 733: if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
! 734: *ptpa=pa; return t[4]?t_COMPLEX:t_REAL;
! 735: }
! 736: if (t[6]||t[11]||t[15])
! 737: {
! 738: *ptpol=pol; *ptp=p;
! 739: i = t[15]? t_POLMOD: (t[11]? t_QUAD: t_COMPLEX);
! 740: return typs(i, t_INTMOD);
! 741: }
! 742: if (t[7]||t[13]||t[16])
! 743: {
! 744: *ptpol=pol; *ptp=p; *ptpa=pa;
! 745: i = t[16]? t_POLMOD: (t[13]? t_QUAD: t_COMPLEX);
! 746: return typs(i, t_PADIC);
! 747: }
! 748: if (t[4]||t[9]||t[14])
! 749: {
! 750: *ptpol=pol;
! 751: i = t[14]? t_POLMOD: (t[9]? t_QUAD: t_COMPLEX);
! 752: return typs(i, t_INT);
! 753: }
! 754: if (t[3]) { *ptp=p; return t_INTMOD; }
! 755: if (t[8]) { *ptp=p; *ptpa=pa; return t_PADIC; }
! 756: return t_INT;
! 757: }
! 758: #undef LT
! 759:
! 760: GEN
! 761: factor0(GEN x,long flag)
! 762: {
! 763: long tx=typ(x);
! 764:
! 765: if (flag<0) return factor(x);
! 766: if (is_matvec_t(tx)) return gboundfact(x,flag);
! 767: if (tx==t_INT || tx==t_FRAC || tx==t_FRACN) return boundfact(x,flag);
! 768: err(talker,"partial factorization is not meaningful here");
! 769: return NULL; /* not reached */
! 770: }
! 771:
! 772: static GEN
! 773: poldeg1(long v, GEN x0, GEN x1)
! 774: {
! 775: GEN x = cgetg(4,t_POL);
! 776: x[1] = evalsigne(1) | evalvarn(v) | evallgef(4);
! 777: x[2] = (long)x0; x[3] = (long)x1; return normalizepol(x);
! 778: }
! 779:
! 780: GEN
! 781: factor(GEN x)
! 782: {
! 783: long tx=typ(x),lx,av,tetpil,i,j,pa,v,r1;
! 784: GEN y,p,p1,p2,p3,p4,p5,pol;
! 785:
! 786: if (is_matvec_t(tx))
! 787: {
! 788: lx=lg(x); y=cgetg(lx,tx);
! 789: for (i=1; i<lx; i++) y[i]=(long)factor((GEN)x[i]);
! 790: return y;
! 791: }
! 792: if (gcmp0(x))
! 793: {
! 794: y=cgetg(3,t_MAT);
! 795: p1=cgetg(2,t_COL); y[1]=(long)p1; p1[1]=zero;
! 796: p2=cgetg(2,t_COL); y[2]=(long)p2; p2[1]=un;
! 797: return y;
! 798: }
! 799: switch(tx)
! 800: {
! 801: case t_INT: return decomp(x);
! 802:
! 803: case t_FRACN:
! 804: av=avma; x=gred(x); /* fall through */
! 805: case t_FRAC:
! 806: if (tx==t_FRAC) av=avma;
! 807: p1 = decomp((GEN)x[1]);
! 808: p2 = decomp((GEN)x[2]);
! 809: p4 = concatsp((GEN)p1[1], (GEN)p2[1]);
! 810: p5 = concatsp((GEN)p1[2], gneg_i((GEN)p2[2]));
! 811: p3 = sindexsort(p4); lx = lg(p3);
! 812: tetpil=avma; y=cgetg(3,t_MAT);
! 813: p1 = cgetg(lx,t_COL); y[1]=(long)p1;
! 814: p2 = cgetg(lx,t_COL); y[2]=(long)p2;
! 815: for (i=1; i<lx; i++)
! 816: {
! 817: p1[i] = licopy((GEN)p4[p3[i]]);
! 818: p2[i] = licopy((GEN)p5[p3[i]]);
! 819: }
! 820: return gerepile(av,tetpil,y);
! 821:
! 822: case t_POL:
! 823: tx=poltype(x,&p,&pol,&pa);
! 824: switch(tx)
! 825: {
! 826: case 0: err(impl,"factor for general polynomials");
! 827: case t_INT: return factpol(x,0,1);
! 828: case t_INTMOD: return factmod(x,p);
! 829:
! 830: case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
! 831: av=avma; p1=roots(x,pa); tetpil=avma;
! 832: p2=cgetg(lx,t_COL);
! 833: for (i=1; i<lx; i++)
! 834: p2[i] = (long)poldeg1(v, gneg((GEN)p1[i]), gun);
! 835: y[1]=lpile(av,tetpil,p2);
! 836: p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
! 837: y[2]=(long)p3; return y;
! 838:
! 839: case t_REAL: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
! 840: av=avma; p1=roots(x,pa); tetpil=avma;
! 841: for(r1=1; r1<lx; r1++)
! 842: if (signe(gmael(p1,r1,2))) break;
! 843: lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
! 844: for(i=1; i<r1; i++)
! 845: p2[i] = (long)poldeg1(v, negr(gmael(p1,i,1)), gun);
! 846: for( ; i<lx; i++)
! 847: {
! 848: GEN a = (GEN) p1[2*i-r1];
! 849: p = cgetg(5, t_POL); p2[i] = (long)p;
! 850: p[1] = evalsigne(1) | evalvarn(v) | evallgef(5);
! 851: p[2] = lnorm(a);
! 852: p[3] = lmul2n((GEN)a[1],1); setsigne(p[3],-signe(p[3]));
! 853: p[4] = un;
! 854: }
! 855: y[1]=lpile(av,tetpil,p2);
! 856: p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
! 857: y[2]=(long)p3; return y;
! 858:
! 859: case t_PADIC: return factorpadic4(x,p,pa);
! 860:
! 861: default:
! 862: {
! 863: long killv;
! 864: av=avma; x = dummycopy(x); lx=lgef(x);
! 865: pol = dummycopy(pol);
! 866: v = manage_var(4,NULL);
! 867: for(i=2; i<lx; i++)
! 868: {
! 869: p1=(GEN)x[i];
! 870: switch(typ(p1))
! 871: {
! 872: case t_QUAD: p1++;
! 873: case t_COMPLEX:
! 874: p2 = cgetg(3, t_POLMOD); x[i] = (long) p2;
! 875: p2[1] = (long)pol;
! 876: p2[2] = (long)poldeg1(v, (GEN)p1[1],(GEN)p1[2]);
! 877: }
! 878: }
! 879: killv = (avma != (long)pol);
! 880: if (killv) setvarn(pol, fetch_var());
! 881: switch (typ2(tx))
! 882: {
! 883: case t_INT: p1 = polfnf(x,pol); break;
! 884: case t_INTMOD: p1 = factmod9(x,p,pol); break;
! 885: default: err(impl,"factor of general polynomial");
! 886: }
! 887: switch (typ1(tx))
! 888: {
! 889: case t_POLMOD:
! 890: if (killv) delete_var();
! 891: return gerepileupto(av,p1);
! 892: case t_COMPLEX: p5 = gi; break;
! 893: case t_QUAD: p5=cgetg(4,t_QUAD);
! 894: p5[1]=(long)pol; p5[2]=zero; p5[3]=un;
! 895: break;
! 896: default: err(impl,"factor of general polynomial");
! 897: }
! 898: p2=(GEN)p1[1];
! 899: for(i=1; i<lg(p2); i++)
! 900: {
! 901: p3=(GEN)p2[i];
! 902: for(j=2; j<lgef(p3); j++)
! 903: {
! 904: p4=(GEN)p3[j];
! 905: if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5);
! 906: }
! 907: }
! 908: if (killv) delete_var();
! 909: tetpil=avma; y=cgetg(3,t_MAT);
! 910: y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]);
! 911: return gerepile(av,tetpil,y);
! 912: }
! 913: }
! 914:
! 915: case t_RFRACN:
! 916: av=avma; x=gred_rfrac(x); /* fall through */
! 917: case t_RFRAC:
! 918: if (tx==t_RFRAC) av=avma;
! 919: p1=factor((GEN)x[1]);
! 920: p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]);
! 921: tetpil=avma; y=cgetg(3,t_MAT);
! 922: y[1]=lconcat((GEN)p1[1],(GEN)p2[1]);
! 923: y[2]=lconcat((GEN)p1[2],p3);
! 924: return gerepile(av,tetpil,y);
! 925: }
! 926: err(impl,"general factorization");
! 927: return NULL; /* not reached */
! 928: }
! 929: #undef typ1
! 930: #undef typ2
! 931: #undef typs
! 932: #undef tsh
! 933:
! 934: GEN
! 935: divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
! 936: {
! 937: long i,k,lx = lg(x);
! 938:
! 939: if (lx == 1) return gun;
! 940: if (lx == 2) return gcopy((GEN)x[1]);
! 941: x = dummycopy(x); k = lx;
! 942: while (k > 2)
! 943: {
! 944: if (DEBUGLEVEL>7)
! 945: fprintferr("prod: remaining objects %ld\n",k-1);
! 946: lx = k; k = 1;
! 947: for (i=1; i<lx-1; i+=2)
! 948: x[k++] = (long)mul((GEN)x[i],(GEN)x[i+1]);
! 949: if (i < lx) x[k++] = x[i];
! 950: }
! 951: return (GEN)x[1];
! 952: }
! 953:
! 954: static GEN static_nf;
! 955:
! 956: static GEN
! 957: myidealmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }
! 958:
! 959: static GEN
! 960: myidealpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }
! 961:
! 962: GEN
! 963: factorback(GEN fa, GEN nf)
! 964: {
! 965: long av=avma,k,l,lx;
! 966: GEN ex, x;
! 967: GEN (*_mul)(GEN,GEN);
! 968: GEN (*_pow)(GEN,GEN);
! 969:
! 970: if (typ(fa)!=t_MAT || lg(fa)!=3)
! 971: err(talker,"incorrect factorisation in factorback");
! 972: ex=(GEN)fa[2]; fa=(GEN)fa[1];
! 973: lx = lg(fa); if (lx == 1) return gun;
! 974: x = cgetg(lx,t_VEC);
! 975: if (nf)
! 976: {
! 977: static_nf = nf;
! 978: _mul = &myidealmul;
! 979: _pow = &myidealpow;
! 980: }
! 981: else
! 982: {
! 983: _mul = &gmul;
! 984: _pow = &powgi;
! 985: }
! 986: for (l=1,k=1; k<lx; k++)
! 987: if (signe(ex[k]))
! 988: x[l++] = (long)_pow((GEN)fa[k],(GEN)ex[k]);
! 989: setlg(x,l);
! 990: return gerepileupto(av, divide_conquer_prod(x, _mul));
! 991: }
! 992:
! 993: GEN
! 994: gisirreducible(GEN x)
! 995: {
! 996: long av=avma, tx = typ(x),l,i;
! 997: GEN y;
! 998:
! 999: if (is_matvec_t(tx))
! 1000: {
! 1001: l=lg(x); y=cgetg(l,tx);
! 1002: for (i=1; i<l; i++) y[i]=(long)gisirreducible((GEN)x[i]);
! 1003: return y;
! 1004: }
! 1005: if (is_intreal_t(tx) || is_frac_t(tx)) return gzero;
! 1006: if (tx!=t_POL) err(notpoler,"gisirreducible");
! 1007: l=lgef(x); if (l<=3) return gzero;
! 1008: y=factor(x); avma=av;
! 1009: return (lgef(gcoeff(y,1,1))==l)?gun:gzero;
! 1010: }
! 1011:
! 1012: /*******************************************************************/
! 1013: /* */
! 1014: /* PGCD GENERAL */
! 1015: /* */
! 1016: /*******************************************************************/
! 1017:
! 1018: GEN
! 1019: gcd0(GEN x, GEN y, long flag)
! 1020: {
! 1021: switch(flag)
! 1022: {
! 1023: case 0: return ggcd(x,y);
! 1024: case 1: return modulargcd(x,y);
! 1025: case 2: return srgcd(x,y);
! 1026: default: err(flagerr,"gcd");
! 1027: }
! 1028: return NULL; /* not reached */
! 1029: }
! 1030:
! 1031: /* x is a COMPLEX or a QUAD */
! 1032: static GEN
! 1033: triv_cont_gcd(GEN x, GEN y)
! 1034: {
! 1035: long av = avma, tetpil;
! 1036: GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2])
! 1037: : ggcd((GEN)x[2],(GEN)x[3]);
! 1038: tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
! 1039: }
! 1040:
! 1041: static GEN
! 1042: cont_gcd(GEN x, GEN y)
! 1043: {
! 1044: long av = avma, tetpil;
! 1045: GEN p1 = content(x);
! 1046:
! 1047: tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
! 1048: }
! 1049:
! 1050: /* y is a PADIC, x a rational number or an INTMOD */
! 1051: static GEN
! 1052: padic_gcd(GEN x, GEN y)
! 1053: {
! 1054: long v = ggval(x,(GEN)y[2]), w = valp(y);
! 1055: if (w < v) v = w;
! 1056: return gpuigs((GEN)y[2], v);
! 1057: }
! 1058:
! 1059: #define fix_frac(z) if (signe(z[2])<0)\
! 1060: {\
! 1061: setsigne(z[1],-signe(z[1]));\
! 1062: setsigne(z[2],1);\
! 1063: }
! 1064:
! 1065: GEN
! 1066: ggcd(GEN x, GEN y)
! 1067: {
! 1068: long l,av,tetpil,i,vx,vy, tx = typ(x), ty = typ(y);
! 1069: GEN p1,z;
! 1070:
! 1071: if (tx>ty) { p1=x; x=y; y=p1; l=tx; tx=ty; ty=l; }
! 1072: if (is_matvec_t(ty))
! 1073: {
! 1074: l=lg(y); z=cgetg(l,ty);
! 1075: for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]);
! 1076: return z;
! 1077: }
! 1078: if (gcmp0(x)) return gcopy(y);
! 1079: if (gcmp0(y)) return gcopy(x);
! 1080: if (is_const_t(tx))
! 1081: {
! 1082: if (ty == t_FRACN)
! 1083: {
! 1084: if (tx==t_INTMOD)
! 1085: {
! 1086: av=avma; y = gred(y); tetpil=avma;
! 1087: return gerepile(av,tetpil,ggcd(x,y));
! 1088: }
! 1089: ty = t_FRAC;
! 1090: }
! 1091: if (tx == t_FRACN) tx = t_FRAC;
! 1092: if (ty == tx) switch(tx)
! 1093: {
! 1094: case t_INT:
! 1095: return mppgcd(x,y);
! 1096:
! 1097: case t_INTMOD: z=cgetg(3,t_INTMOD);
! 1098: if (egalii((GEN)x[1],(GEN)y[1]))
! 1099: { copyifstack(x[1],z[1]); }
! 1100: else
! 1101: z[1]=lmppgcd((GEN)x[1],(GEN)y[1]);
! 1102: if (gcmp1((GEN)z[1])) z[2]=zero;
! 1103: else
! 1104: {
! 1105: av=avma; p1=mppgcd((GEN)z[1],(GEN)x[2]);
! 1106: if (!gcmp1(p1))
! 1107: {
! 1108: tetpil=avma;
! 1109: p1=gerepile(av,tetpil,mppgcd(p1,(GEN)y[2]));
! 1110: }
! 1111: z[2]=(long)p1;
! 1112: }
! 1113: return z;
! 1114:
! 1115: case t_FRAC: z=cgetg(3,t_FRAC);
! 1116: (void)new_chunk(lgefint(x[2])+lgefint(y[2]));
! 1117: p1 = divii((GEN)y[2], mppgcd((GEN)x[2], (GEN)y[2]));
! 1118: avma = (long)z;
! 1119: z[2] = lmulii(p1, (GEN)x[2]);
! 1120: z[1] = (long)mppgcd((GEN)x[1], (GEN)y[1]);
! 1121: fix_frac(z); return z;
! 1122:
! 1123: case t_COMPLEX:
! 1124: av=avma; p1=gdiv(x,y);
! 1125: if (gcmp0((GEN)p1[2]))
! 1126: {
! 1127: p1=(GEN)p1[1];
! 1128: switch(typ(p1))
! 1129: {
! 1130: case t_INT:
! 1131: avma=av; return gcopy(y);
! 1132: case t_FRAC: case t_FRACN:
! 1133: tetpil=avma; return gerepile(av,tetpil,gdiv(y,(GEN)p1[2]));
! 1134: default: avma=av; return gun;
! 1135: }
! 1136: }
! 1137: avma=av;
! 1138: if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(y);
! 1139: p1=gdiv(y,x); avma=av;
! 1140: if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(x);
! 1141: return triv_cont_gcd(y,x);
! 1142:
! 1143: case t_PADIC:
! 1144: if (!egalii((GEN)x[2],(GEN)y[2])) return gun;
! 1145: return gpuigs((GEN)y[2],min(valp(y),valp(x)));
! 1146:
! 1147: case t_QUAD:
! 1148: av=avma; p1=gdiv(x,y);
! 1149: if (gcmp0((GEN)p1[3]))
! 1150: {
! 1151: p1=(GEN)p1[2];
! 1152: if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
! 1153: tetpil=avma; return gerepile(av,tetpil, gdiv(y,(GEN)p1[2]));
! 1154: }
! 1155: avma=av;
! 1156: if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(y);
! 1157: p1=gdiv(y,x); avma=av;
! 1158: if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(x);
! 1159: return triv_cont_gcd(y,x);
! 1160:
! 1161: default: return gun; /* t_REAL */
! 1162: }
! 1163: if (is_const_t(ty)) switch(tx)
! 1164: {
! 1165: case t_INT:
! 1166: switch(ty)
! 1167: {
! 1168: case t_INTMOD: z=cgetg(3,t_INTMOD);
! 1169: copyifstack(y[1],z[1]); av=avma;
! 1170: p1=mppgcd((GEN)y[1],(GEN)y[2]);
! 1171: if (!gcmp1(p1))
! 1172: { tetpil=avma; p1=gerepile(av,tetpil,mppgcd(x,p1)); }
! 1173: z[2]=(long)p1; return z;
! 1174:
! 1175: case t_FRAC: z=cgetg(3,t_FRAC);
! 1176: z[1]=lmppgcd(x,(GEN)y[1]);
! 1177: z[2]=licopy((GEN)y[2]); return z;
! 1178:
! 1179: case t_COMPLEX: case t_QUAD:
! 1180: return triv_cont_gcd(y,x);
! 1181:
! 1182: case t_PADIC:
! 1183: return padic_gcd(x,y);
! 1184:
! 1185: default: /* t_REAL */
! 1186: return gun;
! 1187: }
! 1188:
! 1189: case t_INTMOD:
! 1190: switch(ty)
! 1191: {
! 1192: case t_FRAC:
! 1193: av=avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); tetpil=avma;
! 1194: if (!gcmp1(p1))
! 1195: err(talker,"non invertible fraction in a gcd with INTMOD");
! 1196: return gerepile(av,tetpil, ggcd((GEN)y[1],x));
! 1197:
! 1198: case t_COMPLEX: case t_QUAD:
! 1199: return triv_cont_gcd(y,x);
! 1200:
! 1201: case t_PADIC:
! 1202: return padic_gcd(x,y);
! 1203: }
! 1204:
! 1205: case t_FRAC:
! 1206: switch(ty)
! 1207: {
! 1208: case t_COMPLEX: case t_QUAD:
! 1209: return triv_cont_gcd(y,x);
! 1210:
! 1211: case t_PADIC:
! 1212: return padic_gcd(x,y);
! 1213: }
! 1214:
! 1215: case t_COMPLEX: /* ty = PADIC or QUAD */
! 1216: return triv_cont_gcd(x,y);
! 1217:
! 1218: case t_PADIC: /* ty = QUAD */
! 1219: return triv_cont_gcd(y,x);
! 1220:
! 1221: default: /* tx = t_REAL */
! 1222: return gun;
! 1223: }
! 1224: return cont_gcd(y,x);
! 1225: }
! 1226:
! 1227: vx=gvar9(x); vy=gvar9(y);
! 1228: if (vy<vx) return cont_gcd(y,x);
! 1229: if (vx<vy) return cont_gcd(x,y);
! 1230: switch(tx)
! 1231: {
! 1232: case t_POLMOD:
! 1233: switch(ty)
! 1234: {
! 1235: case t_POLMOD: z=cgetg(3,t_POLMOD);
! 1236: if (gegal((GEN)x[1],(GEN)y[1]))
! 1237: { copyifstack(x[1],z[1]); }
! 1238: else
! 1239: z[1] = lgcd((GEN)x[1],(GEN)y[1]);
! 1240: if (lgef(z[1])<=3) z[2]=zero;
! 1241: else
! 1242: {
! 1243: av=avma; p1=ggcd((GEN)z[1],(GEN)x[2]);
! 1244: if (lgef(p1)>3)
! 1245: {
! 1246: tetpil=avma;
! 1247: p1=gerepile(av,tetpil,ggcd(p1,(GEN)y[2]));
! 1248: }
! 1249: z[2]=(long)p1;
! 1250: }
! 1251: return z;
! 1252:
! 1253: case t_POL: z=cgetg(3,t_POLMOD);
! 1254: copyifstack(x[1],z[1]); av=avma;
! 1255: p1=ggcd((GEN)x[1],(GEN)x[2]);
! 1256: if (lgef(p1)>3)
! 1257: { tetpil=avma; p1=gerepile(av,tetpil,ggcd(y,p1)); }
! 1258: z[2]=(long)p1; return z;
! 1259:
! 1260: case t_RFRAC:
! 1261: av=avma; p1=ggcd((GEN)x[1],(GEN)y[2]); tetpil=avma;
! 1262: if (!gcmp1(p1))
! 1263: err(talker,"non invertible fraction in a gcd with POLMOD");
! 1264: return gerepile(av,tetpil,ggcd((GEN)y[1],x));
! 1265:
! 1266: case t_RFRACN:
! 1267: av=avma; p1=gred_rfrac(y); tetpil=avma;
! 1268: return gerepile(av,tetpil,ggcd(p1,x));
! 1269: }
! 1270: break;
! 1271:
! 1272: case t_POL:
! 1273: switch(ty)
! 1274: {
! 1275: case t_POL:
! 1276: return srgcd(x,y);
! 1277:
! 1278: case t_SER:
! 1279: return gpuigs(polx[vx],min(valp(y),gval(x,vx)));
! 1280:
! 1281: case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty);
! 1282: z[1]=lgcd(x,(GEN)y[1]);
! 1283: z[2]=lcopy((GEN)y[2]); return z;
! 1284: }
! 1285: break;
! 1286:
! 1287: case t_SER:
! 1288: switch(ty)
! 1289: {
! 1290: case t_SER:
! 1291: return gpuigs(polx[vx],min(valp(x),valp(y)));
! 1292:
! 1293: case t_RFRAC: case t_RFRACN:
! 1294: return gpuigs(polx[vx],min(valp(x),gval(y,vx)));
! 1295: }
! 1296: break;
! 1297:
! 1298: case t_RFRAC: case t_RFRACN: z=cgetg(3,ty);
! 1299: if (!is_rfrac_t(ty))
! 1300: err(talker,"forbidden gcd rational function with vector/matrix");
! 1301: p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2]));
! 1302: tetpil = avma;
! 1303: z[2] = lpile((long)z,tetpil,gmul(p1, (GEN)x[2]));
! 1304: z[1] = lgcd((GEN)x[1], (GEN)y[1]); return z;
! 1305: }
! 1306: err(talker,"gcd vector/matrix with a forbidden type");
! 1307: return NULL; /* not reached */
! 1308: }
! 1309:
! 1310: GEN
! 1311: glcm(GEN x, GEN y)
! 1312: {
! 1313: long av=avma,tx,ty=typ(y),i,l;
! 1314: GEN p1,p2,z;
! 1315:
! 1316: if (is_matvec_t(ty))
! 1317: {
! 1318: l=lg(y); z=cgetg(l,ty);
! 1319: for (i=1; i<l; i++) z[i]=(long)glcm(x,(GEN)y[i]);
! 1320: return z;
! 1321: }
! 1322: tx=typ(x);
! 1323: if (is_matvec_t(tx))
! 1324: {
! 1325: l=lg(x); z=cgetg(l,tx);
! 1326: for (i=1; i<l; i++) z[i]=(long)glcm((GEN)x[i],y);
! 1327: return z;
! 1328: }
! 1329: if (gcmp0(x)) return gzero;
! 1330: if (tx==t_INT && ty==t_INT)
! 1331: {
! 1332: p1 = mppgcd(x,y); if (is_pm1(p1)) { avma = av; return mulii(x,y); }
! 1333: p2 = mulii(divii(y,p1), x);
! 1334: if (signe(p2)<0) setsigne(p2,1);
! 1335: return gerepileupto(av, p2);
! 1336: }
! 1337: p1 = ggcd(x,y); if (gcmp1(p1)) {avma = av; return gmul(x,y); }
! 1338: p2 = gmul(gdiv(y,p1), x);
! 1339: if (typ(p2)==t_INT && signe(p2)<0) setsigne(p2,1);
! 1340: if (typ(p2)==t_POL)
! 1341: {
! 1342: p1=leading_term(p2);
! 1343: if (typ(p1)==t_INT && signe(p1)<0) p2=gneg(p2);
! 1344: }
! 1345: return gerepileupto(av,p2);
! 1346: }
! 1347:
! 1348: static GEN
! 1349: polgcdnun(GEN x, GEN y)
! 1350: {
! 1351: long av1, av = avma, lim = stack_lim(av,1);
! 1352: GEN p1, yorig = y;
! 1353:
! 1354: for(;;)
! 1355: {
! 1356: av1 = avma; p1 = gres(x,y);
! 1357: if (gcmp0(p1))
! 1358: {
! 1359: avma = av1;
! 1360: return (y==yorig)? gcopy(y): gerepileupto(av,y);
! 1361: }
! 1362: x = y; y = p1;
! 1363: if (low_stack(lim,stack_lim(av,1)))
! 1364: {
! 1365: GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y;
! 1366: if(DEBUGMEM>1) err(warnmem,"polgcdnun");
! 1367: gerepilemanysp(av,av1,gptr,2);
! 1368: }
! 1369: }
! 1370: }
! 1371:
! 1372: /* return 1 if coeff explosion is not possible */
! 1373: static int
! 1374: issimplefield(GEN x)
! 1375: {
! 1376: long lx,i;
! 1377: switch(typ(x))
! 1378: {
! 1379: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
! 1380: return 1;
! 1381: case t_POL:
! 1382: lx=lgef(x);
! 1383: for (i=2; i<lx; i++)
! 1384: if (!issimplefield((GEN)x[i])) return 0;
! 1385: return 1;
! 1386: case t_COMPLEX: case t_POLMOD:
! 1387: return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]);
! 1388: }
! 1389: return 0;
! 1390: }
! 1391:
! 1392: static int
! 1393: isrational(GEN x)
! 1394: {
! 1395: long i;
! 1396: for (i=lgef(x)-1; i>1; i--)
! 1397: switch(typ(x[i]))
! 1398: {
! 1399: case t_INT:
! 1400: case t_FRAC: break;
! 1401: default: return 0;
! 1402: }
! 1403: return 1;
! 1404: }
! 1405:
! 1406: static int
! 1407: isinexactfield(GEN x)
! 1408: {
! 1409: long lx,i;
! 1410: switch(typ(x))
! 1411: {
! 1412: case t_REAL: case t_PADIC: case t_SER:
! 1413: return 1;
! 1414: case t_POL:
! 1415: lx=lgef(x);
! 1416: for (i=2; i<lx; i++)
! 1417: if (!isinexactfield((GEN)x[i])) return 0;
! 1418: return 1;
! 1419: case t_COMPLEX: case t_POLMOD:
! 1420: return isinexactfield((GEN)x[1]) || isinexactfield((GEN)x[2]);
! 1421: }
! 1422: return 0;
! 1423: }
! 1424:
! 1425: static GEN
! 1426: gcdmonome(GEN x, GEN y)
! 1427: {
! 1428: long tetpil,av=avma, lx=lgef(x), v=varn(x), e=gval(y,v);
! 1429: GEN p1,p2;
! 1430:
! 1431: if (lx-3<e) e=lx-3;
! 1432: p1=ggcd(leading_term(x),content(y)); p2=gpuigs(polx[v],e);
! 1433: tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));
! 1434: }
! 1435:
! 1436: /***********************************************************************/
! 1437: /** **/
! 1438: /** BEZOUT GENERAL **/
! 1439: /** **/
! 1440: /***********************************************************************/
! 1441:
! 1442: static GEN
! 1443: polinvinexact(GEN x, GEN y)
! 1444: {
! 1445: long i,dx=lgef(x)-3,dy=lgef(y)-3,lz=dx+dy, av=avma, tetpil;
! 1446: GEN v,z;
! 1447:
! 1448: z=cgetg(dy+2,t_POL); z[1]=y[1];
! 1449: v=cgetg(lz+1,t_COL);
! 1450: for (i=1; i<lz; i++) v[i]=zero;
! 1451: v[lz]=un; v=gauss(sylvestermatrix(y,x),v);
! 1452: for (i=2; i<dy+2; i++) z[i]=v[lz-i+2];
! 1453: z = normalizepol_i(z, dy+2); tetpil = avma;
! 1454: return gerepile(av,tetpil,gcopy(z));
! 1455: }
! 1456:
! 1457: static GEN
! 1458: polinvmod(GEN x, GEN y)
! 1459: {
! 1460: long av,av1,tx,vx=varn(x),vy=varn(y);
! 1461: GEN u,v,d,p1;
! 1462:
! 1463: while (vx!=vy)
! 1464: {
! 1465: if (vx > vy)
! 1466: {
! 1467: d=cgetg(3,t_RFRAC);
! 1468: d[1]=(long)polun[vx];
! 1469: d[2]=lcopy(x); return d;
! 1470: }
! 1471: if (lgef(x)!=3) err(talker,"non-invertible polynomial in polinvmod");
! 1472: x=(GEN)x[2]; vx=gvar(x);
! 1473: }
! 1474: tx=typ(x);
! 1475: if (tx==t_POL)
! 1476: {
! 1477: if (isinexactfield(x) || isinexactfield(y))
! 1478: return polinvinexact(x,y);
! 1479:
! 1480: av=avma; d=subresext(x,y,&u,&v);
! 1481: if (gcmp0(d)) err(talker,"non-invertible polynomial in polinvmod");
! 1482: if (typ(d)==t_POL && varn(d)==vx)
! 1483: {
! 1484: if (lgef(d)>3) err(talker,"non-invertible polynomial in polinvmod");
! 1485: d=(GEN)d[2];
! 1486: }
! 1487: av1=avma; return gerepile(av,av1,gdiv(u,d));
! 1488: }
! 1489: if (!is_rfrac_t(tx)) err(typeer,"polinvmod");
! 1490: av=avma; p1=gmul((GEN)x[1],polinvmod((GEN)x[2],y));
! 1491: av1=avma; return gerepile(av,av1,gmod(p1,y));
! 1492: }
! 1493:
! 1494: GEN
! 1495: gbezout(GEN x, GEN y, GEN *u, GEN *v)
! 1496: {
! 1497: long tx=typ(x),ty=typ(y);
! 1498:
! 1499: if (tx==t_INT && ty==t_INT) return bezout(x,y,u,v);
! 1500: if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) err(typeer,"gbezout");
! 1501: return bezoutpol(x,y,u,v);
! 1502: }
! 1503:
! 1504: GEN
! 1505: vecbezout(GEN x, GEN y)
! 1506: {
! 1507: GEN z=cgetg(4,t_VEC);
! 1508: z[3]=(long)gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
! 1509: return z;
! 1510: }
! 1511:
! 1512: GEN
! 1513: vecbezoutres(GEN x, GEN y)
! 1514: {
! 1515: GEN z=cgetg(4,t_VEC);
! 1516: z[3]=(long)subresext(x,y,(GEN*)(z+1),(GEN*)(z+2));
! 1517: return z;
! 1518: }
! 1519:
! 1520: /*******************************************************************/
! 1521: /* */
! 1522: /* CONTENU ET PARTIE PRIMITIVE */
! 1523: /* */
! 1524: /*******************************************************************/
! 1525:
! 1526: GEN
! 1527: content(GEN x)
! 1528: {
! 1529: long av,tetpil,lx,i,tx=typ(x);
! 1530: GEN p1,p2;
! 1531:
! 1532: if (is_scalar_t(tx))
! 1533: return tx==t_POLMOD? content((GEN)x[2]): gcopy(x);
! 1534: av = avma;
! 1535: switch(tx)
! 1536: {
! 1537: case t_RFRAC: case t_RFRACN:
! 1538: p1=content((GEN)x[1]);
! 1539: p2=content((GEN)x[2]); tetpil=avma;
! 1540: return gerepile(av,tetpil,gdiv(p1,p2));
! 1541:
! 1542: case t_VEC: case t_COL: case t_MAT:
! 1543: lx = lg(x); if (lx==1) return gun;
! 1544: p1=content((GEN)x[1]);
! 1545: for (i=2; i<lx; i++) p1 = ggcd(p1,content((GEN)x[i]));
! 1546: return gerepileupto(av,p1);
! 1547:
! 1548: case t_POL:
! 1549: if (!signe(x)) return gzero;
! 1550: lx = lgef(x); break;
! 1551: case t_SER:
! 1552: if (!signe(x)) return gzero;
! 1553: lx = lg(x); break;
! 1554: case t_QFR: case t_QFI:
! 1555: lx = 4; break;
! 1556:
! 1557: default: err(typeer,"content");
! 1558: }
! 1559: for (i=lontyp[tx]; i<lx; i++)
! 1560: if (typ(x[i]) != t_INT) break;
! 1561: lx--; p1=(GEN)x[lx];
! 1562: if (i > lx)
! 1563: { /* integer coeffs */
! 1564: while (lx>lontyp[tx])
! 1565: {
! 1566: lx--; p1=mppgcd(p1,(GEN)x[lx]);
! 1567: if (is_pm1(p1)) { avma=av; return gun; }
! 1568: }
! 1569: }
! 1570: else
! 1571: {
! 1572: while (lx>lontyp[tx])
! 1573: {
! 1574: lx--; p1=ggcd(p1,(GEN)x[lx]);
! 1575: }
! 1576: if (isinexactreal(p1)) { avma=av; return gun; }
! 1577: }
! 1578: return av==avma? gcopy(p1): gerepileupto(av,p1);
! 1579: }
! 1580:
! 1581: /*******************************************************************/
! 1582: /* */
! 1583: /* SOUS RESULTANT */
! 1584: /* */
! 1585: /*******************************************************************/
! 1586: GEN
! 1587: gdivexact(GEN x, GEN y)
! 1588: {
! 1589: long i,lx;
! 1590: GEN z;
! 1591: if (gcmp1(y)) return x;
! 1592: switch(typ(x))
! 1593: {
! 1594: case t_INT:
! 1595: if (typ(y)==t_INT) return divii(x,y);
! 1596: if (!signe(x)) return gzero;
! 1597: break;
! 1598: case t_INTMOD:
! 1599: case t_POLMOD: return gmul(x,ginv(y));
! 1600: case t_POL:
! 1601: switch(typ(y))
! 1602: {
! 1603: case t_INTMOD:
! 1604: case t_POLMOD: return gmul(x,ginv(y));
! 1605: case t_POL:
! 1606: if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES_EXACT);
! 1607: }
! 1608: lx = lgef(x); z = cgetg(lx,t_POL);
! 1609: for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);
! 1610: z[1]=x[1]; return z;
! 1611: }
! 1612: if (DEBUGLEVEL) err(warner,"missing case in gdivexact");
! 1613: return gdiv(x,y);
! 1614: }
! 1615:
! 1616: static GEN
! 1617: init_resultant(GEN x, GEN y)
! 1618: {
! 1619: long tx,ty;
! 1620: if (gcmp0(x) || gcmp0(y)) return gzero;
! 1621: tx = typ(x); ty = typ(y);
! 1622: if (is_scalar_t(tx) || is_scalar_t(ty))
! 1623: {
! 1624: if (tx==t_POL) return gpuigs(y,lgef(x)-3);
! 1625: if (ty==t_POL) return gpuigs(x,lgef(y)-3);
! 1626: return gun;
! 1627: }
! 1628: if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall");
! 1629: if (varn(x)==varn(y)) return NULL;
! 1630: return (varn(x)<varn(y))? gpuigs(y,lgef(x)-3): gpuigs(x,lgef(y)-3);
! 1631: }
! 1632:
! 1633: /* return coefficients s.t x = x_0 X^n + ... + x_n */
! 1634: static GEN
! 1635: revpol(GEN x)
! 1636: {
! 1637: long i,n = lgef(x)-3;
! 1638: GEN y = cgetg(n+3,t_POL);
! 1639: y[1] = x[1]; x += 2; y += 2;
! 1640: for (i=0; i<=n; i++) y[i] = x[n-i];
! 1641: return y;
! 1642: }
! 1643:
! 1644: /* assume dx >= dy, y non constant */
! 1645: GEN
! 1646: pseudorem(GEN x, GEN y)
! 1647: {
! 1648: long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,lx, p;
! 1649:
! 1650: if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
! 1651: (void)new_chunk(2);
! 1652: dx=lgef(x)-3; x = revpol(x);
! 1653: dy=lgef(y)-3; y = revpol(y); dz=dx-dy; p = dz+1;
! 1654: av2 = avma; lim = stack_lim(av2,1);
! 1655: for (;;)
! 1656: {
! 1657: x[0] = lneg((GEN)x[0]); p--;
! 1658: for (i=1; i<=dy; i++)
! 1659: x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
! 1660: for ( ; i<=dx; i++)
! 1661: x[i] = lmul((GEN)y[0], (GEN)x[i]);
! 1662: do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0]));
! 1663: if (dx < dy) break;
! 1664: if (low_stack(lim,stack_lim(av2,1)))
! 1665: {
! 1666: if(DEBUGMEM>1) err(warnmem,"pseudorem dx = %ld >= %ld",dx,dy);
! 1667: gerepilemanycoeffs(av2,x,dx+1);
! 1668: }
! 1669: }
! 1670: if (dx < 0) return zeropol(vx);
! 1671: lx = dx+3; x -= 2;
! 1672: x[0]=evaltyp(t_POL) | evallg(lx);
! 1673: x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);
! 1674: x = revpol(x) - 2;
! 1675: return gerepileupto(av, gmul(x, gpowgs((GEN)y[0], p)));
! 1676: }
! 1677:
! 1678: /* Si sol != NULL:
! 1679: * met dans *sol le dernier polynome non nul de la polynomial remainder
! 1680: * sequence si elle a ete effectuee, 0 sinon
! 1681: */
! 1682: GEN
! 1683: subresall(GEN u, GEN v, GEN *sol)
! 1684: {
! 1685: long degq,av,av2,lim,tetpil,dx,dy,du,dv,dr,signh;
! 1686: GEN g,h,r,p1,p2,p3,p4;
! 1687:
! 1688: if (sol) *sol=gzero;
! 1689: if ((r = init_resultant(u,v))) return r;
! 1690:
! 1691: if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v);
! 1692: dx=lgef(u); dy=lgef(v); signh=1;
! 1693: if (dx<dy)
! 1694: {
! 1695: p1=u; u=v; v=p1; du=dx; dx=dy; dy=du;
! 1696: if ((dx&1) == 0 && (dy&1) == 0) signh = -signh;
! 1697: }
! 1698: if (dy==3) return gpowgs((GEN)v[2],dx-3);
! 1699: av=avma;
! 1700: p3=content(u); if (gcmp1(p3)) p3=NULL; else u=gdiv(u,p3);
! 1701: p4=content(v); if (gcmp1(p4)) p4=NULL; else v=gdiv(v,p4);
! 1702: g=gun; h=gun; av2 = avma; lim = stack_lim(av2,1);
! 1703: for(;;)
! 1704: {
! 1705: r = pseudorem(u,v); dr = lgef(r);
! 1706: if (dr==2)
! 1707: {
! 1708: if (sol) {avma = (long)(r+2); *sol=gerepileupto(av,v);} else avma = av;
! 1709: return gzero;
! 1710: }
! 1711: du=lgef(u); dv=lgef(v);
! 1712: degq=du-dv; u=v;
! 1713: p1 = g; g = leading_term(u);
! 1714: switch(degq)
! 1715: {
! 1716: case 0: break;
! 1717: case 1:
! 1718: p1 = gmul(h,p1); h = g; break;
! 1719: default:
! 1720: p1 = gmul(gpuigs(h,degq),p1);
! 1721: h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
! 1722: }
! 1723: if ((du&1) == 0 && (dv&1) == 0) signh = -signh;
! 1724: v = gdivexact(r,p1);
! 1725: if (dr==3) break;
! 1726: if (low_stack(lim,stack_lim(av2,1)))
! 1727: {
! 1728: GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
! 1729: if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);
! 1730: gerepilemany(av2,gptr,4);
! 1731: }
! 1732: }
! 1733: if (dv==4) { tetpil=avma; p2=gcopy((GEN)v[2]); }
! 1734: else
! 1735: {
! 1736: if (dv == 3) err(bugparier,"subres");
! 1737: p1=gpuigs((GEN)v[2],dv-3); p2=gpuigs(h,dv-4);
! 1738: tetpil=avma; p2=gdiv(p1,p2);
! 1739: }
! 1740: if (p3) { p1=gpuigs(p3,dy-3); tetpil=avma; p2=gmul(p2,p1); }
! 1741: if (p4) { p1=gpuigs(p4,dx-3); tetpil=avma; p2=gmul(p2,p1); }
! 1742: if (signh<0) { tetpil=avma; p2=gneg(p2); }
! 1743: {
! 1744: GEN *gptr[2]; gptr[0]=&p2; if (sol) { *sol=gcopy(u); gptr[1]=sol; }
! 1745: gerepilemanysp(av,tetpil,gptr,sol?2:1);
! 1746: return p2;
! 1747: }
! 1748: }
! 1749:
! 1750: static GEN
! 1751: scalar_res(GEN x, GEN y, GEN *U, GEN *V)
! 1752: {
! 1753: long dx=lgef(x)-4;
! 1754: *V=gpuigs(y,dx); *U=gzero; return gmul(y,*V);
! 1755: }
! 1756:
! 1757: /* calcule U et V tel que Ux+By=resultant(x,y) */
! 1758: GEN
! 1759: subresext(GEN x, GEN y, GEN *U, GEN *V)
! 1760: {
! 1761: long degq,av,tetpil,tx,ty,dx,dy,du,dv,dr,signh;
! 1762: GEN z,g,h,r,p1,p2,p3,p4,u,v,lpu,um1,uze, *gptr[2];
! 1763:
! 1764: if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; }
! 1765: tx=typ(x); ty=typ(y);
! 1766: if (is_scalar_t(tx) || is_scalar_t(ty))
! 1767: {
! 1768: if (tx==t_POL) return scalar_res(x,y,U,V);
! 1769: if (ty==t_POL) return scalar_res(y,x,V,U);
! 1770: *U=ginv(x); *V=gzero; return gun;
! 1771: }
! 1772: if (tx!=t_POL || ty!=t_POL) err(typeer,"subresext");
! 1773: if (varn(x)!=varn(y))
! 1774: return (varn(x)<varn(y))? scalar_res(x,y,U,V)
! 1775: : scalar_res(y,x,V,U);
! 1776: dx=lgef(x); dy=lgef(y); signh=1;
! 1777: if (dx<dy)
! 1778: {
! 1779: GEN *W = U; U=V; V=W;
! 1780: du=dx; dx=dy; dy=du; p1=x; x=y; y=p1;
! 1781: if ((dx&1) == 0 && (dy&1) == 0) signh = -signh;
! 1782: }
! 1783: if (dy==3)
! 1784: {
! 1785: *V=gpuigs((GEN)y[2],dx-4);
! 1786: *U=gzero; return gmul(*V,(GEN)y[2]);
! 1787: }
! 1788: av=avma;
! 1789: p3=content(x); if (gcmp1(p3)) {p3 = NULL; u=x; } else u=gdiv(x,p3);
! 1790: p4=content(y); if (gcmp1(p4)) {p4 = NULL; v=y; } else v=gdiv(y,p4);
! 1791: g=gun; h=gun; um1=gun; uze=gzero;
! 1792: for(;;)
! 1793: {
! 1794: du=lgef(u); dv=lgef(v); degq=du-dv;
! 1795: lpu=gpuigs((GEN)v[dv-1],degq+1);
! 1796: p1=gmul(lpu,u); p2=poldivres(p1,v,&r);
! 1797: dr=lgef(r); if (dr==2) { *U=gzero; *V=gzero; avma=av; return gzero; }
! 1798:
! 1799: p1=gsub(gmul(lpu,um1),gmul(p2,uze));
! 1800: um1=uze; uze=p1; u=v;
! 1801: p1 = g; g = leading_term(u);
! 1802: switch(degq)
! 1803: {
! 1804: case 0: break;
! 1805: case 1: p1 = gmul(h,p1); h = g; break;
! 1806: default:
! 1807: p1 = gmul(gpuigs(h,degq),p1);
! 1808: h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
! 1809: }
! 1810: if ((du & 1) == 0 && (dv & 1) == 0) signh= -signh;
! 1811: v=gdivexact(r,p1); uze=gdivexact(uze,p1);
! 1812: if (dr==3) break;
! 1813: }
! 1814:
! 1815: p2=(dv==4)?gun:gpuigs(gdiv((GEN)v[2],h),dv-4);
! 1816: if (p3) p2 = gmul(p2,gpuigs(p3,dy-3));
! 1817: if (p4) p2 = gmul(p2,gpuigs(p4,dx-3));
! 1818: if (signh<0) p2=gneg_i(p2);
! 1819: p3 = p3? gdiv(p2,p3): p2;
! 1820:
! 1821: tetpil=avma; z=gmul((GEN)v[2],p2); uze=gmul(uze,p3);
! 1822: gptr[0]=&z; gptr[1]=&uze; gerepilemanysp(av,tetpil,gptr,2);
! 1823:
! 1824: av=avma; p1 = gadd(z, gneg(gmul(uze,x))); tetpil = avma;
! 1825: if (!poldivis(p1,y,&p1)) err(bugparier,"subresext");
! 1826: *U=uze; *V=gerepile(av,tetpil,p1); return z;
! 1827: }
! 1828:
! 1829: static GEN
! 1830: scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
! 1831: {
! 1832: long v = varn(x);
! 1833: *U=gzero; *V=gdiv(polun[v],y); return polun[v];
! 1834: }
! 1835:
! 1836: static GEN
! 1837: zero_bezout(GEN y, GEN *U, GEN *V)
! 1838: {
! 1839: GEN x=content(y);
! 1840: *U=gzero; *V = gdiv(polun[varn(y)],x); return gmul(y,*V);
! 1841: }
! 1842:
! 1843: /* calcule U et V tel que Ux+Vy=GCD(x,y) par le sous-resultant */
! 1844: GEN
! 1845: bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
! 1846: {
! 1847: long degq,av,tetpil,tx,ty,dx,dy,du,dv,dr;
! 1848: GEN g,h,r,p1,p2,p3,p4,u,v,lpu,um1,uze,vze,xprim,yprim;
! 1849: GEN *gptr[3];
! 1850:
! 1851: if (gcmp0(x)) return zero_bezout(y,U,V);
! 1852: if (gcmp0(y)) return zero_bezout(x,V,U);
! 1853: tx=typ(x); ty=typ(y); av=avma;
! 1854: if (is_scalar_t(tx) || is_scalar_t(ty))
! 1855: {
! 1856: if (tx==t_POL) return scalar_bezout(x,y,U,V);
! 1857: if (ty==t_POL) return scalar_bezout(y,x,V,U);
! 1858: *U=ginv(x); *V=gzero; return polun[0];
! 1859: }
! 1860: if (tx!=t_POL || ty!=t_POL) err(typeer,"bezoutpol");
! 1861: if (varn(x)!=varn(y))
! 1862: return (varn(x)<varn(y))? scalar_bezout(x,y,U,V)
! 1863: : scalar_bezout(y,x,V,U);
! 1864: dx=lgef(x); dy=lgef(y);
! 1865: if (dx<dy)
! 1866: {
! 1867: GEN *W=U; U=V; V=W;
! 1868: p1=x; x=y; y=p1; du=dx; dx=dy; dy=du;
! 1869: }
! 1870: if (dy==3) return scalar_bezout(x,y,U,V);
! 1871:
! 1872: p3=content(x); u=gdiv(x,p3);
! 1873: p4=content(y); v=gdiv(y,p4);
! 1874: xprim=u; yprim=v; g=gun; h=gun; um1=gun; uze=gzero;
! 1875: for(;;)
! 1876: {
! 1877: du=lgef(u); dv=lgef(v); degq=du-dv;
! 1878: lpu=gpuigs((GEN)v[dv-1],degq+1);
! 1879: p1=gmul(lpu,u); p2=poldivres(p1,v,&r);
! 1880: dr=lgef(r); if (dr<=2) break;
! 1881: p1=gsub(gmul(lpu,um1),gmul(p2,uze)); um1=uze; uze=p1;
! 1882:
! 1883: u=v; p1 = g; g = leading_term(u);
! 1884: switch(degq)
! 1885: {
! 1886: case 0: break;
! 1887: case 1:
! 1888: p1 = gmul(h,p1); h = g; break;
! 1889: default:
! 1890: p1 = gmul(gpuigs(h,degq), p1);
! 1891: h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
! 1892: }
! 1893: v=gdiv(r,p1); uze=gdiv(uze,p1);
! 1894: if (dr==3) break;
! 1895: }
! 1896: if (!poldivis(gsub(v,gmul(uze,xprim)),yprim, &vze))
! 1897: err(warner,"non-exact computation in bezoutpol");
! 1898: uze=gdiv(uze,p3); vze=gdiv(vze,p4); p1=ginv(content(v));
! 1899:
! 1900: tetpil=avma; uze=gmul(uze,p1); vze=gmul(vze,p1); p1=gmul(v,p1);
! 1901: gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&p1;
! 1902: gerepilemanysp(av,tetpil,gptr,3);
! 1903: *U=uze; *V=vze; return p1;
! 1904: }
! 1905:
! 1906:
! 1907:
! 1908: /*******************************************************************/
! 1909: /* */
! 1910: /* RESULTANT PAR L'ALGORITHME DE DUCOS */
! 1911: /* */
! 1912: /*******************************************************************/
! 1913: GEN addshiftw(GEN x, GEN y, long d);
! 1914:
! 1915: static GEN
! 1916: reductum(GEN P)
! 1917: {
! 1918: if (signe(P)==0) return P;
! 1919: return normalizepol_i(dummycopy(P),lgef(P)-1);
! 1920: }
! 1921:
! 1922: static GEN
! 1923: Lazard(GEN x, GEN y, long n)
! 1924: {
! 1925: long a, b;
! 1926: GEN c;
! 1927:
! 1928: if (n==1) return x;
! 1929: a=1; while (n >= (b=a+a)) a=b;
! 1930: c=x; n-=a;
! 1931: while (a>1)
! 1932: {
! 1933: a>>=1; c=gdivexact(gsqr(c),y);
! 1934: if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
! 1935: }
! 1936: return c;
! 1937: }
! 1938:
! 1939: static GEN
! 1940: Lazard2(GEN F, GEN x, GEN y, long n)
! 1941: {
! 1942: if (n<=1) return F;
! 1943: x = Lazard(x,y,n-1); return gdivexact(gmul(x,F),y);
! 1944: }
! 1945:
! 1946: static GEN
! 1947: addshift(GEN x, GEN y)
! 1948: {
! 1949: long v = varn(x);
! 1950: if (!signe(x)) return y;
! 1951: x = addshiftw(x,y,1);
! 1952: setvarn(x,v); return x;
! 1953: }
! 1954:
! 1955: static GEN
! 1956: nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
! 1957: {
! 1958: GEN p0, q0, z0, H, A;
! 1959: long p, q, j, av, lim, v = varn(P);
! 1960:
! 1961: z0 = leading_term(Z);
! 1962: p = degree(P); p0 = leading_term(P); P = reductum(P);
! 1963: q = degree(Q); q0 = leading_term(Q); Q = reductum(Q);
! 1964:
! 1965: av = avma; lim = stack_lim(av,1);
! 1966: H = gneg(reductum(Z));
! 1967: A = gmul((GEN)P[q+2],H);
! 1968: for (j = q+1; j < p; j++)
! 1969: {
! 1970: H = (degree(H) == q-1) ?
! 1971: addshift(reductum(H), gdivexact(gmul(gneg((GEN)H[q+1]),Q), q0)) :
! 1972: addshift(H, zeropol(v));
! 1973: A = gadd(A,gmul((GEN)P[j+2],H));
! 1974: if (low_stack(lim,stack_lim(av,1)))
! 1975: {
! 1976: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&H;
! 1977: if(DEBUGMEM>1) err(warnmem,"nextSousResultant j = %ld/%ld",j,p);
! 1978: gerepilemany(av,gptr,2);
! 1979: }
! 1980: }
! 1981: P = normalizepol_i(P, q+2);
! 1982: A = gdivexact(gadd(A,gmul(z0,P)), p0);
! 1983: A = (degree(H) == q-1) ?
! 1984: gadd(gmul(q0,addshift(reductum(H),A)), gmul(gneg((GEN)H[q+1]),Q)) :
! 1985: gmul(q0, addshift(H,A));
! 1986: return gdivexact(A, ((p-q)&1)? s: gneg(s));
! 1987: }
! 1988:
! 1989: GEN
! 1990: resultantducos(GEN P, GEN Q)
! 1991: {
! 1992: long delta, av=avma, tetpil, lim = stack_lim(av,1);
! 1993: GEN Z, s;
! 1994:
! 1995: if ((Z = init_resultant(P,Q))) return Z;
! 1996: delta = degree(P) - degree(Q);
! 1997: if (delta < 0)
! 1998: {
! 1999: Z = ((degree(P)&1) && (degree(Q)&1)) ? gneg(Q) : Q;
! 2000: Q = P; P = Z; delta = -delta;
! 2001: }
! 2002: s = gun;
! 2003: if (degree(Q) > 0)
! 2004: {
! 2005: s = gpuigs(leading_term(Q),delta);
! 2006: Z = Q;
! 2007: Q = pseudorem(P, gneg(Q));
! 2008: P = Z;
! 2009: while(degree(Q) > 0)
! 2010: {
! 2011: if (low_stack(lim,stack_lim(av,1)))
! 2012: {
! 2013: GEN *gptr[2]; gptr[0]=&P; gptr[1]=&Q;
! 2014: if(DEBUGMEM>1) err(warnmem,"resultantducos, deg Q = %ld",degree(Q));
! 2015: gerepilemany(av,gptr,2); s = leading_term(P);
! 2016: }
! 2017: delta = degree(P) - degree(Q);
! 2018: Z = Lazard2(Q, leading_term(Q), s, delta);
! 2019: Q = nextSousResultant(P, Q, Z, s);
! 2020: P = Z;
! 2021: s = leading_term(P);
! 2022: }
! 2023: }
! 2024: if (!signe(Q)) { avma = av; return gzero; }
! 2025: s = Lazard(leading_term(Q), s, degree(P));
! 2026: tetpil = avma; return gerepile(av,tetpil,gcopy(s));
! 2027: }
! 2028:
! 2029: /*******************************************************************/
! 2030: /* */
! 2031: /* RESULTANT PAR MATRICE DE SYLVESTER */
! 2032: /* */
! 2033: /*******************************************************************/
! 2034: GEN
! 2035: sylvestermatrix_i(GEN x, GEN y)
! 2036: {
! 2037: long d,dx,dy,i,j;
! 2038: GEN p1,p2;
! 2039:
! 2040: dx=lgef(x)-3; dy=lgef(y)-3; d=dx+dy;
! 2041: p1=cgetg(d+1,t_MAT);
! 2042: for (j=1; j<=dy; j++)
! 2043: {
! 2044: p2=cgetg(d+1,t_COL); p1[j]=(long)p2;
! 2045: for (i=1; i<j; i++) p2[i]=zero;
! 2046: for ( ; i<=j+dx; i++) p2[i]=x[dx-i+j+2];
! 2047: for ( ; i<=d; i++) p2[i]=zero;
! 2048: }
! 2049: for (j=1; j<=dx; j++)
! 2050: {
! 2051: p2=cgetg(d+1,t_COL); p1[j+dy]=(long)p2;
! 2052: for (i=1; i<j; i++) p2[i]=zero;
! 2053: for ( ; i<=j+dy; i++) p2[i]=y[dy-i+j+2];
! 2054: for ( ; i<=d; i++) p2[i]=zero;
! 2055: }
! 2056: return p1;
! 2057: }
! 2058:
! 2059: GEN
! 2060: sylvestermatrix(GEN x, GEN y)
! 2061: {
! 2062: long i,j,lx;
! 2063: if (typ(x)!=t_POL || typ(y)!=t_POL) err(typeer,"sylvestermatrix");
! 2064: if (varn(x) != varn(y))
! 2065: err(talker,"not the same variables in sylvestermatrix");
! 2066: x = sylvestermatrix_i(x,y); lx = lg(x);
! 2067: for (i=1; i<lx; i++)
! 2068: for (j=1; j<lx; j++) coeff(x,i,j) = lcopy(gcoeff(x,i,j));
! 2069: return x;
! 2070: }
! 2071:
! 2072: GEN
! 2073: resultant2(GEN x, GEN y)
! 2074: {
! 2075: long av;
! 2076: GEN r;
! 2077: if ((r = init_resultant(x,y))) return r;
! 2078: av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
! 2079: }
! 2080:
! 2081: static GEN
! 2082: fix_pol(GEN x, long v, long *mx)
! 2083: {
! 2084: long vx;
! 2085: GEN p1;
! 2086:
! 2087: if (typ(x) == t_POL)
! 2088: {
! 2089: vx = varn(x);
! 2090: if (vx)
! 2091: {
! 2092: if (v>=vx) return gsubst(x,v,polx[0]);
! 2093: p1 = cgetg(3,t_POL);
! 2094: p1[1] = evalvarn(0)|evalsigne(signe(x))|evallgef(3);
! 2095: p1[2] = (long)x; return p1;
! 2096: }
! 2097: if (v)
! 2098: {
! 2099: *mx = 1;
! 2100: return gsubst(gsubst(x,0,polx[MAXVARN]),v,polx[0]);
! 2101: }
! 2102: }
! 2103: return x;
! 2104: }
! 2105:
! 2106: /* resultant of x and y with respect to variable v, or with respect to their
! 2107: * main variable if v < 0.
! 2108: */
! 2109: GEN
! 2110: polresultant0(GEN x, GEN y, long v, long flag)
! 2111: {
! 2112: long av = avma, m = 0;
! 2113:
! 2114: if (v >= 0)
! 2115: {
! 2116: x = fix_pol(x,v, &m);
! 2117: y = fix_pol(y,v, &m);
! 2118: }
! 2119: switch(flag)
! 2120: {
! 2121: case 0: x=subresall(x,y,NULL); break;
! 2122: case 1: x=resultant2(x,y); break;
! 2123: case 2: x=resultantducos(x,y); break;
! 2124: default: err(flagerr,"polresultant");
! 2125: }
! 2126: if (m) x = gsubst(x,MAXVARN,polx[0]);
! 2127: return gerepileupto(av,x);
! 2128: }
! 2129:
! 2130: /*******************************************************************/
! 2131: /* */
! 2132: /* P.G.C.D PAR L'ALGORITHME DU SOUS RESULTANT */
! 2133: /* */
! 2134: /*******************************************************************/
! 2135:
! 2136: GEN
! 2137: srgcd(GEN x, GEN y)
! 2138: {
! 2139: long av,av1,tetpil,tx=typ(x),ty=typ(y),dx,dy,vx,lim;
! 2140: GEN d,g,h,p1,p2,u,v;
! 2141:
! 2142: if (!signe(y)) return gcopy(x);
! 2143: if (!signe(x)) return gcopy(y);
! 2144: if (is_scalar_t(tx) || is_scalar_t(ty)) return gun;
! 2145: if (tx!=t_POL || ty!=t_POL) err(typeer,"srgcd");
! 2146: vx=varn(x); if (vx!=varn(y)) return gun;
! 2147: if (ismonome(x)) return gcdmonome(x,y);
! 2148: if (ismonome(y)) return gcdmonome(y,x);
! 2149: if (isrational(x) && isrational(y)) return modulargcd(x,y);
! 2150:
! 2151: av=avma;
! 2152: if (issimplefield(x) || issimplefield(y)) x = polgcdnun(x,y);
! 2153: else
! 2154: {
! 2155: dx=lgef(x); dy=lgef(y);
! 2156: if (dx<dy) { p1=x; x=y; y=p1; tx=dx; dx=dy; dy=tx; }
! 2157: p1=content(x); p2=content(y); d=ggcd(p1,p2);
! 2158:
! 2159: tetpil=avma; d=gmul(d,polun[vx]);
! 2160: if (dy==3) return gerepile(av,tetpil,d);
! 2161:
! 2162: av1=avma; lim=stack_lim(av1,1);
! 2163: u=gdiv(x,p1); v=gdiv(y,p2); g=h=gun;
! 2164: for(;;)
! 2165: {
! 2166: GEN r = pseudorem(u,v);
! 2167: long degq, du, dv, dr=lgef(r);
! 2168:
! 2169: if (dr<=3)
! 2170: {
! 2171: if (gcmp0(r)) break;
! 2172: avma=av1; return gerepile(av,tetpil,d);
! 2173: }
! 2174: if (DEBUGLEVEL > 9) fprintferr("srgcd: dr = %ld\n", dr);
! 2175: du=lgef(u); dv=lgef(v); degq=du-dv; u=v;
! 2176: switch(degq)
! 2177: {
! 2178: case 0:
! 2179: v = gdiv(r,g);
! 2180: g = leading_term(u);
! 2181: break;
! 2182: case 1:
! 2183: v = gdiv(r,gmul(h,g));
! 2184: h = g = leading_term(u);
! 2185: break;
! 2186: default:
! 2187: v = gdiv(r,gmul(gpuigs(h,degq),g));
! 2188: g = leading_term(u);
! 2189: h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
! 2190: }
! 2191: if (low_stack(lim, stack_lim(av1,1)))
! 2192: {
! 2193: GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
! 2194: if(DEBUGMEM>1) err(warnmem,"srgcd");
! 2195: gerepilemany(av1,gptr,4);
! 2196: }
! 2197: }
! 2198: p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);
! 2199: x = gmul(d,v);
! 2200: }
! 2201:
! 2202: if (typ(x)!=t_POL) x = gmul(polun[vx],x);
! 2203: else
! 2204: {
! 2205: p1=leading_term(x); ty=typ(p1);
! 2206: if ((is_frac_t(ty) || is_intreal_t(ty)) && gsigne(p1)<0) x = gneg(x);
! 2207: }
! 2208: return gerepileupto(av,x);
! 2209: }
! 2210:
! 2211: GEN qf_disc(GEN x, GEN y, GEN z);
! 2212:
! 2213: GEN
! 2214: poldisc0(GEN x, long v)
! 2215: {
! 2216: long av,tx=typ(x),i;
! 2217: GEN z,p1,p2;
! 2218:
! 2219: switch(tx)
! 2220: {
! 2221: case t_POL:
! 2222: if (gcmp0(x)) return gzero;
! 2223: av = avma; i = 0;
! 2224: if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
! 2225: p1 = subres(x, derivpol(x));
! 2226: p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
! 2227: if ((lgef(x)-3) & 2) p1 = gneg(p1);
! 2228: if (i) p1 = gsubst(p1, MAXVARN, polx[0]);
! 2229: return gerepileupto(av,p1);
! 2230:
! 2231: case t_COMPLEX:
! 2232: return stoi(-4);
! 2233:
! 2234: case t_QUAD: case t_POLMOD:
! 2235: return poldisc0((GEN)x[1], v);
! 2236:
! 2237: case t_QFR: case t_QFI:
! 2238: av = avma; return gerepileuptoint(av, qf_disc(x,NULL,NULL));
! 2239:
! 2240: case t_VEC: case t_COL: case t_MAT:
! 2241: i=lg(x); z=cgetg(i,tx);
! 2242: for (i--; i; i--) z[i]=(long)poldisc0((GEN)x[i], v);
! 2243: return z;
! 2244: }
! 2245: err(typeer,"discsr");
! 2246: return NULL; /* not reached */
! 2247: }
! 2248:
! 2249: GEN
! 2250: discsr(GEN x)
! 2251: {
! 2252: return poldisc0(x, -1);
! 2253: }
! 2254:
! 2255: GEN
! 2256: reduceddiscsmith(GEN pol)
! 2257: {
! 2258: long av=avma,tetpil,i,j,n;
! 2259: GEN polp,alpha,p1,m;
! 2260:
! 2261: if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith");
! 2262: n=lgef(pol)-3;
! 2263: if (n<=0) err(constpoler,"reduceddiscsmith");
! 2264: check_pol_int(pol);
! 2265: if (!gcmp1((GEN)pol[n+2]))
! 2266: err(talker,"non-monic polynomial in poldiscreduced");
! 2267: polp = derivpol(pol); alpha = polx[varn(pol)];
! 2268: m=cgetg(n+1,t_MAT);
! 2269: for (j=1; j<=n; j++)
! 2270: {
! 2271: p1=cgetg(n+1,t_COL); m[j]=(long)p1;
! 2272: for (i=1; i<=n; i++) p1[i]=(long)truecoeff(polp,i-1);
! 2273: if (j<n) polp = gres(gmul(alpha,polp), pol);
! 2274: }
! 2275: tetpil=avma; return gerepile(av,tetpil,smith(m));
! 2276: }
! 2277:
! 2278: /***********************************************************************/
! 2279: /** **/
! 2280: /** ALGORITHME DE STURM **/
! 2281: /** (number of real roots of x in ]a,b]) **/
! 2282: /** **/
! 2283: /***********************************************************************/
! 2284:
! 2285: /* if a (resp. b) is NULL, set it to -oo (resp. +oo) */
! 2286: long
! 2287: sturmpart(GEN x, GEN a, GEN b)
! 2288: {
! 2289: long av = avma,sl,sr,s,t,r1;
! 2290: GEN g,h,u,v;
! 2291:
! 2292: if (typ(x) != t_POL) err(typeer,"sturm");
! 2293: if (gcmp0(x)) err(zeropoler,"sturm");
! 2294: s=lgef(x); if (s==3) return 0;
! 2295:
! 2296: sl = gsigne(leading_term(x));
! 2297: if (s==4)
! 2298: {
! 2299: s = b? gsigne(poleval(x,b)): sl;
! 2300: t = a? gsigne(poleval(x,a)): -sl;
! 2301: avma = av; return (s == t)? 0: 1;
! 2302: }
! 2303: u=gdiv(x,content(x)); v=derivpol(x);
! 2304: v=gdiv(v,content(v));
! 2305: g=gun; h=gun;
! 2306: s = b? gsigne(poleval(u,b)): sl;
! 2307: t = a? gsigne(poleval(u,a)): ((lgef(u)&1)? sl: -sl);
! 2308: r1=0;
! 2309: sr = b? gsigne(poleval(v,b)): s;
! 2310: if (sr)
! 2311: {
! 2312: if (!s) s=sr;
! 2313: else if (sr!=s) { s= -s; r1--; }
! 2314: }
! 2315: sr = a? gsigne(poleval(v,a)): -t;
! 2316: if (sr)
! 2317: {
! 2318: if (!t) t=sr;
! 2319: else if (sr!=t) { t= -t; r1++; }
! 2320: }
! 2321: for(;;)
! 2322: {
! 2323: GEN p1, r = pseudorem(u,v);
! 2324: long du=lgef(u), dv=lgef(v), dr=lgef(r), degq=du-dv;
! 2325:
! 2326: if (dr<=2) err(talker,"not a squarefree polynomial in sturm");
! 2327: if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
! 2328: sl = gsigne((GEN) r[dr-1]);
! 2329: sr = b? gsigne(poleval(r,b)): sl;
! 2330: if (sr)
! 2331: {
! 2332: if (!s) s=sr;
! 2333: else if (sr!=s) { s= -s; r1--; }
! 2334: }
! 2335: sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
! 2336: if (sr)
! 2337: {
! 2338: if (!t) t=sr;
! 2339: else if (sr!=t) { t= -t; r1++; }
! 2340: }
! 2341: if (dr==3) { avma=av; return r1; }
! 2342:
! 2343: u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
! 2344: switch(degq)
! 2345: {
! 2346: case 0: break;
! 2347: case 1:
! 2348: p1 = gmul(h,p1); h = g; break;
! 2349: default:
! 2350: p1 = gmul(gpuigs(h,degq),p1);
! 2351: h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
! 2352: }
! 2353: v = gdiv(r,p1);
! 2354: }
! 2355: }
! 2356:
! 2357: /*******************************************************************/
! 2358: /* */
! 2359: /* POLYNOME QUADRATIQUE ASSOCIE A UN DISCRIMINANT */
! 2360: /* */
! 2361: /*******************************************************************/
! 2362:
! 2363: GEN
! 2364: quadpoly0(GEN x, long v)
! 2365: {
! 2366: long res,l,tetpil,i,sx, tx = typ(x);
! 2367: GEN y,p1;
! 2368:
! 2369: if (is_matvec_t(tx))
! 2370: {
! 2371: l=lg(x); y=cgetg(l,tx);
! 2372: for (i=1; i<l; i++) y[i]=(long)quadpoly0((GEN)x[i],v);
! 2373: return y;
! 2374: }
! 2375: if (tx!=t_INT) err(arither1);
! 2376: if (v < 0) v = 0;
! 2377: sx=signe(x);
! 2378: if (!sx) err(talker,"zero discriminant in quadpoly");
! 2379: y=cgetg(5,t_POL);
! 2380: y[1]=evalsigne(1) | evalvarn(v) | evallgef(5); y[4]=un;
! 2381: res=mod4(x); if (sx<0 && res) res=4-res;
! 2382: if (res>1) err(funder2,"quadpoly");
! 2383:
! 2384: l=avma; p1=shifti(x,-2); setsigne(p1,-signe(p1));
! 2385: y[2] = (long) p1;
! 2386: if (!res) y[3] = zero;
! 2387: else
! 2388: {
! 2389: if (sx<0) { tetpil=avma; y[2] = lpile(l,tetpil,addsi(1,p1)); }
! 2390: y[3] = lnegi(gun);
! 2391: }
! 2392: return y;
! 2393: }
! 2394:
! 2395: GEN
! 2396: quadpoly(GEN x)
! 2397: {
! 2398: return quadpoly0(x,-1);
! 2399: }
! 2400:
! 2401: GEN
! 2402: quadgen(GEN x)
! 2403: {
! 2404: GEN y=cgetg(4,t_QUAD);
! 2405: y[1]=lquadpoly(x); y[2]=zero; y[3]=un; return y;
! 2406: }
! 2407:
! 2408: /*******************************************************************/
! 2409: /* */
! 2410: /* INVERSE MODULO GENERAL */
! 2411: /* */
! 2412: /*******************************************************************/
! 2413:
! 2414: GEN
! 2415: ginvmod(GEN x, GEN y)
! 2416: {
! 2417: long tx=typ(x);
! 2418:
! 2419: switch(typ(y))
! 2420: {
! 2421: case t_POL:
! 2422: if (tx==t_POL) return polinvmod(x,y);
! 2423: if (is_scalar_t(tx)) return ginv(x);
! 2424: break;
! 2425: case t_INT:
! 2426: if (tx==t_INT) return mpinvmod(x,y);
! 2427: if (tx==t_POL) return gzero;
! 2428: }
! 2429: err(typeer,"ginvmod");
! 2430: return NULL; /* not reached */
! 2431: }
! 2432:
! 2433: /***********************************************************************/
! 2434: /** **/
! 2435: /** NEWTON POLYGON **/
! 2436: /** **/
! 2437: /***********************************************************************/
! 2438:
! 2439: /* assume leading coeff of x is non-zero */
! 2440: GEN
! 2441: newtonpoly(GEN x, GEN p)
! 2442: {
! 2443: GEN y;
! 2444: long n,ind,a,b,c,u1,u2,r1,r2;
! 2445: long *vval, num[] = {evaltyp(t_INT) | m_evallg(3), 0, 0};
! 2446:
! 2447: if (typ(x)!=t_POL) err(typeer,"newtonpoly");
! 2448: n=lgef(x)-3; if (n<=0) { y=cgetg(1,t_VEC); return y; }
! 2449: y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
! 2450: vval = (long *) gpmalloc(sizeof(long)*(n+1));
! 2451: for (a=0; a<=n; a++) vval[a] = ggval((GEN)x[a],p);
! 2452: for (a=0, ind=1; a<n; a++)
! 2453: {
! 2454: if (vval[a] != VERYBIGINT) break;
! 2455: y[ind++] = lstoi(VERYBIGINT);
! 2456: }
! 2457: for (b=a+1; b<=n; a=b, b=a+1)
! 2458: {
! 2459: while (vval[b] == VERYBIGINT) b++;
! 2460: u1=vval[a]-vval[b]; u2=b-a;
! 2461: for (c=b+1; c<=n; c++)
! 2462: {
! 2463: if (vval[c] == VERYBIGINT) continue;
! 2464: r1=vval[a]-vval[c]; r2=c-a;
! 2465: if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
! 2466: }
! 2467: while (ind<=b) { affsi(u1,num); y[ind++] = ldivgs(num,u2); }
! 2468: }
! 2469: free(vval); return y;
! 2470: }
! 2471:
! 2472: int cmp_pol(GEN x, GEN y);
! 2473:
! 2474: /* Factor polynomial a on the number field defined by polynomial t */
! 2475: GEN
! 2476: polfnf(GEN a, GEN t)
! 2477: {
! 2478: GEN x0, y,p1,p2,u,g,fa,n,unt;
! 2479: long av=avma,tetpil,lx,v,i,k,vt;
! 2480:
! 2481: if (typ(a)!=t_POL || typ(t)!=t_POL) err(typeer,"polfnf");
! 2482: if (gcmp0(a)) return gcopy(a);
! 2483: vt=varn(t); v=varn(a);
! 2484: if (vt<=v)
! 2485: err(talker,"polynomial variable must be of higher priority than number field variable\nin factornf");
! 2486: u = gdiv(a, ggcd(a,derivpol(a)));
! 2487: unt = gmodulsg(1,t); u = gmul(unt,u);
! 2488: g = lift(u);
! 2489: for (k=-1; ; k++)
! 2490: {
! 2491: n = gsub(polx[MAXVARN], gmulsg(k,polx[vt]));
! 2492: n = subres(t, poleval(g,n));
! 2493: if (lgef(ggcd(n,derivpol(n))) == 3) break;
! 2494: }
! 2495: fa=factor(n); fa=(GEN)fa[1]; lx=lg(fa);
! 2496: y=cgetg(3,t_MAT);
! 2497: p1=cgetg(lx,t_COL); y[1]=(long)p1;
! 2498: p2=cgetg(lx,t_COL); y[2]=(long)p2;
! 2499: x0 = gadd(polx[v],gmulsg(k,gmodulcp(polx[vt],t)));
! 2500: for (i=1; i<lx; i++)
! 2501: {
! 2502: GEN b, pro = ggcd(u, gmul(unt, poleval((GEN)fa[i], x0)));
! 2503: long e;
! 2504:
! 2505: p1[i] = (typ(pro)==t_POL)? ldiv(pro,leading_term(pro)): (long)pro;
! 2506: if (gcmp1((GEN)p1[i])) err(talker,"reducible modulus in factornf");
! 2507: e=0; while (poldivis(a,(GEN)p1[i], &b)) { a = b; e++; }
! 2508: p2[i] = lstoi(e);
! 2509: }
! 2510: (void)sort_factor(y, cmp_pol);
! 2511: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
! 2512: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to polarit2.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