Annotation of OpenXM_contrib/pari/src/basemath/polarit2.c, Revision 1.1.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