Annotation of OpenXM_contrib/pari-2.2/src/basemath/polarit2.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: polarit2.c,v 1.120 2001/10/01 12:11:31 karim Exp $
2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /***********************************************************************/
17: /** **/
18: /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
19: /** (second part) **/
20: /** **/
21: /***********************************************************************/
22: #include "pari.h"
23:
24: #define swap(a,b) { GEN _x = a; a = b; b = _x; }
25: #define lswap(a,b) { long _x = a; a = b; b = _x; }
26: #define pswap(a,b) { GEN *_x = a; a = b; b = _x; }
27: #define both_even(a,b) ((((a)|(b))&1) == 0)
28:
29: GEN gassoc_proto(GEN f(GEN,GEN),GEN,GEN);
30:
31: GEN
32: polsym(GEN x, long n)
33: {
34: long av1,av2,dx=degpol(x),i,k;
35: GEN s,y,x_lead;
36:
37: if (n<0) err(impl,"polsym of a negative n");
38: if (typ(x) != t_POL) err(typeer,"polsym");
39: if (!signe(x)) err(zeropoler,"polsym");
40: y=cgetg(n+2,t_COL); y[1]=lstoi(dx);
41: x_lead=(GEN)x[dx+2]; if (gcmp1(x_lead)) x_lead=NULL;
42: for (k=1; k<=n; k++)
43: {
44: av1=avma; s = (dx>=k)? gmulsg(k,(GEN)x[dx+2-k]): gzero;
45: for (i=1; i<k && i<=dx; i++)
46: s = gadd(s,gmul((GEN)y[k-i+1],(GEN)x[dx+2-i]));
47: if (x_lead) s = gdiv(s,x_lead);
48: av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
49: }
50: return y;
51: }
52:
53: static int (*polcmp_coeff_cmp)(GEN,GEN);
54:
55: /* assume x and y are polynomials in the same variable whose coeffs can be
56: * compared (used to sort polynomial factorizations)
57: */
58: static int
59: polcmp(GEN x, GEN y)
60: {
61: long i, lx = lgef(x), ly = lgef(y);
62: int fl;
63: #if 0 /* for efficiency */
64: if (typ(x) != t_POL || typ(y) != t_POL)
65: err(talker,"not a polynomial in polcmp");
66: #endif
67: if (lx > ly) return 1;
68: if (lx < ly) return -1;
69: for (i=lx-1; i>1; i--)
70: if ((fl = polcmp_coeff_cmp((GEN)x[i], (GEN)y[i]))) return fl;
71: return 0;
72: }
73:
74: /* sort polynomial factorization so that factors appear by decreasing
75: * degree, sorting coefficients according to cmp. In place */
76: GEN
77: sort_factor(GEN y, int (*cmp)(GEN,GEN))
78: {
79: int (*old)(GEN,GEN) = polcmp_coeff_cmp;
80: polcmp_coeff_cmp = cmp;
81: (void)sort_factor_gen(y,polcmp);
82: polcmp_coeff_cmp = old; return y;
83: }
84:
85: /* sort generic factorisation */
86: GEN
87: sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))
88: {
89: long n,i, av = avma;
90: GEN a,b,A,B,w;
91: a = (GEN)y[1]; n = lg(a); A = new_chunk(n);
92: b = (GEN)y[2]; B = new_chunk(n);
93: w = gen_sort(a, cmp_IND | cmp_C, cmp);
94: for (i=1; i<n; i++) { A[i] = a[w[i]]; B[i] = b[w[i]]; }
95: for (i=1; i<n; i++) { a[i] = A[i]; b[i] = B[i]; }
96: avma = av; return y;
97: }
98:
99: /* for internal use */
100: GEN
101: centermod_i(GEN x, GEN p, GEN ps2)
102: {
103: long av,i,lx;
104: GEN y,p1;
105:
106: if (!ps2) ps2 = shifti(p,-1);
107: switch(typ(x))
108: {
109: case t_INT:
110: y=modii(x,p);
111: if (cmpii(y,ps2)>0) return subii(y,p);
112: return y;
113:
114: case t_POL: lx=lgef(x);
115: y=cgetg(lx,t_POL); y[1]=x[1];
116: for (i=2; i<lx; i++)
117: {
118: av=avma; p1=modii((GEN)x[i],p);
119: if (cmpii(p1,ps2)>0) p1=subii(p1,p);
120: y[i]=lpileupto(av,p1);
121: }
122: return y;
123:
124: case t_COL: lx=lg(x);
125: y=cgetg(lx,t_COL);
126: for (i=1; i<lx; i++)
127: {
128: p1=modii((GEN)x[i],p);
129: if (cmpii(p1,ps2)>0) p1=subii(p1,p);
130: y[i]=(long)p1;
131: }
132: return y;
133: }
134: return x;
135: }
136:
137: GEN
138: centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
139:
140: /* assume same variables, y normalized, non constant */
141: static GEN
142: polidivis(GEN x, GEN y, GEN bound)
143: {
144: long dx,dy,dz,i,j,av, vx = varn(x);
145: GEN z,p1,y_lead;
146:
147: dy=degpol(y);
148: dx=degpol(x);
149: dz=dx-dy; if (dz<0) return NULL;
150: z=cgetg(dz+3,t_POL);
151: x += 2; y += 2; z += 2;
152: y_lead = (GEN)y[dy];
153: if (gcmp1(y_lead)) y_lead = NULL;
154:
155: p1 = (GEN)x[dx];
156: z[dz]=y_lead? ldivii(p1,y_lead): licopy(p1);
157: for (i=dx-1; i>=dy; i--)
158: {
159: av=avma; p1=(GEN)x[i];
160: for (j=i-dy+1; j<=i && j<=dz; j++)
161: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
162: if (y_lead) { p1 = gdiv(p1,y_lead); if (typ(p1)!=t_INT) return NULL; }
163: if (absi_cmp(p1, bound) > 0) return NULL;
164: p1 = gerepileupto(av,p1);
165: z[i-dy] = (long)p1;
166: }
167: av = avma;
168: for (;; i--)
169: {
170: p1 = (GEN)x[i];
171: /* we always enter this loop at least once */
172: for (j=0; j<=i && j<=dz; j++)
173: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
174: if (!gcmp0(p1)) return NULL;
175: avma = av;
176: if (!i) break;
177: }
178: z -= 2;
179: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
180: return z;
181: }
182:
183: static long
184: min_deg(long jmax,ulong tabbit[])
185: {
186: long j, k = 1, j1 = 2;
187:
188: for (j=jmax; j>=0; j--)
189: {
190: for ( ; k<=15; k++)
191: {
192: if (tabbit[j]&j1) return ((jmax-j)<<4)+k;
193: j1<<=1;
194: }
195: k = 0; j1 = 1;
196: }
197: return (jmax<<4)+15;
198: }
199:
200: /* tabkbit is a bit vector (only lowest 32 bits of each word are used
201: * on 64bit architecture): reading from right to left, bit i+1 is set iff
202: * degree i is attainable from the factorisation mod p.
203: *
204: * record N modular factors of degree d. */
205: static void
206: record_factors(long N, long d, long jmax, ulong *tabkbit, ulong *tmp)
207: {
208: long n,j, a = d>>4, b = d&15; /* d = 16 a + b */
209: ulong *tmp2 = tmp - a;
210:
211: for (n = 1; n <= N; n++)
212: {
213: ulong pro, rem = 0;
214: for (j=jmax; j>=a; j--)
215: {
216: pro = tabkbit[j] << b;
217: tmp2[j] = (pro&0xffff) + rem; rem = pro >> 16;
218: }
219: for (j=jmax-a; j>=0; j--) tabkbit[j] |= tmp[j];
220: }
221: }
222:
223: /***********************************************************************/
224: /** **/
225: /** QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL) **/
226: /** **/
227: /***********************************************************************/
228: /* Setup for divide/conquer quadratic Hensel lift
229: * a = set of k t_POL in Z[X] corresponding to factors over Fp
230: * V = set of products of factors built as follows
231: * 1) V[1..k] = initial a
232: * 2) iterate:
233: * append to V the two smallest factors (minimal degree) in a, remove them
234: * from a and replace them by their product [net loss for a = 1 factor]
235: *
236: * W = bezout coeffs W[i]V[i] + W[i+1]V[i+1] = 1
237: *
238: * link[i] = -j if V[i] = a[j]
239: * j if V[i] = V[j] * V[j+1]
240: * Arrays (link, V, W) pre-allocated for 2k - 2 elements */
241: #if 0 /* restricted to Fp */
242: static void
243: BuildTree(GEN link, GEN V, GEN W, GEN a, GEN p)
244: {
245: long k = lg(a)-1;
246: long i, j, s, minp, mind;
247:
248: for (i=1; i<=k; i++) { V[i] = a[i]; link[i] = -i; }
249:
250: for (j=1; j <= 2*k-5; j+=2,i++)
251: {
252: minp = j;
253: mind = degpol(V[j]);
254: for (s=j+1; s<i; s++)
255: if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
256:
257: lswap(V[j], V[minp]);
258: lswap(link[j], link[minp]);
259:
260: minp = j+1;
261: mind = degpol(V[j+1]);
262: for (s=j+2; s<i; s++)
263: if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
264:
265: lswap(V[j+1], V[minp]);
266: lswap(link[j+1], link[minp]);
267:
268: V[i] = (long)FpX_mul((GEN)V[j], (GEN)V[j+1], p);
269: link[i] = j;
270: }
271:
272: for (j=1; j <= 2*k-3; j+=2)
273: {
274: GEN d, u, v;
275: d = FpX_extgcd((GEN)V[j], (GEN)V[j+1], p, &u, &v);
276: if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");
277: d = (GEN)d[2];
278: if (!gcmp1(d))
279: {
280: d = mpinvmod(d, p);
281: u = FpX_Fp_mul(u, d, p);
282: v = FpX_Fp_mul(v, d, p);
283: }
284: W[j] = (long)u;
285: W[j+1] = (long)v;
286: }
287: }
288: #endif
289:
290: /* same over Fp[Y] / (T) [copy-paste] */
291: static void
292: BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
293: {
294: long k = lg(a)-1;
295: long i, j, s, minp, mind;
296:
297: for (i=1; i<=k; i++) { V[i] = a[i]; link[i] = -i; }
298:
299: for (j=1; j <= 2*k-5; j+=2,i++)
300: {
301: minp = j;
302: mind = degpol(V[j]);
303: for (s=j+1; s<i; s++)
304: if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
305:
306: lswap(V[j], V[minp]);
307: lswap(link[j], link[minp]);
308:
309: minp = j+1;
310: mind = degpol(V[j+1]);
311: for (s=j+2; s<i; s++)
312: if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
313:
314: lswap(V[j+1], V[minp]);
315: lswap(link[j+1], link[minp]);
316:
317: V[i] = (long)FpXQX_mul((GEN)V[j], (GEN)V[j+1], T, p);
318: link[i] = j;
319: }
320:
321: for (j=1; j <= 2*k-3; j+=2)
322: {
323: GEN d, u, v;
324: d = FpXQX_extgcd((GEN)V[j], (GEN)V[j+1], T, p, &u, &v);
325: if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");
326: d = (GEN)d[2];
327: if (!gcmp1(d))
328: {
329: d = FpXQ_inv(d, T, p);
330: u = FpXQX_FpXQ_mul(u, d, T, p);
331: v = FpXQX_FpXQ_mul(v, d, T, p);
332: }
333: W[j] = (long)u;
334: W[j+1] = (long)v;
335: }
336: }
337:
338: /* au + bv = 1 (p0), ab = f (p0). Lift mod p1 = p0 pd (<= p0^2).
339: * If noinv is set, don't lift the inverses u and v */
340: static void
341: HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, GEN p0, int noinv)
342: {
343: ulong av = avma;
344: long space = lgef(f) * (lgefint(pd) + lgefint(p0) - 2);
345: GEN a2,b2,g,z,s,t;
346: GEN a = (GEN)V[j], b = (GEN)V[j+1];
347: GEN u = (GEN)W[j], v = (GEN)W[j+1];
348:
349: if (T) space *= lgef(T);
350:
351: (void)new_chunk(space); /* HACK */
352: g = gadd(f, gneg_i(gmul(a,b)));
353: if (T) g = FpXQX_red(g, T, mulii(p0,pd));
354: g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);
355: z = FpXQX_mul(v,g, T,pd);
356: t = FpXQX_divres(z,a, T,pd, &s);
357: t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);
358: t = gmul(t,p0);
359: s = gmul(s,p0);
360: avma = av;
361:
362: /* already reduced mod p1 = pd p0 */
363: a2 = gadd(a,s); V[j] = (long)a2;
364: b2 = gadd(b,t); V[j+1] = (long)b2;
365: if (noinv) return;
366:
367: av = avma;
368: (void)new_chunk(space); /* HACK */
369: g = gadd(gmul(u,a2), gmul(v,b2));
370: g = gadd(gneg_i(g), gun);
371:
372: if (T) g = FpXQX_red(g, T, mulii(p0,pd));
373: g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);
374: z = FpXQX_mul(v,g, T,pd);
375: t = FpXQX_divres(z,a, T,pd, &s);
376: t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);
377: t = gmul(t,p0);
378: s = gmul(s,p0);
379: avma = av;
380:
381: u = gadd(u,t); W[j] = (long)u;
382: v = gadd(v,s); W[j+1] = (long)v;
383: }
384:
385: /* v list of factors, w list of inverses. f = v[j] v[j+1]
386: * Lift v[j] and v[j+1] mod p0 pd (possibly mod T), then all their divisors */
387: static void
388: RecTreeLift(GEN link, GEN v, GEN w, GEN T, GEN pd, GEN p0, GEN f, long j, int noinv)
389: {
390: if (j < 0) return;
391:
392: HenselLift(v, w, j, f, T, pd, p0, noinv);
393: RecTreeLift(link, v, w, T, pd, p0, (GEN)v[j] , link[j ], noinv);
394: RecTreeLift(link, v, w, T, pd, p0, (GEN)v[j+1], link[j+1], noinv);
395: }
396:
397: /* lift from p^{e0} to p^{e1} */
398: static void
399: TreeLift(GEN link, GEN v, GEN w, GEN T, GEN p, long e0, long e1, GEN f, int noinv)
400: {
401: GEN p0 = gpowgs(p, e0);
402: GEN pd = gpowgs(p, e1-e0);
403: RecTreeLift(link, v, w, T, pd, p0, f, lg(v)-2, noinv);
404: }
405:
406: /* a = modular factors of f mod (p,T) [possibly T=NULL], lift to precision p^e0
407: * flag = 0: standard.
408: * flag = 1: return TreeLift structure
409: * flag = 2: a = TreeLift structure, go on lifting, as flag = 1 otherwise */
410: static GEN
411: MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, int flag)
412: {
413: long l, i, e = e0, k = lg(a) - 1;
414: GEN E, v, w, link;
415:
416: if (k < 2 || e < 1) err(talker, "MultiLift: bad args");
417: if (e == 1) return a;
418: if (typ(a[1]) == t_INT) flag = 2;
419: else if (flag == 2) flag = 1;
420:
421: E = cgetg(BITS_IN_LONG, t_VECSMALL);
422: l = 0; E[++l] = e;
423: while (e > 1) { e = (e+1)/2; E[++l] = e; }
424:
425: if (DEBUGLEVEL > 3) timer2();
426:
427: if (flag != 2)
428: {
429: v = cgetg(2*k - 2 + 1, t_VEC);
430: w = cgetg(2*k - 2 + 1, t_VEC);
431: link=cgetg(2*k - 2 + 1, t_VECSMALL);
432: BuildTree(link, v, w, a, T, p);
433: if (DEBUGLEVEL > 3) msgtimer("building tree");
434: }
435: else
436: {
437: e = itos((GEN)a[1]);
438: link = (GEN)a[2];
439: v = (GEN)a[3];
440: w = (GEN)a[4];
441: }
442:
443: for (i = l; i > 1; i--) {
444: if (E[i-1] < e) continue;
445: TreeLift(link, v, w, T, p, E[i], E[i-1], f, (flag == 0) && (i == 2));
446: if (DEBUGLEVEL > 3) msgtimer("lifting to prec %ld", E[i-1]);
447: }
448:
449: if (flag)
450: {
451: E = cgetg(4, t_VEC);
452: E[1] = lstoi(e0);
453: E[2] = (long)link;
454: E[3] = (long)v;
455: E[4] = (long)w;
456: }
457: else
458: {
459: E = cgetg(k+1, t_VEC);
460: for (i = 1; i <= 2*k-2; i++)
461: {
462: long t = link[i];
463: if (t < 0) E[-t] = v[i];
464: }
465: }
466: return E;
467: }
468:
469: /* Q list of (coprime, monic) factors of pol mod (T,p). Lift mod p^e = pe.
470: * T may be NULL */
471: GEN
472: hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pe, long e)
473: {
474: if (lg(Q) == 2) { GEN d = cgetg(2, t_VEC); d[1] = (long)pol; return d; }
475: pol = FpXQX_normalize(pol, T, pe);
476: return MultiLift(pol, Q, T, p, e, 0);
477: }
478:
479: /* U = NULL treated as 1 */
480: static void
481: BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
482: {
483: GEN Q, R;
484: if (j < 0) return;
485:
486: Q = FpX_mul((GEN)v[j], (GEN)w[j], pe);
487: if (U)
488: {
489: Q = FpXQ_mul(Q, U, f, pe);
490: R = FpX_sub(U, Q, pe);
491: }
492: else
493: R = FpX_Fp_add(FpX_neg(Q, pe), gun, pe);
494: w[j+1] = (long)Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
495: w[j ] = (long)R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
496: BezoutPropagate(link, v, w, pe, R, f, link[j ]);
497: BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
498: }
499:
500: /* as above, but return the Bezout coefficients for the lifted modular factors
501: * U[i] = 1 mod Qlift[i]
502: * 0 mod Qlift[j], j != i */
503: GEN
504: bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
505: {
506: GEN E, link, v, w, pe;
507: long i, k = lg(Q)-1;
508: if (k == 1) { GEN d = cgetg(2, t_VEC); d[1] = (long)pol; return d; }
509: pe = gpowgs(p, e);
510: pol = FpX_normalize(pol, pe);
511: E = MultiLift(pol, Q, NULL, p, e, 1);
512: link = (GEN) E[2];
513: v = (GEN) E[3];
514: w = (GEN) E[4];
515: BezoutPropagate(link, v, w, pe, NULL, pol, lg(v) - 2);
516: E = cgetg(k+1, t_VEC);
517: for (i = 1; i <= 2*k-2; i++)
518: {
519: long t = link[i];
520: if (t < 0) E[-t] = w[i];
521: }
522: return gcopy(E);
523: }
524:
525: /* Front-end for hensel_lift_fact:
526: lift the factorization of pol mod p given by fct to p^exp (if possible) */
527: GEN
528: polhensellift(GEN pol, GEN fct, GEN p, long exp)
529: {
530: GEN p1, p2;
531: long av = avma, i, j, l;
532:
533: /* we check the arguments */
534: if (typ(pol) != t_POL)
535: err(talker, "not a polynomial in polhensellift");
536: if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3))
537: err(talker, "not a factorization in polhensellift");
538: if (typ(p) != t_INT || !isprime(p))
539: err(talker, "not a prime number in polhensellift");
540: if (exp < 1)
541: err(talker, "not a positive exponent in polhensellift");
542:
543: p1 = lift(fct); /* make sure the coeffs are integers and not intmods */
544: l = lg(p1) - 1;
545: for (i=1; i<=l; i++)
546: {
547: p2 = (GEN)p1[i];
548: if (typ(p2) != t_POL)
549: {
550: if (typ(p2) != t_INT)
551: err(talker, "not an integral factorization in polhensellift");
552: p1[i] = (long)scalarpol(p2, varn(pol));
553: }
554: }
555:
556: /* then we check that pol \equiv \prod f ; f \in fct mod p */
557: p2 = (GEN)p1[1];
558: for (j = 2; j <= l; j++) p2 = FpX_mul(p2, (GEN)p1[j], p);
559: if (!gcmp0(FpX_sub(pol, p2, p)))
560: err(talker, "not a correct factorization in polhensellift");
561:
562: /* finally we check that the elements of fct are coprime mod p */
563: if (!FpX_is_squarefree(pol, p))
564: {
565: for (i = 1; i <= l; i++)
566: for (j = 1; j < i; j++)
567: if (lgef(FpX_gcd((GEN)p1[i], (GEN)p1[j], p)) != 3)
568: err(talker, "polhensellift: factors %Z and %Z are not coprime",
569: p1[i], p1[j]);
570: }
571: return gerepilecopy(av, hensel_lift_fact(pol,p1,NULL,p,gpowgs(p,exp),exp));
572: }
573:
574: #if 0
575: /* cf Beauzamy et al: upper bound for
576: * lc(x) * [2^(5/8) / pi^(3/8)] e^(1/4n) 2^(n/2) sqrt([x]_2)/ n^(3/8)
577: * where [x]_2 = sqrt(\sum_i=0^n x[i]^2 / binomial(n,i)). One factor has
578: * all coeffs less than then bound */
579: static GEN
580: two_factor_bound(GEN x)
581: {
582: long av = avma, i,j, n = lgef(x) - 3;
583: GEN *invbin, c, r = cgetr(3), z;
584:
585: x += 2; invbin = (GEN*)new_chunk(n+1);
586: z = realun(3); /* invbin[i] = 1 / binomial(n, i) */
587: for (i=0,j=n; j >= i; i++,j--)
588: {
589: invbin[i] = invbin[j] = z;
590: z = divrs(mulrs(z, i+1), n-i);
591: }
592: z = invbin[0]; /* = 1 */
593: for (i=0; i<=n; i++)
594: {
595: c = (GEN)x[i]; if (!signe(c)) continue;
596: affir(c, r);
597: z = addrr(z, mulrr(gsqr(r), invbin[i]));
598: }
599: z = shiftr(mpsqrt(z), n);
600: z = divrr(z, dbltor(pow((double)n, 0.75)));
601: z = grndtoi(mpsqrt(z), &i);
602: z = mulii(z, absi((GEN)x[n]));
603: return gerepileuptoint(av, shifti(z, 1));
604: }
605: #endif
606:
607: static GEN
608: uniform_Mignotte_bound(GEN x)
609: {
610: long e, n = degpol(x);
611: GEN p1, N2 = mpsqrt(QuickNormL2(x,DEFAULTPREC));
612:
613: p1 = grndtoi(gmul2n(N2, n), &e);
614: if (e>=0) p1 = addii(p1, shifti(gun, e));
615: return p1;
616: }
617:
618: /* all factors have coeffs less than the bound */
619: static GEN
620: all_factor_bound(GEN x)
621: {
622: long i, n = lgef(x) - 3;
623: GEN t, z = gzero;
624: for (i=2; i<=n+2; i++) z = addii(z,sqri((GEN)x[i]));
625: t = absi((GEN)x[n+2]);
626: z = addii(t, addsi(1, racine(z)));
627: z = mulii(z, binome(stoi(n-1), n>>1));
628: return shifti(mulii(t,z),1);
629: }
630:
631: /* x defined mod p^a, return mods ( (x - mods(x, p^b)) / p^b , p^(b-a) ) */
632: static GEN
633: TruncTrace(GEN x, GEN pb, GEN pa_b, GEN pa_bs2, GEN pbs2)
634: {
635: GEN r, q = dvmdii(x, pb, &r);
636: if (cmpii(r, pbs2) > 0) q = addis(q,1);
637: if (pa_bs2 && cmpii(q,pa_bs2) > 0) q = subii(q,pa_b);
638: return q;
639: }
640:
641: /* Naive recombination of modular factors: combine up to maxK modular
642: * factors, degree <= klim and divisible by hint
643: *
644: * target = polynomial we want to factor
645: * famod = array of modular factors. Product should be congruent to
646: * target/lc(target) modulo p^a
647: * All rational factors are bounded by p^b, b <= a and p^(b-a) < 2^(BIL-1) */
648: static GEN
649: cmbf(GEN target, GEN famod, GEN p, long b, long a,
650: long maxK, long klim,long hint)
651: {
652: long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
653: long spa_b, spa_bs2;
654: GEN lc, lctarget, pa = gpowgs(p,a), pas2 = shifti(pa,-1);
655: GEN trace = cgetg(lfamod+1, t_VECSMALL);
656: GEN ind = cgetg(lfamod+1, t_VECSMALL);
657: GEN degpol = cgetg(lfamod+1, t_VECSMALL);
658: GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
659: GEN listmod = cgetg(lfamod+1, t_COL);
660: GEN fa = cgetg(lfamod+1, t_COL);
661: GEN res = cgetg(3, t_VEC);
662: GEN bound = all_factor_bound(target);
663:
664: if (maxK < 0) maxK = lfamod-1;
665:
666: lc = absi(leading_term(target));
667: lctarget = gmul(lc,target);
668:
669: {
670: GEN pa_b,pa_bs2,pb,pbs2;
671: pa_b = gpowgs(p, a-b); /* make sure p^(a-b) < 2^(BIL-1) */
672: while (is_bigint(pa_b)) { b++; pa_b = divii(pa_b, p); }
673: pa_bs2 = shifti(pa_b,-1);
674: pb= gpowgs(p, b); pbs2 = shifti(pb,-1);
675: for (i=1; i <= lfamod; i++)
676: {
677: GEN T, p1 = (GEN)famod[i];
678: degpol[i] = degpol(p1);
679: if (!gcmp1(lc))
680: T = modii(mulii(lc, (GEN)p1[degpol[i]+1]), pa);
681: else
682: T = (GEN)p1[degpol[i]+1]; /* d-1 term */
683: trace[i] = itos( TruncTrace(T, pb,pa_b,NULL,pbs2) );
684: }
685: spa_b = pa_b[2]; /* < 2^(BIL-1) */
686: spa_bs2 = pa_bs2[2]; /* < 2^(BIL-1) */
687: }
688: degsofar[0] = 0; /* sentinel */
689:
690: /* ind runs through strictly increasing sequences of length K,
691: * 1 <= ind[i] <= lfamod */
692: nextK:
693: if (K > maxK || 2*K > lfamod) goto END;
694: if (DEBUGLEVEL > 4)
695: fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K));
696: setlg(ind, K+1); ind[1] = 1;
697: i = 1; curdeg = degpol[ind[1]];
698: for(;;)
699: { /* try all combinations of K factors */
700: for (j = i; j < K; j++)
701: {
702: degsofar[j] = curdeg;
703: ind[j+1] = ind[j]+1; curdeg += degpol[ind[j+1]];
704: }
705: if (curdeg <= klim && curdeg % hint == 0) /* trial divide */
706: {
707: GEN y, q, cont, list;
708: ulong av;
709: long t;
710:
711: /* d - 1 test, overflow is not a problem (correct mod 2^BIL) */
712: for (t=trace[ind[1]],i=2; i<=K; i++)
713: t = addssmod(t, trace[ind[i]], spa_b);
714: if (t > spa_bs2) t -= spa_b;
715: if (labs(t) > ((K+1)>>1))
716: {
717: if (DEBUGLEVEL>6) fprintferr(".");
718: goto NEXT;
719: }
720: av = avma;
721:
722: /* check trailing coeff */
723: y = lc;
724: for (i=1; i<=K; i++)
725: y = centermod_i(mulii(y, gmael(famod,ind[i],2)), pa, pas2);
726: if (!signe(y) || resii((GEN)lctarget[2], y) != gzero)
727: {
728: if (DEBUGLEVEL>6) fprintferr("T");
729: avma = av; goto NEXT;
730: }
731: y = lc; /* full computation */
732: for (i=1; i<=K; i++)
733: y = centermod_i(gmul(y, (GEN)famod[ind[i]]), pa, pas2);
734:
735: /* y is the candidate factor */
736: if (! (q = polidivis(lctarget,y,bound)) )
737: {
738: if (DEBUGLEVEL>6) fprintferr("*");
739: avma = av; goto NEXT;
740: }
741: /* found a factor */
742: cont = content(y);
743: if (signe(leading_term(y)) < 0) cont = negi(cont);
744: y = gdiv(y, cont);
745:
746: list = cgetg(K+1, t_VEC);
747: listmod[cnt] = (long)list;
748: for (i=1; i<=K; i++) list[i] = famod[ind[i]];
749:
750: fa[cnt++] = (long)y;
751: /* fix up target */
752: target = gdiv(q, leading_term(y));
753: for (i=j=k=1; i <= lfamod; i++)
754: { /* remove used factors */
755: if (j <= K && i == ind[j]) j++;
756: else
757: {
758: famod[k] = famod[i];
759: trace[k] = trace[i];
760: degpol[k] = degpol[i]; k++;
761: }
762: }
763: lfamod -= K;
764: if (lfamod < 2*K) goto END;
765: i = 1; curdeg = degpol[ind[1]];
766: bound = all_factor_bound(target);
767: lc = absi(leading_term(target));
768: lctarget = gmul(lc,target);
769: if (DEBUGLEVEL > 2)
770: {
771: fprintferr("\n"); msgtimer("to find factor %Z",y);
772: fprintferr("remaining modular factor(s): %ld\n", lfamod);
773: }
774: continue;
775: }
776:
777: NEXT:
778: for (i = K+1;;)
779: {
780: if (--i == 0) { K++; goto nextK; }
781: if (++ind[i] <= lfamod - K + i)
782: {
783: curdeg = degsofar[i-1] + degpol[ind[i]];
784: if (curdeg <= klim) break;
785: }
786: }
787: }
788: END:
789: if (degpol(target) > 0)
790: { /* leftover factor */
791: if (signe(leading_term(target)) < 0) target = gneg_i(target);
792:
793: setlg(famod, lfamod+1);
794: listmod[cnt] = (long)dummycopy(famod);
795: fa[cnt++] = (long)target;
796: }
797: if (DEBUGLEVEL>6) fprintferr("\n");
798: setlg(listmod, cnt); setlg(fa, cnt);
799: res[1] = (long)fa;
800: res[2] = (long)listmod; return res;
801: }
802:
803: void
804: factor_quad(GEN x, GEN res, long *ptcnt)
805: {
806: GEN a = (GEN)x[4], b = (GEN)x[3], c = (GEN)x[2], d, u, z1, z2, t;
807: GEN D = subii(sqri(b), shifti(mulii(a,c), 2));
808: long v, cnt = *ptcnt;
809:
810: if (!carrecomplet(D, &d)) { res[cnt++] = (long)x; *ptcnt = cnt; return; }
811:
812: t = shifti(negi(addii(b, d)), -1);
813: z1 = gdiv(t, a); u = denom(z1);
814: z2 = gdiv(addii(t, d), a);
815: v = varn(x);
816: res[cnt++] = lmul(u, gsub(polx[v], z1)); u = divii(a, u);
817: res[cnt++] = lmul(u, gsub(polx[v], z2)); *ptcnt = cnt;
818: }
819:
820: /* y > 1 and B integers. Let n such that y^(n-1) <= B < y^n.
821: * Return e = max(n,1), set *ptq = y^e if non-NULL */
822: long
823: logint(GEN B, GEN y, GEN *ptq)
824: {
825: ulong av = avma;
826: long e,i,fl;
827: GEN q,pow2, r = y;
828:
829: if (typ(B) != t_INT) B = ceil_safe(B);
830: if (expi(B) <= (expi(y) << 6)) /* e small, be naive */
831: {
832: for (e=1; cmpii(r, B) < 0; e++) r = mulii(r,y);
833: goto END;
834: }
835: /* binary splitting: compute bits of e one by one */
836: /* compute pow2[i] = y^(2^i) [i < very crude upper bound for log_2(n)] */
837: pow2 = new_chunk(bit_accuracy(lgefint(B)));
838: pow2[0] = (long)y;
839: for (i=0,q=r;; )
840: {
841: fl = cmpii(r,B); if (fl >= 0) break;
842: q = r; r = sqri(q);
843: i++; pow2[i] = (long)r;
844: }
845: if (i == 0) { e = 1; goto END; } /* y <= B */
846:
847: for (i--, e=1<<i;;)
848: { /* y^e = q < B <= r = q * y^(2^i) */
849: if (!fl) break; /* B = r */
850: /* q < B < r */
851: if (--i < 0) { if (fl > 0) e++; break; }
852: r = mulii(q, (GEN)pow2[i]);
853: fl = cmpii(r, B);
854: if (fl <= 0) { e += (1<<i); q = r; }
855: }
856: if (fl <= 0) { e++; r = mulii(r,y); }
857: END:
858: if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av;
859: return e;
860: }
861:
862: /* recombination of modular factors: van Hoeij's algorithm */
863:
864: /* compute Newton sums of P (i-th powers of roots, i=1..n)
865: * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
866: * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0) */
867: GEN
868: polsym_gen(GEN P, GEN y0, long n, GEN N)
869: {
870: long av1,av2,dP=degpol(P),i,k,m;
871: GEN s,y,P_lead;
872:
873: if (n<0) err(impl,"polsym of a negative n");
874: if (typ(P) != t_POL) err(typeer,"polsym");
875: if (!signe(P)) err(zeropoler,"polsym");
876: y = cgetg(n+2,t_COL);
877: if (y0)
878: {
879: if (typ(y0) != t_COL) err(typeer,"polsym_gen");
880: m = lg(y0)-1;
881: for (i=1; i<=m; i++) y[i] = lcopy((GEN)y0[i]);
882: }
883: else
884: {
885: m = 1;
886: y[1] = lstoi(dP);
887: }
888: P += 2; /* strip codewords */
889:
890: P_lead=(GEN)P[dP]; if (gcmp1(P_lead)) P_lead=NULL;
891: if (N && P_lead)
892: if (!invmod(P_lead,N,&P_lead))
893: err(talker,"polsyn: non-invertible leading coeff: %Z", P_lead);
894: for (k=m; k<=n; k++)
895: {
896: av1=avma; s = (dP>=k)? gmulsg(k,(GEN)P[dP-k]): gzero;
897: for (i=1; i<k && i<=dP; i++)
898: s = gadd(s, gmul((GEN)y[k-i+1],(GEN)P[dP-i]));
899: if (N)
900: {
901: s = modii(s, N);
902: if (P_lead) { s = mulii(s,P_lead); s = modii(s,N); }
903: }
904: else
905: if (P_lead) s = gdiv(s,P_lead);
906: av2=avma; y[k+1]=lpile(av1,av2,gneg(s));
907: }
908: return y;
909: }
910:
911: /* return integer y such that all roots of P are less than y */
912: static GEN
913: root_bound(GEN P)
914: {
915: GEN P0 = dummycopy(P), lP = absi(leading_term(P)), x,y,z;
916: long k,d = degpol(P);
917:
918: setlgef(P0, d+2); /* P = lP x^d + P0 */
919: P = P0+2; /* strip codewords */
920: for (k=0; k<d; k++) P[k] = labsi((GEN)P[k]);
921:
922: x = y = gun;
923: for (k=0; ; k++)
924: {
925: if (cmpii(poleval(P0,y), mulii(lP, gpowgs(y, d))) < 0) break;
926: x = y; y = shifti(y,1);
927: }
928: for(;;)
929: {
930: z = shifti(addii(x,y), -1);
931: if (egalii(x,z)) break;
932: if (cmpii(poleval(P0,z), mulii(lP, gpowgs(z, d))) < 0)
933: y = z;
934: else
935: x = z;
936: }
937: return y;
938: }
939:
940: extern GEN gscal(GEN x,GEN y);
941: extern GEN sqscal(GEN x);
942:
943: /* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the
944: * vector of the (f_i . f_i)*/
945: GEN
946: gram_schmidt(GEN e, GEN *ptB)
947: {
948: long i,j,lx = lg(e);
949: GEN f = dummycopy(e), B, iB;
950:
951: B = cgetg(lx, t_VEC);
952: iB= cgetg(lx, t_VEC);
953:
954: for (i=1;i<lx;i++)
955: {
956: GEN p1 = NULL;
957: ulong av;
958: B[i] = (long)sqscal((GEN)f[i]);
959: iB[i]= linv((GEN)B[i]); av = avma;
960: for (j=1; j<i; j++)
961: {
962: GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]);
963: GEN p2 = gmul(mu, (GEN)f[j]);
964: p1 = p1? gadd(p1,p2): p2;
965: }
966: p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i];
967: f[i] = (long)p1;
968: }
969: *ptB = B; return f;
970: }
971:
972: extern GEN hnfall_i(GEN A, GEN *ptB, long remove);
973:
974: GEN
975: special_pivot(GEN x)
976: {
977: GEN t, H = hnfall_i(x, NULL, 1);
978: long i,j, l = lg(H), h = lg(H[1]);
979: for (i=1; i<h; i++)
980: {
981: int fl = 0;
982: for (j=1; j<l; j++)
983: {
984: t = gcoeff(H,i,j);
985: if (signe(t))
986: {
987: if (!is_pm1(t) || fl) return NULL;
988: fl = 1;
989: }
990: }
991: }
992: return H;
993: }
994:
995: /* x matrix: compute a bound for | \sum e_i x_i | ^ 2, e_i = 0,1 */
996: static GEN
997: my_norml2(GEN x)
998: {
999: long i,j, co = lg(x), li;
1000: GEN S = gzero;
1001: if (co == 1) return S;
1002: li = lg(x[1]);
1003: for (i=1; i<li; i++)
1004: {
1005: GEN p = gzero;
1006: GEN m = gzero;
1007: for (j=1; j<co; j++)
1008: {
1009: GEN u = gcoeff(x,i,j);
1010: long s = gsigne(u);
1011: if (s < 0) m = gadd(m,u);
1012: else if (s > 0) p = gadd(p,u);
1013: }
1014: if (m != gzero) m = gneg(m);
1015: if (gcmp(p,m) > 0) m = p;
1016:
1017: S = gadd(S, gsqr(m));
1018: }
1019: return S;
1020: }
1021:
1022: /* each entry < 2^e */
1023: static long
1024: init_padic_prec(long e, int BitPerFactor, long n0, double LOGp2)
1025: {
1026: /* long b, goodb = (long) ((e - 0.175 * n0 * n0) * LOGp2); */
1027: long b, goodb = (long) (((e - BitPerFactor*n0)) * LOGp2);
1028: b = (long) ((e - 32)*LOGp2); if (b < goodb) goodb = b;
1029: /* b = (long) ((e - Max*32)*LOGp2); if (b > goodb) goodb = b; */
1030: return goodb;
1031: }
1032:
1033: extern GEN sindexrank(GEN x);
1034: extern GEN vconcat(GEN Q1, GEN Q2);
1035: extern GEN gauss_intern(GEN a, GEN b);
1036:
1037: /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
1038: * van Hoeij's knapsack
1039: *
1040: * P = monic squarefree in Z[X].
1041: * famod = array of (lifted) modular factors mod p^a
1042: * bound = Mignotte bound for the size of divisors of P (sor the sup norm)
1043: * previously recombined all set of factors with less than rec elts
1044: */
1045: GEN
1046: LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
1047: {
1048: long s = 2; /* # of traces added at each step */
1049: long BitPerFactor = 3; /* nb bits in p^(a-b) / modular factor */
1050: long i,j,C,r,S,tmax,n0,n,dP = degpol(P);
1051: double logp = log(gtodouble(p)), LOGp2 = LOG2/logp;
1052: double b0 = log((double)dP*2) / logp;
1053: double k = gtodouble(glog(root_bound(P), DEFAULTPREC)) / logp;
1054: GEN y, T, T2, TT, BL, m, u, norm, target, M, piv, list;
1055: GEN run = realun(DEFAULTPREC);
1056: ulong av,av2, lim;
1057: int did_recompute_famod = 0;
1058:
1059: n0 = n = r = lg(famod) - 1;
1060: list = cgetg(n0+1, t_COL);
1061:
1062: av = avma; lim = stack_lim(av, 1);
1063: TT = cgetg(n0+1, t_VEC);
1064: T = cgetg(n0+1, t_MAT);
1065: for (i=1; i<=n0; i++)
1066: {
1067: TT[i] = 0;
1068: T [i] = lgetg(s+1, t_COL);
1069: }
1070: BL = idmat(n0);
1071: /* tmax = current number of traces used (and computed so far)
1072: * S = number of traces used at the round's end = tmax + s */
1073: for(tmax = 0;; tmax = S)
1074: {
1075: GEN pas2, pa_b, BE;
1076: long b, goodb;
1077: double Nx;
1078:
1079: if (DEBUGLEVEL>3)
1080: fprintferr("LLL_cmbf: %ld potential factors (tmax = %ld)\n", r, tmax);
1081: av2 = avma;
1082: if (tmax > 0)
1083: { /* bound small vector in terms of a modified L2 norm of a
1084: * left inverse of BL */
1085: GEN z = gauss_intern(BL,NULL); /* 1/BL */
1086: if (!z) /* not maximal rank */
1087: {
1088: avma = av2;
1089: BL = hnfall_i(BL,NULL,1);
1090: r = lg(BL)-1; z = invmat(BL);
1091: av2 = avma;
1092: }
1093: Nx = gtodouble(my_norml2(gmul(run, z)));
1094: avma = av2;
1095: }
1096: else
1097: Nx = (double)n0; /* first time: BL = id */
1098: C = (long)sqrt(s*n0*n0/4. / Nx);
1099: if (C == 0) C = 1; /* frequent after a few iterations */
1100: M = dbltor((Nx * C*C + s*n0*n0/4.) * 1.00001);
1101:
1102: S = tmax + s;
1103: b = (long)ceil(b0 + S*k);
1104: if (a <= b)
1105: {
1106: a = (long)ceil(b + 3*s*k) + 1; /* enough for 3 more rounds */
1107: pa = gpowgs(p,a);
1108: famod = hensel_lift_fact(P,famod,NULL,p,pa,a);
1109: /* recompute old Newton sums to new precision */
1110: for (i=1; i<=n0; i++)
1111: TT[i] = (long)polsym_gen((GEN)famod[i], NULL, tmax, pa);
1112: did_recompute_famod = 1;
1113: }
1114: for (i=1; i<=n0; i++)
1115: {
1116: GEN p1 = (GEN)T[i];
1117: GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], S, pa);
1118: TT[i] = (long)p2;
1119: p2 += 1+tmax; /* ignore traces number 0...tmax */
1120: for (j=1; j<=s; j++) p1[j] = p2[j];
1121: }
1122:
1123: av2 = avma;
1124: T2 = gmod( gmul(T, BL), pa );
1125: goodb = init_padic_prec(gexpo(T2), BitPerFactor, n0, LOGp2);
1126: if (goodb > b) b = goodb;
1127: if (DEBUGLEVEL>2)
1128: fprintferr("LLL_cmbf: (a, b) = (%ld, %ld), expo(T) = %ld\n",
1129: a,b,gexpo(T2));
1130: pa_b = gpowgs(p, a-b);
1131: {
1132: GEN pb = gpowgs(p, b);
1133: GEN pa_bs2 = shifti(pa_b,-1);
1134: GEN pbs2 = shifti(pb,-1);
1135: for (i=1; i<=r; i++)
1136: {
1137: GEN p1 = (GEN)T2[i];
1138: for (j=1; j<=s; j++)
1139: p1[j] = (long)TruncTrace((GEN)p1[j],pb,pa_b,pa_bs2,pbs2);
1140: }
1141: }
1142: if (gcmp0(T2)) { avma = av2; continue; }
1143:
1144: BE = cgetg(s+1, t_MAT);
1145: for (i=1; i<=s; i++)
1146: {
1147: BE[i] = (long)zerocol(r+s);
1148: coeff(BE, i+r, i) = (long)pa_b;
1149: }
1150: m = concatsp( vconcat( gscalmat(stoi(C), r), T2 ), BE );
1151: /* [ C 0 ]
1152: * m = [ ] square matrix
1153: * [ T2 p^(a-b) I_S ] T2 = T * BL truncated
1154: */
1155: u = lllgramintern(gram_matrix(m), 4, 0, 0);
1156: m = gmul(m,u);
1157: (void)gram_schmidt(gmul(run,m), &norm);
1158: for (i=r+s; i>0; i--)
1159: if (cmprr((GEN)norm[i], M) < 0) break;
1160: if (i > r)
1161: { /* no progress */
1162: avma = av2; BitPerFactor += 2;
1163: if (DEBUGLEVEL>2)
1164: fprintferr("LLL_cmbf: increasing BitPerFactor = %ld\n", BitPerFactor);
1165: #if 0
1166: s++; for (i=1; i<=n; i++) T[i] = lgetg(s+1, t_COL);
1167: if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: increasing s = %ld\n", s);
1168: #endif
1169: continue;
1170: }
1171:
1172: n = r; r = i;
1173: if (r <= 1)
1174: {
1175: if (r == 0) err(bugparier,"LLL_cmbf [no factor]");
1176: if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: 1 factor\n");
1177: list[1] = (long)P; setlg(list,2); return list;
1178: }
1179:
1180: setlg(u, r+1);
1181: for (i=1; i<=r; i++) setlg(u[i], n+1);
1182: BL = gerepileupto(av2, gmul(BL, u));
1183: if (low_stack(lim, stack_lim(av,1)))
1184: {
1185: GEN *gptr[4]; gptr[0]=&BL; gptr[1]=&TT; gptr[2]=&T; gptr[3]=&famod;
1186: if(DEBUGMEM>1) err(warnmem,"LLL_cmbf");
1187: gerepilemany(av, gptr, did_recompute_famod? 4: 3);
1188: }
1189: if (r*rec >= n0) continue;
1190:
1191: av2 = avma;
1192: piv = special_pivot(BL);
1193: if (DEBUGLEVEL) fprintferr("special_pivot output:\n%Z\n",piv);
1194: if (!piv) { avma = av2; continue; }
1195:
1196: pas2 = shifti(pa,-1); target = P;
1197: for (i=1; i<=r; i++)
1198: {
1199: GEN p1 = (GEN)piv[i];
1200: if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);
1201:
1202: y = gun;
1203: for (j=1; j<=n0; j++)
1204: if (signe(p1[j]))
1205: y = centermod_i(gmul(y, (GEN)famod[j]), pa, pas2);
1206:
1207: /* y is the candidate factor */
1208: if (! (target = polidivis(target,y,bound)) ) break;
1209: list[i] = (long)y;
1210: }
1211: if (i > r)
1212: {
1213: if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: %ld factors\n", r);
1214: setlg(list,r+1); return list;
1215: }
1216: }
1217: }
1218:
1219: extern GEN primitive_pol_to_monic(GEN pol, GEN *ptlead);
1220:
1221: /* P(hx), in place. Assume P in Z[X], h in Z */
1222: void
1223: rescale_pol_i(GEN P, GEN h)
1224: {
1225: GEN hi = gun;
1226: long i, l = lgef(P);
1227: for (i=3; i<l; i++)
1228: {
1229: hi = mulii(hi,h);
1230: P[i] = lmulii((GEN)P[i], hi);
1231: }
1232: }
1233:
1234: /* P(hx) in Fp[X], in place. Assume P in Z[X], h in Z */
1235: void
1236: FpX_rescale_i(GEN P, GEN h, GEN p)
1237: {
1238: GEN hi = gun;
1239: long i, l = lgef(P);
1240: for (i=3; i<l; i++)
1241: {
1242: hi = modii(mulii(hi,h), p);
1243: P[i] = lmodii(mulii((GEN)P[i], hi), p);
1244: }
1245: }
1246:
1247: /* use van Hoeij's knapsack algorithm */
1248: static GEN
1249: combine_factors(GEN a, GEN famod, GEN p, long klim, long hint)
1250: {
1251: GEN B = uniform_Mignotte_bound(a), res,y,lt,L,pe,pE,listmod,p1;
1252: long i,E,e,l, maxK = 3, nft = lg(famod)-1;
1253:
1254: e = logint(B, p, &pe);
1255: if (DEBUGLEVEL > 4) fprintferr("Mignotte bound: %Z^%ld\n", p,e);
1256:
1257: { /* overlift for the d-1 test */
1258: int over = (int) (32 * LOG2 / log((double)itos(p)));
1259: pE = mulii(pe, gpowgs(p,over)); E = e + over;
1260: }
1261: famod = hensel_lift_fact(a,famod,NULL,p,pE,E);
1262:
1263: if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */
1264: else
1265: {
1266: lt = leading_term(a);
1267: if (!is_pm1(lt) && nft < 13) maxK = -1;
1268: }
1269: L = cmbf(a, famod, p, e, E, maxK, klim, hint);
1270:
1271: res = (GEN)L[1];
1272: listmod = (GEN)L[2]; l = lg(listmod);
1273: famod = (GEN)listmod[l-1];
1274: if (maxK >= 0 && lg(famod)-1 > 2*maxK)
1275: {
1276: a = (GEN)res[l-1];
1277: lt = leading_term(a);
1278: if (signe(lt) < 0) { a = gneg_i(a); lt = leading_term(a); }
1279: if (DEBUGLEVEL > 4) fprintferr("last factor still to be checked\n");
1280: if (gcmp1(lt))
1281: lt = NULL;
1282: else
1283: {
1284: GEN invlt, invLT;
1285: if (DEBUGLEVEL > 4) fprintferr("making it monic\n");
1286: a = primitive_pol_to_monic(a, <);
1287: B = uniform_Mignotte_bound(a);
1288: e = logint(B, p, &pe);
1289:
1290: invlt = mpinvmod(lt,p);
1291: for (i = 1; i<lg(famod); i++)
1292: { /* renormalize modular factors */
1293: p1 = (GEN)famod[i]; FpX_rescale_i(p1, invlt, p);
1294: invLT = powmodulo(lt, stoi(degpol(p1)), p);
1295: famod[i] = (long)FpX_Fp_mul(p1, invLT, p);
1296: }
1297: famod = hensel_lift_fact(a,famod,NULL,p,pe,e);
1298: }
1299: setlg(res, l-1); /* remove last elt (possibly unfactored) */
1300: L = LLL_cmbf(a, famod, p, pe, B, e, maxK);
1301: if (lt)
1302: {
1303: for (i=1; i<lg(L); i++)
1304: {
1305: y = (GEN)L[i]; rescale_pol_i(y, lt);
1306: L[i] = (long)primitive_part(y, &p1);
1307: }
1308: }
1309: res = concatsp(res, L);
1310: }
1311: return res;
1312: }
1313:
1314: extern long split_berlekamp(GEN *t, GEN pp, GEN pps2);
1315: #define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),(0),NULL)
1316:
1317: /* assume degree(a) > 0, a(0) != 0, and a squarefree */
1318: static GEN
1319: squff(GEN a, long hint)
1320: {
1321: GEN PolX,lead,res,prime,primes2,famod,y,g,z,w,*tabd,*tabdnew;
1322: long da = degpol(a), va = varn(a);
1323: long klim,chosenp,p,nfacp,lbit,j,d,e,np,nmax,nf,nft;
1324: ulong *tabbit, *tabkbit, *tmp, av = avma;
1325: byteptr pt=diffptr;
1326: const int MAXNP = max(5, (long)sqrt(da));
1327:
1328: if (hint <= 0) hint = 1;
1329: if (DEBUGLEVEL > 2) timer2();
1330: lbit = (da>>4)+1; nmax = da+1; klim = da>>1;
1331: chosenp = 0;
1332: tabd = NULL;
1333: tabdnew = (GEN*)new_chunk(nmax);
1334: tabbit = (ulong*)new_chunk(lbit);
1335: tabkbit = (ulong*)new_chunk(lbit);
1336: tmp = (ulong*)new_chunk(lbit);
1337: prime = icopy(gun);
1338: lead = (GEN)a[da+2]; PolX = u_Fp_FpX(polx[0],0, 2);
1339: for (p = np = 0; np < MAXNP; )
1340: {
1341: ulong av0 = avma;
1342: p += *pt++; if (!*pt) err(primer1);
1343: if (!smodis(lead,p)) continue;
1344: z = u_Fp_FpX(a,0, p);
1345: if (!u_FpX_is_squarefree(z, p)) { avma = av0; continue ; }
1346:
1347: for (j=0; j<lbit-1; j++) tabkbit[j] = 0;
1348: tabkbit[j] = 1;
1349: d = 0; e = da; nfacp = 0;
1350: w = PolX; prime[2] = p;
1351: while (d < (e>>1))
1352: {
1353: long lgg;
1354: d++; w = u_FpXQ_pow(w, prime, z, p);
1355: g = u_FpX_gcd(z, u_FpX_sub(w, PolX, p), p);
1356: lgg = degpol(g);
1357: if (lgg == 0) g = NULL;
1358: else
1359: {
1360: z = u_FpX_div(z, g, p); e = degpol(z);
1361: w = u_FpX_rem(w, z, p);
1362: lgg /= d; nfacp += lgg;
1363: if (DEBUGLEVEL>5)
1364: fprintferr(" %3ld factor%s of degree %3ld\n", lgg, lgg==1?"":"s",d);
1365: record_factors(lgg, d, lbit-1, tabkbit, tmp);
1366: }
1367: tabdnew[d] = g;
1368: }
1369: if (e > 0)
1370: {
1371: if (DEBUGLEVEL>5) fprintferr(" %3ld factor of degree %3ld\n",1,e);
1372: tabdnew[e] = z; nfacp++;
1373: record_factors(1, e, lbit-1, tabkbit, tmp);
1374: }
1375:
1376: if (np)
1377: for (j=0; j<lbit; j++) tabbit[j] &= tabkbit[j];
1378: else
1379: for (j=0; j<lbit; j++) tabbit[j] = tabkbit[j];
1380: if (DEBUGLEVEL > 4)
1381: fprintferr("...tried prime %3ld (%-3ld factor%s). Time = %ld\n",
1382: p, nfacp, nfacp==1?"": "s", timer2());
1383: if (min_deg(lbit-1,tabbit) > klim) { avma = av; return _col(a); }
1384: if (nfacp < nmax)
1385: {
1386: nmax = nfacp; tabd = tabdnew; chosenp = p;
1387: for (j=d+1; j<e; j++) tabd[j] = NULL;
1388: tabdnew = (GEN*)new_chunk(da);
1389: }
1390: else avma = av0;
1391: np++;
1392: }
1393: prime[2] = chosenp; primes2 = shifti(prime, -1);
1394: nf = nmax; nft = 1;
1395: y = cgetg(nf+1,t_COL); famod = cgetg(nf+1,t_COL);
1396: for (d = 1; nft <= nf; d++)
1397: {
1398: g = tabd[d];
1399: if (g)
1400: {
1401: long n = degpol(g)/d;
1402: famod[nft] = (long)small_to_pol(u_FpX_normalize(g, chosenp),va);
1403: if (n > 1 && n != split_berlekamp((GEN*)(famod+nft),prime,primes2))
1404: err(bugparier,"squff: wrong numbers of factors");
1405: nft += n;
1406: }
1407: }
1408: if (DEBUGLEVEL > 4) msgtimer("splitting mod p = %ld",chosenp);
1409: res = combine_factors(a, famod, prime, da-1, hint);
1410: return gerepilecopy(av, res);
1411: }
1412:
1413: /* A(X^d) --> A(X) */
1414: GEN
1415: poldeflate_i(GEN x0, long d)
1416: {
1417: GEN z, y, x;
1418: long i,id, dy, dx = degpol(x0);
1419: if (d == 1) return x0;
1420: if (dx < 0) return zeropol(varn(x0));
1421: dy = dx/d;
1422: y = cgetg(dy+3, t_POL);
1423: y[1] = evalsigne(1)|evaldeg(dy)|evalvarn(varn(x0));
1424: z = y + 2;
1425: x = x0+ 2;
1426: for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1427: return y;
1428: }
1429:
1430: long
1431: checkdeflate(GEN x)
1432: {
1433: long d = 0, i, lx = lgef(x);
1434: for (i=3; i<lx; i++)
1435: if (!gcmp0((GEN)x[i])) { d = cgcd(d,i-2); if (d == 1) break; }
1436: return d;
1437: }
1438:
1439: GEN
1440: gdeflate(GEN x, long v, long d)
1441: {
1442: long i, lx, tx = typ(x);
1443: GEN z;
1444: if (is_scalar_t(tx)) return gcopy(x);
1445: if (tx == t_POL)
1446: {
1447: long vx = varn(x);
1448: ulong av;
1449: if (vx < v)
1450: {
1451: lx = lgef(x);
1452: z = cgetg(lx,t_POL); z[1] = x[1];
1453: for (i=2; i<lx; i++) z[i] = (long)gdeflate((GEN)x[i],v,d);
1454: return z;
1455: }
1456: if (vx > v) return gcopy(x);
1457: av = avma;
1458: if (checkdeflate(x) % d != 0)
1459: err(talker,"impossible substitution in gdeflate");
1460: return gerepilecopy(av, poldeflate_i(x,d));
1461: }
1462: if (tx == t_RFRAC)
1463: {
1464: z = cgetg(3, t_RFRAC);
1465: z[1] = (long)gdeflate((GEN)x[1],v,d);
1466: z[2] = (long)gdeflate((GEN)x[2],v,d);
1467: return z;
1468: }
1469: if (is_matvec_t(tx))
1470: {
1471: lx = lg(x); z = cgetg(lx, tx);
1472: for (i=1; i<lx; i++) z[i] = (long)gdeflate((GEN)x[i],v,d);
1473: return z;
1474: }
1475: err(typeer,"gdeflate");
1476: return NULL; /* not reached */
1477: }
1478:
1479: /* set *m to the largest d such that x0 = A(X^d); return A */
1480: GEN
1481: poldeflate(GEN x, long *m)
1482: {
1483: *m = checkdeflate(x);
1484: return poldeflate_i(x, *m);
1485: }
1486:
1487: /* return x0(X^d) */
1488: GEN
1489: polinflate(GEN x0, long d)
1490: {
1491: long i, id, dy, dx = degpol(x0);
1492: GEN x = x0 + 2, z, y;
1493: dy = dx*d;
1494: y = cgetg(dy+3, t_POL);
1495: y[1] = evalsigne(1)|evaldeg(dy)|evalvarn(varn(x0));
1496: z = y + 2;
1497: for (i=0; i<=dy; i++) z[i] = zero;
1498: for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1499: return y;
1500: }
1501:
1502: GEN
1503: squff2(GEN x, long hint)
1504: {
1505: GEN L;
1506: long m;
1507: x = poldeflate(x, &m);
1508: L = squff(x, hint);
1509: if (m > 1)
1510: {
1511: GEN e, v, fa = decomp(stoi(m));
1512: long i,j,k, l;
1513:
1514: e = (GEN)fa[2]; k = 0;
1515: fa= (GEN)fa[1]; l = lg(fa);
1516: for (i=1; i<l; i++)
1517: {
1518: e[i] = itos((GEN)e[i]);
1519: k += e[i];
1520: }
1521: v = cgetg(k+1, t_VECSMALL); k = 1;
1522: for (i=1; i<l; i++)
1523: for (j=1; j<=e[i]; j++) v[k++] = itos((GEN)fa[i]);
1524: for (k--; k; k--)
1525: {
1526: GEN L2 = cgetg(1,t_VEC);
1527: for (i=1; i < lg(L); i++)
1528: L2 = concatsp(L2, squff(polinflate((GEN)L[i], v[k]), hint));
1529: L = L2;
1530: }
1531: }
1532: return L;
1533: }
1534:
1535: /* Factor x in Z[t]. Assume all factors have degree divisible by hint */
1536: GEN
1537: factpol(GEN x, long hint)
1538: {
1539: GEN fa,p1,d,t,v,w, y = cgetg(3,t_MAT);
1540: long av=avma,av2,lx,vx,i,j,k,ex,nbfac,zval;
1541:
1542: if (typ(x)!=t_POL) err(notpoler,"factpol");
1543: if (!signe(x)) err(zeropoler,"factpol");
1544:
1545: ex = nbfac = 0;
1546: p1 = x+2; while (gcmp0((GEN)*p1)) p1++;
1547: zval = p1 - (x + 2);
1548: lx = lgef(x) - zval;
1549: vx = varn(x);
1550: if (zval)
1551: {
1552: x = cgetg(lx, t_POL); p1 -= 2;
1553: for (i=2; i<lx; i++) x[i] = p1[i];
1554: x[1] = evalsigne(1)|evalvarn(vx)|evallgef(lx);
1555: nbfac++;
1556: }
1557: /* now x(0) != 0 */
1558: if (lx==3) { fa = NULL;/* for lint */ goto END; }
1559: p1 = cgetg(1,t_VEC); fa=cgetg(lx,t_VEC);
1560: for (i=1; i<lx; i++) fa[i] = (long)p1;
1561: d=content(x); if (gsigne(leading_term(x)) < 0) d = gneg_i(d);
1562: if (!gcmp1(d)) x=gdiv(x,d);
1563: if (lx==4) { nbfac++; ex++; fa[1] = (long)concatsp(p1,x); goto END; }
1564:
1565: w=derivpol(x); t=modulargcd(x,w);
1566: if (!gcmp1(t)) { x=gdeuc(x,t); w=gdeuc(w,t); }
1567: k=1;
1568: while (k)
1569: {
1570: ex++; w=gadd(w, gneg_i(derivpol(x))); k=signe(w);
1571: if (k) { t=modulargcd(x,w); x=gdeuc(x,t); w=gdeuc(w,t); } else t=x;
1572: if (degpol(t) > 0)
1573: {
1574: fa[ex] = (long)squff2(t,hint);
1575: nbfac += lg(fa[ex])-1;
1576: }
1577: }
1578: END: av2=avma;
1579: v=cgetg(nbfac+1,t_COL); y[1]=(long)v;
1580: w=cgetg(nbfac+1,t_COL); y[2]=(long)w;
1581: if (zval) { v[1]=lpolx[vx]; w[1]=lstoi(zval); k=1; } else k=0;
1582: for (i=1; i<=ex; i++)
1583: for (j=1; j<lg(fa[i]); j++)
1584: {
1585: k++; v[k]=lcopy(gmael(fa,i,j)); w[k]=lstoi(i);
1586: }
1587: gerepilemanyvec(av,av2,y+1,2);
1588: return sort_factor(y, cmpii);
1589: }
1590:
1591: /***********************************************************************/
1592: /** **/
1593: /** FACTORIZATION **/
1594: /** **/
1595: /***********************************************************************/
1596: #define LT 17
1597: #define assign_or_fail(x,y) {\
1598: if (y==NULL) y=x; else if (!gegal(x,y)) return 0;\
1599: }
1600: #define tsh 6
1601: #define typs(x,y) ((x << tsh) | y)
1602: #define typ1(x) (x >> tsh)
1603: #define typ2(x) (x & ((1<<tsh)-1))
1604:
1605: static long
1606: poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)
1607: {
1608: long t[LT]; /* code for 0,1,2,3,61,62,63,67,7,81,82,83,86,87,91,93,97 */
1609: long tx = typ(x),lx,i,j,s,pa=BIGINT;
1610: GEN pcx=NULL, p=NULL,pol=NULL,p1,p2;
1611:
1612: if (is_scalar_t(tx))
1613: {
1614: if (tx == t_POLMOD) return 0;
1615: x = scalarpol(x,0);
1616: }
1617: for (i=2; i<LT; i++) t[i]=0; /* t[0..1] unused */
1618: lx = lgef(x);
1619: for (i=2; i<lx; i++)
1620: {
1621: p1=(GEN)x[i];
1622: switch(typ(p1))
1623: {
1624: case t_INT: case t_FRAC: case t_FRACN:
1625: break;
1626: case t_REAL:
1627: s = precision(p1); if (s < pa) pa = s;
1628: t[2]=1; break;
1629: case t_INTMOD:
1630: assign_or_fail((GEN)p1[1],p);
1631: t[3]=1; break;
1632: case t_COMPLEX:
1633: if (!pcx)
1634: {
1635: pcx = cgetg(5,t_POL); /* x^2 + 1 */
1636: pcx[1] = evalsigne(1)|evalvarn(0)|m_evallgef(5),
1637: pcx[2]=pcx[4]=un; pcx[3]=zero;
1638: }
1639: for (j=1; j<=2; j++)
1640: {
1641: p2 = (GEN)p1[j];
1642: switch(typ(p2))
1643: {
1644: case t_INT: case t_FRAC: case t_FRACN:
1645: assign_or_fail(pcx,pol);
1646: t[4]=1; break;
1647: case t_REAL:
1648: s = precision(p2); if (s < pa) pa = s;
1649: t[5]=1; break;
1650: case t_INTMOD:
1651: assign_or_fail((GEN)p2[1],p);
1652: assign_or_fail(pcx,pol);
1653: t[6]=1; break;
1654: case t_PADIC:
1655: s = precp(p2) + valp(p2); if (s < pa) pa = s;
1656: assign_or_fail((GEN)p2[2],p);
1657: assign_or_fail(pcx,pol);
1658: t[7]=1; break;
1659: default: return 0;
1660: }
1661: }
1662: break;
1663: case t_PADIC:
1664: s = precp(p1) + valp(p1); if (s < pa) pa = s;
1665: assign_or_fail((GEN)p1[2],p);
1666: t[8]=1; break;
1667: case t_QUAD:
1668: for (j=2; j<=3; j++)
1669: {
1670: p2 = (GEN)p1[j];
1671: switch(typ(p2))
1672: {
1673: case t_INT: case t_FRAC: case t_FRACN:
1674: assign_or_fail((GEN)p1[1],pol);
1675: t[9]=1; break;
1676: case t_REAL:
1677: s = precision(p2); if (s < pa) pa = s;
1678: if (gsigne(discsr((GEN)p1[1]))>0) t[10]=1; else t[12]=1;
1679: break;
1680: case t_INTMOD:
1681: assign_or_fail((GEN)p2[1],p);
1682: assign_or_fail((GEN)p1[1],pol);
1683: t[11]=1; break;
1684: case t_PADIC:
1685: s = precp(p2) + valp(p2); if (s < pa) pa = s;
1686: assign_or_fail((GEN)p2[2],p);
1687: assign_or_fail((GEN)p1[1],pol);
1688: t[13]=1; break;
1689: default: return 0;
1690: }
1691: }
1692: break;
1693: case t_POLMOD:
1694: assign_or_fail((GEN)p1[1],pol);
1695: for (j=1; j<=2; j++)
1696: {
1697: GEN pbis = NULL, polbis = NULL;
1698: long pabis;
1699: switch(poltype((GEN)p1[j],&pbis,&polbis,&pabis))
1700: {
1701: case t_INT: t[14]=1; break;
1702: case t_INTMOD: t[15]=1; break;
1703: case t_PADIC: t[16]=1; if (pabis<pa) pa=pabis; break;
1704: default: return 0;
1705: }
1706: if (pbis) assign_or_fail(pbis,p);
1707: if (polbis) assign_or_fail(polbis,pol);
1708: }
1709: break;
1710: default: return 0;
1711: }
1712: }
1713: if (t[5]||t[12])
1714: {
1715: if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
1716: *ptpa=pa; return t_COMPLEX;
1717: }
1718: if (t[2]||t[10])
1719: {
1720: if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
1721: *ptpa=pa; return t[4]?t_COMPLEX:t_REAL;
1722: }
1723: if (t[6]||t[11]||t[15])
1724: {
1725: *ptpol=pol; *ptp=p;
1726: i = t[15]? t_POLMOD: (t[11]? t_QUAD: t_COMPLEX);
1727: return typs(i, t_INTMOD);
1728: }
1729: if (t[7]||t[13]||t[16])
1730: {
1731: *ptpol=pol; *ptp=p; *ptpa=pa;
1732: i = t[16]? t_POLMOD: (t[13]? t_QUAD: t_COMPLEX);
1733: return typs(i, t_PADIC);
1734: }
1735: if (t[4]||t[9]||t[14])
1736: {
1737: *ptpol=pol;
1738: i = t[14]? t_POLMOD: (t[9]? t_QUAD: t_COMPLEX);
1739: return typs(i, t_INT);
1740: }
1741: if (t[3]) { *ptp=p; return t_INTMOD; }
1742: if (t[8]) { *ptp=p; *ptpa=pa; return t_PADIC; }
1743: return t_INT;
1744: }
1745: #undef LT
1746:
1747: GEN
1748: factor0(GEN x,long flag)
1749: {
1750: long tx=typ(x);
1751:
1752: if (flag<0) return factor(x);
1753: if (is_matvec_t(tx)) return gboundfact(x,flag);
1754: if (tx==t_INT || tx==t_FRAC || tx==t_FRACN) return boundfact(x,flag);
1755: err(talker,"partial factorization is not meaningful here");
1756: return NULL; /* not reached */
1757: }
1758:
1759: /* assume f and g coprime integer factorizations */
1760: GEN
1761: merge_factor_i(GEN f, GEN g)
1762: {
1763: GEN h;
1764: if (lg(f) == 1) return g;
1765: if (lg(g) == 1) return f;
1766: h = cgetg(3,t_MAT);
1767: h[1] = (long)concatsp((GEN)f[1], (GEN)g[1]);
1768: h[2] = (long)concatsp((GEN)f[2], (GEN)g[2]);
1769: return sort_factor_gen(h, cmpii);
1770: }
1771:
1772: GEN
1773: factor(GEN x)
1774: {
1775: long tx=typ(x),lx,av,tetpil,i,j,pa,v,r1;
1776: GEN y,p,p1,p2,p3,p4,p5,pol;
1777:
1778: if (is_matvec_t(tx))
1779: {
1780: lx=lg(x); y=cgetg(lx,tx);
1781: for (i=1; i<lx; i++) y[i]=(long)factor((GEN)x[i]);
1782: return y;
1783: }
1784: if (gcmp0(x))
1785: {
1786: y=cgetg(3,t_MAT);
1787: p1=cgetg(2,t_COL); y[1]=(long)p1; p1[1]=lcopy(x);
1788: p2=cgetg(2,t_COL); y[2]=(long)p2; p2[1]=un;
1789: return y;
1790: }
1791: av = avma;
1792: switch(tx)
1793: {
1794: case t_INT: return decomp(x);
1795:
1796: case t_FRACN:
1797: x = gred(x); /* fall through */
1798: case t_FRAC:
1799: p1 = decomp((GEN)x[1]);
1800: p2 = decomp((GEN)x[2]); p2[2] = (long)gneg_i((GEN)p2[2]);
1801: return gerepilecopy(av, merge_factor_i(p1,p2));
1802:
1803: case t_POL:
1804: tx=poltype(x,&p,&pol,&pa);
1805: switch(tx)
1806: {
1807: case 0: err(impl,"factor for general polynomials");
1808: case t_INT: return factpol(x,1);
1809: case t_INTMOD: return factmod(x,p);
1810:
1811: case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
1812: av = avma; p1 = roots(x,pa); tetpil = avma;
1813: p2=cgetg(lx,t_COL);
1814: for (i=1; i<lx; i++)
1815: p2[i] = (long)deg1pol_i(gun, gneg((GEN)p1[i]), v);
1816: y[1]=lpile(av,tetpil,p2);
1817: p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
1818: y[2]=(long)p3; return y;
1819:
1820: case t_REAL: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);
1821: av=avma; p1=roots(x,pa); tetpil=avma;
1822: for(r1=1; r1<lx; r1++)
1823: if (signe(gmael(p1,r1,2))) break;
1824: lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
1825: for(i=1; i<r1; i++)
1826: p2[i] = (long)deg1pol_i(gun, negr(gmael(p1,i,1)), v);
1827: for( ; i<lx; i++)
1828: {
1829: GEN a = (GEN) p1[2*i-r1];
1830: p = cgetg(5, t_POL); p2[i] = (long)p;
1831: p[1] = evalsigne(1) | evalvarn(v) | evallgef(5);
1832: p[2] = lnorm(a);
1833: p[3] = lmul2n((GEN)a[1],1); setsigne(p[3],-signe(p[3]));
1834: p[4] = un;
1835: }
1836: y[1]=lpile(av,tetpil,p2);
1837: p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;
1838: y[2]=(long)p3; return y;
1839:
1840: case t_PADIC: return factorpadic4(x,p,pa);
1841:
1842: default:
1843: {
1844: long killv;
1845: x = dummycopy(x); lx=lgef(x);
1846: pol = dummycopy(pol);
1847: v = manage_var(4,NULL);
1848: for(i=2; i<lx; i++)
1849: {
1850: p1=(GEN)x[i];
1851: switch(typ(p1))
1852: {
1853: case t_QUAD: p1++;
1854: case t_COMPLEX:
1855: p2 = cgetg(3, t_POLMOD); x[i] = (long) p2;
1856: p2[1] = (long)pol;
1857: p2[2] = (long)deg1pol_i((GEN)p1[2], (GEN)p1[1], v);
1858: }
1859: }
1860: killv = (avma != (ulong)pol);
1861: if (killv) setvarn(pol, fetch_var());
1862: switch (typ2(tx))
1863: {
1864: case t_INT: p1 = polfnf(x,pol); break;
1865: case t_INTMOD: p1 = factmod9(x,p,pol); break;
1866: default: err(impl,"factor of general polynomial");
1867: return NULL; /* not reached */
1868: }
1869: switch (typ1(tx))
1870: {
1871: case t_POLMOD:
1872: if (killv) delete_var();
1873: return gerepileupto(av,p1);
1874: case t_COMPLEX: p5 = gi; break;
1875: case t_QUAD: p5=cgetg(4,t_QUAD);
1876: p5[1]=(long)pol; p5[2]=zero; p5[3]=un;
1877: break;
1878: default: err(impl,"factor of general polynomial");
1879: return NULL; /* not reached */
1880: }
1881: p2=(GEN)p1[1];
1882: for(i=1; i<lg(p2); i++)
1883: {
1884: p3=(GEN)p2[i];
1885: for(j=2; j<lgef(p3); j++)
1886: {
1887: p4=(GEN)p3[j];
1888: if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5);
1889: }
1890: }
1891: if (killv) delete_var();
1892: tetpil=avma; y=cgetg(3,t_MAT);
1893: y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]);
1894: return gerepile(av,tetpil,y);
1895: }
1896: }
1897:
1898: case t_RFRACN:
1899: x=gred_rfrac(x); /* fall through */
1900: case t_RFRAC:
1901: p1=factor((GEN)x[1]);
1902: p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]);
1903: tetpil=avma; y=cgetg(3,t_MAT);
1904: y[1]=lconcat((GEN)p1[1],(GEN)p2[1]);
1905: y[2]=lconcat((GEN)p1[2],p3);
1906: return gerepile(av,tetpil,y);
1907: }
1908: err(talker,"can't factor %Z",x);
1909: return NULL; /* not reached */
1910: }
1911: #undef typ1
1912: #undef typ2
1913: #undef typs
1914: #undef tsh
1915:
1916: GEN
1917: divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
1918: {
1919: long i,k,lx = lg(x);
1920:
1921: if (lx == 1) return gun;
1922: if (lx == 2) return gcopy((GEN)x[1]);
1923: x = dummycopy(x); k = lx;
1924: while (k > 2)
1925: {
1926: if (DEBUGLEVEL>7)
1927: fprintferr("prod: remaining objects %ld\n",k-1);
1928: lx = k; k = 1;
1929: for (i=1; i<lx-1; i+=2)
1930: x[k++] = (long)mul((GEN)x[i],(GEN)x[i+1]);
1931: if (i < lx) x[k++] = x[i];
1932: }
1933: return (GEN)x[1];
1934: }
1935:
1936: static GEN static_nf;
1937:
1938: static GEN
1939: myidealmulred(GEN x, GEN y) { return idealmulred(static_nf, x, y, 0); }
1940:
1941: static GEN
1942: myidealpowred(GEN x, GEN n) { return idealpowred(static_nf, x, n, 0); }
1943:
1944: static GEN
1945: myidealmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }
1946:
1947: static GEN
1948: myidealpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }
1949:
1950: GEN
1951: factorback_i(GEN fa, GEN nf, int red)
1952: {
1953: long av=avma,k,l,lx,t=typ(fa);
1954: GEN ex, x;
1955: GEN (*_mul)(GEN,GEN);
1956: GEN (*_pow)(GEN,GEN);
1957: if (nf)
1958: {
1959: static_nf = nf;
1960: if (red)
1961: {
1962: _mul = &myidealmulred;
1963: _pow = &myidealpowred;
1964: }
1965: else
1966: {
1967: _mul = &myidealmul;
1968: _pow = &myidealpow;
1969: }
1970: }
1971: else
1972: {
1973: _mul = &gmul;
1974: _pow = &powgi;
1975: }
1976: if ( t == t_VEC || t == t_COL)
1977: return gerepileupto(av, divide_conquer_prod(fa, _mul));
1978: if (t!=t_MAT || lg(fa)!=3)
1979: err(talker,"not a factorisation in factorback");
1980: ex=(GEN)fa[2]; fa=(GEN)fa[1];
1981: lx = lg(fa); if (lx == 1) return gun;
1982: x = cgetg(lx,t_VEC);
1983: for (l=1,k=1; k<lx; k++)
1984: if (signe(ex[k]))
1985: x[l++] = (long)_pow((GEN)fa[k],(GEN)ex[k]);
1986: setlg(x,l);
1987: return gerepileupto(av, divide_conquer_prod(x, _mul));
1988: }
1989:
1990: GEN
1991: factorback(GEN fa, GEN nf)
1992: {
1993: return factorback_i(fa,nf,0);
1994: }
1995:
1996: GEN
1997: gisirreducible(GEN x)
1998: {
1999: long av=avma, tx = typ(x),l,i;
2000: GEN y;
2001:
2002: if (is_matvec_t(tx))
2003: {
2004: l=lg(x); y=cgetg(l,tx);
2005: for (i=1; i<l; i++) y[i]=(long)gisirreducible((GEN)x[i]);
2006: return y;
2007: }
2008: if (is_intreal_t(tx) || is_frac_t(tx)) return gzero;
2009: if (tx!=t_POL) err(notpoler,"gisirreducible");
2010: l=lgef(x); if (l<=3) return gzero;
2011: y=factor(x); avma=av;
2012: return (lgef(gcoeff(y,1,1))==l)?gun:gzero;
2013: }
2014:
2015: /*******************************************************************/
2016: /* */
2017: /* GENERIC GCD */
2018: /* */
2019: /*******************************************************************/
2020:
2021: GEN
2022: gcd0(GEN x, GEN y, long flag)
2023: {
2024: switch(flag)
2025: {
2026: case 0: return gassoc_proto(ggcd,x,y);
2027: case 1: return gassoc_proto(modulargcd,x,y);
2028: case 2: return gassoc_proto(srgcd,x,y);
2029: default: err(flagerr,"gcd");
2030: }
2031: return NULL; /* not reached */
2032: }
2033:
2034: /* x is a COMPLEX or a QUAD */
2035: static GEN
2036: triv_cont_gcd(GEN x, GEN y)
2037: {
2038: long av = avma, tetpil;
2039: GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2])
2040: : ggcd((GEN)x[2],(GEN)x[3]);
2041: tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
2042: }
2043:
2044: static GEN
2045: cont_gcd(GEN x, GEN y)
2046: {
2047: long av = avma, tetpil;
2048: GEN p1 = content(x);
2049:
2050: tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
2051: }
2052:
2053: /* y is a PADIC, x a rational number or an INTMOD */
2054: static GEN
2055: padic_gcd(GEN x, GEN y)
2056: {
2057: long v = ggval(x,(GEN)y[2]), w = valp(y);
2058: if (w < v) v = w;
2059: return gpuigs((GEN)y[2], v);
2060: }
2061:
2062: #define fix_frac(z) if (signe(z[2])<0)\
2063: {\
2064: setsigne(z[1],-signe(z[1]));\
2065: setsigne(z[2],1);\
2066: }
2067:
2068: GEN
2069: ggcd(GEN x, GEN y)
2070: {
2071: long l,av,tetpil,i,vx,vy, tx = typ(x), ty = typ(y);
2072: GEN p1,z;
2073:
2074: if (tx>ty) { swap(x,y); lswap(tx,ty); }
2075: if (is_matvec_t(ty))
2076: {
2077: l=lg(y); z=cgetg(l,ty);
2078: for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]);
2079: return z;
2080: }
2081: if (is_noncalc_t(tx) || is_noncalc_t(ty)) err(operf,"g",tx,ty);
2082: if (gcmp0(x)) return gcopy(y);
2083: if (gcmp0(y)) return gcopy(x);
2084: if (is_const_t(tx))
2085: {
2086: if (ty == t_FRACN)
2087: {
2088: if (tx==t_INTMOD)
2089: {
2090: av=avma; y = gred(y); tetpil=avma;
2091: return gerepile(av,tetpil,ggcd(x,y));
2092: }
2093: ty = t_FRAC;
2094: }
2095: if (tx == t_FRACN) tx = t_FRAC;
2096: if (ty == tx) switch(tx)
2097: {
2098: case t_INT:
2099: return mppgcd(x,y);
2100:
2101: case t_INTMOD: z=cgetg(3,t_INTMOD);
2102: if (egalii((GEN)x[1],(GEN)y[1]))
2103: { copyifstack(x[1],z[1]); }
2104: else
2105: z[1] = lmppgcd((GEN)x[1],(GEN)y[1]);
2106: if (gcmp1((GEN)z[1])) z[2] = zero;
2107: else
2108: {
2109: av = avma; p1 = mppgcd((GEN)z[1],(GEN)x[2]);
2110: if (!is_pm1(p1)) p1 = gerepileuptoint(av, mppgcd(p1,(GEN)y[2]));
2111: z[2] = (long)p1;
2112: }
2113: return z;
2114:
2115: case t_FRAC: z=cgetg(3,t_FRAC);
2116: z[1] = (long)mppgcd((GEN)x[1], (GEN)y[1]);
2117: z[2] = (long)mpppcm((GEN)x[2], (GEN)y[2]);
2118: return z;
2119:
2120: case t_COMPLEX:
2121: av=avma; p1=gdiv(x,y);
2122: if (gcmp0((GEN)p1[2]))
2123: {
2124: p1 = (GEN)p1[1];
2125: switch(typ(p1))
2126: {
2127: case t_INT:
2128: avma=av; return gcopy(y);
2129: case t_FRAC: case t_FRACN:
2130: tetpil=avma; return gerepile(av,tetpil,gdiv(y,(GEN)p1[2]));
2131: default: avma=av; return gun;
2132: }
2133: }
2134: if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) {avma=av; return gcopy(y);}
2135: p1 = ginv(p1); avma=av;
2136: if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(x);
2137: return triv_cont_gcd(y,x);
2138:
2139: case t_PADIC:
2140: if (!egalii((GEN)x[2],(GEN)y[2])) return gun;
2141: vx = valp(x);
2142: vy = valp(y); return gpuigs((GEN)y[2], min(vy,vx));
2143:
2144: case t_QUAD:
2145: av=avma; p1=gdiv(x,y);
2146: if (gcmp0((GEN)p1[3]))
2147: {
2148: p1=(GEN)p1[2];
2149: if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
2150: tetpil=avma; return gerepile(av,tetpil, gdiv(y,(GEN)p1[2]));
2151: }
2152: if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) {avma=av; return gcopy(y);}
2153: p1 = ginv(p1); avma=av;
2154: if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(x);
2155: return triv_cont_gcd(y,x);
2156:
2157: default: return gun; /* t_REAL */
2158: }
2159: if (is_const_t(ty)) switch(tx)
2160: {
2161: case t_INT:
2162: switch(ty)
2163: {
2164: case t_INTMOD: z = cgetg(3,t_INTMOD);
2165: copyifstack(y[1],z[1]); av=avma;
2166: p1 = mppgcd((GEN)y[1],(GEN)y[2]);
2167: if (!is_pm1(p1)) p1 = gerepileuptoint(av, mppgcd(x,p1));
2168: z[2] = (long)p1; return z;
2169:
2170: case t_FRAC: z = cgetg(3,t_FRAC);
2171: z[1] = lmppgcd(x,(GEN)y[1]);
2172: z[2] = licopy((GEN)y[2]); return z;
2173:
2174: case t_COMPLEX: case t_QUAD:
2175: return triv_cont_gcd(y,x);
2176:
2177: case t_PADIC:
2178: return padic_gcd(x,y);
2179:
2180: default: return gun; /* t_REAL */
2181: }
2182:
2183: case t_INTMOD:
2184: switch(ty)
2185: {
2186: case t_FRAC:
2187: av = avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); avma = av;
2188: if (!is_pm1(p1)) err(operi,"g",tx,ty);
2189: return ggcd((GEN)y[1], x);
2190:
2191: case t_COMPLEX: case t_QUAD:
2192: return triv_cont_gcd(y,x);
2193:
2194: case t_PADIC:
2195: return padic_gcd(x,y);
2196: }
2197:
2198: case t_FRAC:
2199: switch(ty)
2200: {
2201: case t_COMPLEX: case t_QUAD:
2202: return triv_cont_gcd(y,x);
2203:
2204: case t_PADIC:
2205: return padic_gcd(x,y);
2206: }
2207:
2208: case t_COMPLEX: /* ty = PADIC or QUAD */
2209: return triv_cont_gcd(x,y);
2210:
2211: case t_PADIC: /* ty = QUAD */
2212: return triv_cont_gcd(y,x);
2213:
2214: default: return gun; /* tx = t_REAL */
2215: }
2216: return cont_gcd(y,x);
2217: }
2218:
2219: vx = gvar9(x); vy = gvar9(y);
2220: if (vy<vx) return cont_gcd(y,x);
2221: if (vx<vy) return cont_gcd(x,y);
2222: switch(tx)
2223: {
2224: case t_POLMOD:
2225: switch(ty)
2226: {
2227: case t_POLMOD: z=cgetg(3,t_POLMOD);
2228: if (gegal((GEN)x[1],(GEN)y[1]))
2229: { copyifstack(x[1],z[1]); }
2230: else
2231: z[1] = lgcd((GEN)x[1],(GEN)y[1]);
2232: if (lgef(z[1])<=3) z[2]=zero;
2233: else
2234: {
2235: av=avma; p1=ggcd((GEN)z[1],(GEN)x[2]);
2236: if (lgef(p1)>3)
2237: {
2238: tetpil=avma;
2239: p1=gerepile(av,tetpil,ggcd(p1,(GEN)y[2]));
2240: }
2241: z[2]=(long)p1;
2242: }
2243: return z;
2244:
2245: case t_POL: z=cgetg(3,t_POLMOD);
2246: copyifstack(x[1],z[1]); av=avma;
2247: p1=ggcd((GEN)x[1],(GEN)x[2]);
2248: if (lgef(p1)>3)
2249: { tetpil=avma; p1=gerepile(av,tetpil,ggcd(y,p1)); }
2250: z[2]=(long)p1; return z;
2251:
2252: case t_RFRAC:
2253: av = avma; p1=ggcd((GEN)x[1],(GEN)y[2]); avma = av;
2254: if (!gcmp1(p1)) err(operi,"g",tx,ty);
2255: return ggcd((GEN)y[1],x);
2256:
2257: case t_RFRACN:
2258: av=avma; p1=gred_rfrac(y); tetpil=avma;
2259: return gerepile(av,tetpil,ggcd(p1,x));
2260: }
2261: break;
2262:
2263: case t_POL:
2264: switch(ty)
2265: {
2266: case t_POL:
2267: return srgcd(x,y);
2268:
2269: case t_SER:
2270: return gpuigs(polx[vx],min(valp(y),gval(x,vx)));
2271:
2272: case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty);
2273: z[1]=lgcd(x,(GEN)y[1]);
2274: z[2]=lcopy((GEN)y[2]); return z;
2275: }
2276: break;
2277:
2278: case t_SER:
2279: switch(ty)
2280: {
2281: case t_SER:
2282: return gpuigs(polx[vx],min(valp(x),valp(y)));
2283:
2284: case t_RFRAC: case t_RFRACN:
2285: return gpuigs(polx[vx],min(valp(x),gval(y,vx)));
2286: }
2287: break;
2288:
2289: case t_RFRAC: case t_RFRACN: z=cgetg(3,ty);
2290: if (!is_rfrac_t(ty)) err(operf,"g",tx,ty);
2291: p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2]));
2292: tetpil = avma;
2293: z[2] = lpile((long)z,tetpil,gmul(p1, (GEN)x[2]));
2294: z[1] = lgcd((GEN)x[1], (GEN)y[1]); return z;
2295: }
2296: err(operf,"g",tx,ty);
2297: return NULL; /* not reached */
2298: }
2299:
2300: GEN glcm0(GEN x, GEN y)
2301: {
2302: return gassoc_proto(glcm,x,y);
2303: }
2304:
2305: GEN
2306: glcm(GEN x, GEN y)
2307: {
2308: long av,tx,ty,i,l;
2309: GEN p1,p2,z;
2310:
2311: ty = typ(y);
2312: if (is_matvec_t(ty))
2313: {
2314: l=lg(y); z=cgetg(l,ty);
2315: for (i=1; i<l; i++) z[i]=(long)glcm(x,(GEN)y[i]);
2316: return z;
2317: }
2318: tx = typ(x);
2319: if (is_matvec_t(tx))
2320: {
2321: l=lg(x); z=cgetg(l,tx);
2322: for (i=1; i<l; i++) z[i]=(long)glcm((GEN)x[i],y);
2323: return z;
2324: }
2325: if (tx==t_INT && ty==t_INT) return mpppcm(x,y);
2326: if (gcmp0(x)) return gzero;
2327:
2328: av = avma;
2329: p1 = ggcd(x,y); if (!gcmp1(p1)) y = gdiv(y,p1);
2330: p2 = gmul(x,y);
2331: switch(typ(p2))
2332: {
2333: case t_INT: if (signe(p2)<0) setsigne(p2,1);
2334: break;
2335: case t_POL:
2336: if (lgef(p2) <= 2) break;
2337: p1=leading_term(p2);
2338: if (typ(p1)==t_INT && signe(p1)<0) p2=gneg(p2);
2339: }
2340: return gerepileupto(av,p2);
2341: }
2342:
2343: static GEN
2344: polgcdnun(GEN x, GEN y)
2345: {
2346: long av1, av = avma, lim = stack_lim(av,1);
2347: GEN p1, yorig = y;
2348:
2349: for(;;)
2350: {
2351: av1 = avma; p1 = gres(x,y);
2352: if (gcmp0(p1))
2353: {
2354: avma = av1;
2355: if (lgef(y) == 3 && isinexactreal(y)) { avma = av; return gun; }
2356: return (y==yorig)? gcopy(y): gerepileupto(av,y);
2357: }
2358: x = y; y = p1;
2359: if (low_stack(lim,stack_lim(av,1)))
2360: {
2361: GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y;
2362: if(DEBUGMEM>1) err(warnmem,"polgcdnun");
2363: gerepilemanysp(av,av1,gptr,2);
2364: }
2365: }
2366: }
2367:
2368: /* return 1 if coeff explosion is not possible */
2369: static int
2370: issimplefield(GEN x)
2371: {
2372: long lx,i;
2373: switch(typ(x))
2374: {
2375: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
2376: return 1;
2377: case t_POL:
2378: lx=lgef(x);
2379: for (i=2; i<lx; i++)
2380: if (issimplefield((GEN)x[i])) return 1;
2381: return 0;
2382: case t_COMPLEX: case t_POLMOD:
2383: return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]);
2384: }
2385: return 0;
2386: }
2387:
2388: int
2389: isrational(GEN x)
2390: {
2391: long i;
2392: for (i=lgef(x)-1; i>1; i--)
2393: switch(typ(x[i]))
2394: {
2395: case t_INT: case t_FRAC: break;
2396: default: return 0;
2397: }
2398: return 1;
2399: }
2400:
2401: static int
2402: can_use_modular_gcd(GEN x, GEN *mod, long *v)
2403: {
2404: GEN p1, p2;
2405: long i;
2406: for (i=lgef(x)-1; i>1; i--)
2407: {
2408: p1 = (GEN)x[i];
2409: switch(typ(p1))
2410: {
2411: case t_INT: case t_FRAC: break;
2412: case t_POLMOD:
2413: p2 = (GEN)p1[1];
2414: if (*mod)
2415: {
2416: if (!isrational((GEN)p1[2])) return 0;
2417: if (!gegal(*mod,p2)) return 0;
2418: } else
2419: {
2420: if (!isrational(p2)) return 0;
2421: *mod = p2; *v = varn(p2);
2422: }
2423: break;
2424: case t_POL:
2425: if (!isrational(p1)) return 0;
2426: if (*v >= 0)
2427: {
2428: if (*v != varn(p1)) return 0;
2429: } else *v = varn(p1);
2430: break;
2431: default: return 0;
2432: }
2433: }
2434: return 1;
2435: }
2436:
2437: static int
2438: isinexactfield(GEN x)
2439: {
2440: long lx,i;
2441: switch(typ(x))
2442: {
2443: case t_REAL: case t_PADIC: case t_SER:
2444: return 1;
2445: case t_POL:
2446: lx=lgef(x); if (lx == 2) return 0;/*true 0 polynomial*/
2447: for (i=2; i<lx; i++)
2448: if (isinexactfield((GEN)x[i])) return 1;
2449: return 0;
2450: case t_COMPLEX: case t_POLMOD:
2451: return isinexactfield((GEN)x[1]) || isinexactfield((GEN)x[2]);
2452: }
2453: return 0;
2454: }
2455:
2456: static GEN
2457: gcdmonome(GEN x, GEN y)
2458: {
2459: long tetpil,av=avma, dx=degpol(x), v=varn(x), e=gval(y,v);
2460: GEN p1,p2;
2461:
2462: if (dx < e) e = dx;
2463: p1=ggcd(leading_term(x),content(y)); p2=gpuigs(polx[v],e);
2464: tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));
2465: }
2466:
2467: /***********************************************************************/
2468: /** **/
2469: /** GENERIC EXTENDED GCD **/
2470: /** **/
2471: /***********************************************************************/
2472:
2473: static GEN
2474: polinvinexact(GEN x, GEN y)
2475: {
2476: ulong av = avma;
2477: long i,dx=degpol(x),dy=degpol(y),lz=dx+dy;
2478: GEN v,z;
2479:
2480: if (dx < 0 || dy < 0) err(talker,"non-invertible polynomial in polinvmod");
2481: z=cgetg(dy+2,t_POL); z[1]=y[1];
2482: v=cgetg(lz+1,t_COL);
2483: for (i=1; i<lz; i++) v[i]=zero;
2484: v[lz]=un; v=gauss(sylvestermatrix(y,x),v);
2485: for (i=2; i<dy+2; i++) z[i]=v[lz-i+2];
2486: z = normalizepol_i(z, dy+2);
2487: return gerepilecopy(av,z);
2488: }
2489:
2490: static GEN
2491: polinvmod(GEN x, GEN y)
2492: {
2493: long av,av1,tx,vx=varn(x),vy=varn(y);
2494: GEN u,v,d,p1;
2495:
2496: while (vx!=vy)
2497: {
2498: if (vx > vy)
2499: {
2500: d=cgetg(3,t_RFRAC);
2501: d[1]=(long)polun[vx];
2502: d[2]=lcopy(x); return d;
2503: }
2504: if (lgef(x)!=3) err(talker,"non-invertible polynomial in polinvmod");
2505: x=(GEN)x[2]; vx=gvar(x);
2506: }
2507: tx=typ(x);
2508: if (tx==t_POL)
2509: {
2510: if (isinexactfield(x) || isinexactfield(y))
2511: return polinvinexact(x,y);
2512:
2513: av=avma; d=subresext(x,y,&u,&v);
2514: if (gcmp0(d)) err(talker,"non-invertible polynomial in polinvmod");
2515: if (typ(d)==t_POL && varn(d)==vx)
2516: {
2517: if (lgef(d)>3) err(talker,"non-invertible polynomial in polinvmod");
2518: d=(GEN)d[2];
2519: }
2520: av1=avma; return gerepile(av,av1,gdiv(u,d));
2521: }
2522: if (!is_rfrac_t(tx)) err(typeer,"polinvmod");
2523: av=avma; p1=gmul((GEN)x[1],polinvmod((GEN)x[2],y));
2524: av1=avma; return gerepile(av,av1,gmod(p1,y));
2525: }
2526:
2527: GEN
2528: gbezout(GEN x, GEN y, GEN *u, GEN *v)
2529: {
2530: long tx=typ(x),ty=typ(y);
2531:
2532: if (tx==t_INT && ty==t_INT) return bezout(x,y,u,v);
2533: if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) err(typeer,"gbezout");
2534: return bezoutpol(x,y,u,v);
2535: }
2536:
2537: GEN
2538: vecbezout(GEN x, GEN y)
2539: {
2540: GEN z=cgetg(4,t_VEC);
2541: z[3]=(long)gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
2542: return z;
2543: }
2544:
2545: GEN
2546: vecbezoutres(GEN x, GEN y)
2547: {
2548: GEN z=cgetg(4,t_VEC);
2549: z[3]=(long)subresext(x,y,(GEN*)(z+1),(GEN*)(z+2));
2550: return z;
2551: }
2552:
2553: /*******************************************************************/
2554: /* */
2555: /* CONTENT / PRIMITIVE PART */
2556: /* */
2557: /*******************************************************************/
2558:
2559: GEN
2560: content(GEN x)
2561: {
2562: long lx,i,tx=typ(x);
2563: ulong av,tetpil;
2564: GEN p1,p2;
2565:
2566: if (is_scalar_t(tx))
2567: return tx==t_POLMOD? content((GEN)x[2]): gcopy(x);
2568: av = avma;
2569: switch(tx)
2570: {
2571: case t_RFRAC: case t_RFRACN:
2572: p1=content((GEN)x[1]);
2573: p2=content((GEN)x[2]); tetpil=avma;
2574: return gerepile(av,tetpil,gdiv(p1,p2));
2575:
2576: case t_VEC: case t_COL: case t_MAT:
2577: lx = lg(x); if (lx==1) return gun;
2578: p1=content((GEN)x[1]);
2579: for (i=2; i<lx; i++) p1 = ggcd(p1,content((GEN)x[i]));
2580: return gerepileupto(av,p1);
2581:
2582: case t_POL:
2583: if (!signe(x)) return gzero;
2584: lx = lgef(x); break;
2585: case t_SER:
2586: if (!signe(x)) return gzero;
2587: lx = lg(x); break;
2588: case t_QFR: case t_QFI:
2589: lx = 4; break;
2590:
2591: default: err(typeer,"content");
2592: return NULL; /* not reached */
2593: }
2594: for (i=lontyp[tx]; i<lx; i++)
2595: if (typ(x[i]) != t_INT) break;
2596: lx--; p1=(GEN)x[lx];
2597: if (i > lx)
2598: { /* integer coeffs */
2599: while (lx>lontyp[tx])
2600: {
2601: lx--; p1=mppgcd(p1,(GEN)x[lx]);
2602: if (is_pm1(p1)) { avma=av; return gun; }
2603: }
2604: }
2605: else
2606: {
2607: while (lx>lontyp[tx])
2608: {
2609: lx--; p1=ggcd(p1,(GEN)x[lx]);
2610: }
2611: if (isinexactreal(p1)) { avma=av; return gun; }
2612: }
2613: return av==avma? gcopy(p1): gerepileupto(av,p1);
2614: }
2615:
2616: GEN
2617: primitive_part(GEN x, GEN *c)
2618: {
2619: *c = content(x);
2620: if (gcmp1(*c)) *c = NULL; else x = gdiv(x,*c);
2621: return x;
2622: }
2623:
2624: /*******************************************************************/
2625: /* */
2626: /* SUBRESULTANT */
2627: /* */
2628: /*******************************************************************/
2629: /* for internal use */
2630: GEN
2631: gdivexact(GEN x, GEN y)
2632: {
2633: long i,lx;
2634: GEN z;
2635: if (gcmp1(y)) return x;
2636: switch(typ(x))
2637: {
2638: case t_INT:
2639: if (typ(y)==t_INT) return diviiexact(x,y);
2640: if (!signe(x)) return gzero;
2641: break;
2642: case t_INTMOD:
2643: case t_POLMOD: return gmul(x,ginv(y));
2644: case t_POL:
2645: switch(typ(y))
2646: {
2647: case t_INTMOD:
2648: case t_POLMOD: return gmul(x,ginv(y));
2649: case t_POL:
2650: if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES_EXACT);
2651: }
2652: lx = lgef(x); z = cgetg(lx,t_POL);
2653: for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);
2654: z[1]=x[1]; return z;
2655: case t_VEC: case t_COL: case t_MAT:
2656: lx = lg(x); z = cgetg(lx, typ(x));
2657: for (i=1; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);
2658: return z;
2659: }
2660: if (DEBUGLEVEL) err(warner,"missing case in gdivexact");
2661: return gdiv(x,y);
2662: }
2663:
2664: static GEN
2665: init_resultant(GEN x, GEN y)
2666: {
2667: long tx,ty;
2668: if (gcmp0(x) || gcmp0(y)) return gzero;
2669: tx = typ(x); ty = typ(y);
2670: if (is_scalar_t(tx) || is_scalar_t(ty))
2671: {
2672: if (tx==t_POL) return gpuigs(y,degpol(x));
2673: if (ty==t_POL) return gpuigs(x,degpol(y));
2674: return gun;
2675: }
2676: if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall");
2677: if (varn(x)==varn(y)) return NULL;
2678: return (varn(x)<varn(y))? gpuigs(y,degpol(x)): gpuigs(x,degpol(y));
2679: }
2680:
2681: /* return coefficients s.t x = x_0 X^n + ... + x_n */
2682: static GEN
2683: revpol(GEN x)
2684: {
2685: long i,n = degpol(x);
2686: GEN y = cgetg(n+3,t_POL);
2687: y[1] = x[1]; x += 2; y += 2;
2688: for (i=0; i<=n; i++) y[i] = x[n-i];
2689: return y;
2690: }
2691:
2692: /* assume dx >= dy, y non constant */
2693: GEN
2694: pseudorem(GEN x, GEN y)
2695: {
2696: long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,lx, p;
2697:
2698: if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
2699: (void)new_chunk(2);
2700: dx=degpol(x); x = revpol(x);
2701: dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;
2702: av2 = avma; lim = stack_lim(av2,1);
2703: for (;;)
2704: {
2705: x[0] = lneg((GEN)x[0]); p--;
2706: for (i=1; i<=dy; i++)
2707: x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
2708: for ( ; i<=dx; i++)
2709: x[i] = lmul((GEN)y[0], (GEN)x[i]);
2710: do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0]));
2711: if (dx < dy) break;
2712: if (low_stack(lim,stack_lim(av2,1)))
2713: {
2714: if(DEBUGMEM>1) err(warnmem,"pseudorem dx = %ld >= %ld",dx,dy);
2715: gerepilemanycoeffs(av2,x,dx+1);
2716: }
2717: }
2718: if (dx < 0) return zeropol(vx);
2719: lx = dx+3; x -= 2;
2720: x[0]=evaltyp(t_POL) | evallg(lx);
2721: x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);
2722: x = revpol(x) - 2;
2723: return gerepileupto(av, gmul(x, gpowgs((GEN)y[0], p)));
2724: }
2725:
2726: extern void gerepilemanycoeffs2(long av, GEN x, long n, GEN y, long o);
2727:
2728: /* assume dx >= dy, y non constant
2729: * Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2730: GEN
2731: pseudodiv(GEN x, GEN y, GEN *ptr)
2732: {
2733: long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,iz,lx,lz,p;
2734: GEN z, r, ypow;
2735:
2736: if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
2737: (void)new_chunk(2);
2738: dx=degpol(x); x = revpol(x);
2739: dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;
2740: lz = dz+3; z = cgetg(lz, t_POL) + 2;
2741: ypow = new_chunk(dz+1);
2742: ypow[0] = un;
2743: for (i=1; i<=dz; i++) ypow[i] = lmul((GEN)ypow[i-1], (GEN)y[0]);
2744: av2 = avma; lim = stack_lim(av2,1);
2745: for (iz=0;;)
2746: {
2747: p--;
2748: z[iz++] = lmul((GEN)x[0], (GEN)ypow[p]);
2749: x[0] = lneg((GEN)x[0]);
2750: for (i=1; i<=dy; i++)
2751: x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));
2752: for ( ; i<=dx; i++)
2753: x[i] = lmul((GEN)y[0], (GEN)x[i]);
2754: x++; dx--;
2755: while (dx >= dy && gcmp0((GEN)x[0])) { x++; dx--; z[iz++] = zero; }
2756: if (dx < dy) break;
2757: if (low_stack(lim,stack_lim(av2,1)))
2758: {
2759: if(DEBUGMEM>1) err(warnmem,"pseudodiv dx = %ld >= %ld",dx,dy);
2760: gerepilemanycoeffs2(av2,x,dx+1, z,iz);
2761: }
2762: }
2763: while (dx >= 0 && gcmp0((GEN)x[0])) { x++; dx--; }
2764: if (dx < 0)
2765: x = zeropol(vx);
2766: else
2767: {
2768: lx = dx+3; x -= 2;
2769: x[0] = evaltyp(t_POL) | evallg(lx);
2770: x[1] = evalsigne(1) | evalvarn(vx) | evallgef(lx);
2771: x = revpol(x) - 2;
2772: }
2773:
2774: z -= 2;
2775: z[0] = evaltyp(t_POL) | evallg(lz);
2776: z[1] = evalsigne(1) | evalvarn(vx) | evallgef(lz);
2777: z = revpol(z) - 2;
2778: r = gmul(x, (GEN)ypow[p]);
2779: {
2780: GEN *gptr[2]; gptr[0] = &z; gptr[1] = &r;
2781: gerepilemany(av,gptr,2);
2782: }
2783: *ptr = r; return z;
2784: }
2785:
2786: /* Return resultant(u,v). If sol != NULL: set *sol to the last non-zero
2787: * polynomial in the prs IF the sequence was computed, and gzero otherwise
2788: */
2789: GEN
2790: subresall(GEN u, GEN v, GEN *sol)
2791: {
2792: ulong av,av2,lim;
2793: long degq,dx,dy,du,dv,dr,signh;
2794: GEN z,g,h,r,p1,p2,cu,cv;
2795:
2796: if (sol) *sol=gzero;
2797: if ((r = init_resultant(u,v))) return r;
2798:
2799: if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v);
2800: dx=lgef(u); dy=lgef(v); signh=1;
2801: if (dx<dy)
2802: {
2803: swap(u,v); lswap(dx,dy);
2804: if (both_even(dx, dy)) signh = -signh;
2805: }
2806: if (dy==3) return gpowgs((GEN)v[2],dx-3);
2807: av = avma;
2808: u = primitive_part(u, &cu);
2809: v = primitive_part(v, &cv);
2810: g = h = gun; av2 = avma; lim = stack_lim(av2,1);
2811: for(;;)
2812: {
2813: r = pseudorem(u,v); dr = lgef(r);
2814: if (dr == 2)
2815: {
2816: if (sol) { avma = (long)(r+2); *sol=gerepileupto(av,v); } else avma = av;
2817: return gzero;
2818: }
2819: du = lgef(u); dv = lgef(v); degq = du-dv;
2820: u = v; p1 = g; g = leading_term(u);
2821: switch(degq)
2822: {
2823: case 0: break;
2824: case 1:
2825: p1 = gmul(h,p1); h = g; break;
2826: default:
2827: p1 = gmul(gpowgs(h,degq),p1);
2828: h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
2829: }
2830: if (both_even(du,dv)) signh = -signh;
2831: v = gdivexact(r,p1);
2832: if (dr==3) break;
2833: if (low_stack(lim,stack_lim(av2,1)))
2834: {
2835: GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
2836: if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);
2837: gerepilemany(av2,gptr,4);
2838: }
2839: }
2840: z = (GEN)v[2];
2841: if (dv > 4) z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4));
2842: if (signh < 0) z = gneg(z); /* z = resultant(ppart(x), ppart(y)) */
2843: p2 = gun;
2844: if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));
2845: if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));
2846: z = gmul(z, p2);
2847:
2848: if (sol) u = gclone(u);
2849: z = gerepileupto(av, z);
2850: if (sol) { *sol = forcecopy(u); gunclone(u); }
2851: return z;
2852: }
2853:
2854: static GEN
2855: scalar_res(GEN x, GEN y, GEN *U, GEN *V)
2856: {
2857: long dx=lgef(x)-4;
2858: *V=gpuigs(y,dx); *U=gzero; return gmul(y,*V);
2859: }
2860:
2861: /* compute U, V s.t Ux + Vy = resultant(x,y) */
2862: GEN
2863: subresext(GEN x, GEN y, GEN *U, GEN *V)
2864: {
2865: ulong av,av2,tetpil,lim;
2866: long degq,tx,ty,dx,dy,du,dv,dr,signh;
2867: GEN q,z,g,h,r,p1,p2,cu,cv,u,v,um1,uze,vze,xprim,yprim;
2868:
2869: if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; }
2870: tx=typ(x); ty=typ(y);
2871: if (is_scalar_t(tx) || is_scalar_t(ty))
2872: {
2873: if (tx==t_POL) return scalar_res(x,y,U,V);
2874: if (ty==t_POL) return scalar_res(y,x,V,U);
2875: *U = ginv(x); *V = gzero; return gun;
2876: }
2877: if (tx!=t_POL || ty!=t_POL) err(typeer,"subresext");
2878: if (varn(x) != varn(y))
2879: return (varn(x)<varn(y))? scalar_res(x,y,U,V)
2880: : scalar_res(y,x,V,U);
2881: dx = lgef(x); dy = lgef(y); signh = 1;
2882: if (dx < dy)
2883: {
2884: pswap(U,V); lswap(dx,dy); swap(x,y);
2885: if (both_even(dx, dy)) signh = -signh;
2886: }
2887: if (dy == 3)
2888: {
2889: *V = gpowgs((GEN)y[2],dx-4);
2890: *U = gzero; return gmul(*V,(GEN)y[2]);
2891: }
2892: av = avma;
2893: u = primitive_part(x, &cu); xprim = u;
2894: v = primitive_part(y, &cv); yprim = v;
2895: g = h = gun; av2 = avma; lim = stack_lim(av2,1);
2896: um1 = gun; uze = gzero;
2897: for(;;)
2898: {
2899: q = pseudodiv(u,v, &r); dr = lgef(r);
2900: if (dr == 2) { *U = *V = gzero; avma = av; return gzero; }
2901:
2902: du = lgef(u); dv = lgef(v); degq = du-dv;
2903: /* lead(v)^(degq + 1) * um1 - q * uze */
2904: p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));
2905: um1 = uze; uze = p1;
2906: u = v; p1 = g; g = leading_term(u);
2907: switch(degq)
2908: {
2909: case 0: break;
2910: case 1: p1 = gmul(h,p1); h = g; break;
2911: default:
2912: p1 = gmul(gpowgs(h,degq),p1);
2913: h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
2914: }
2915: if (both_even(du, dv)) signh = -signh;
2916: v = gdivexact(r,p1);
2917: uze= gdivexact(uze,p1);
2918: if (dr==3) break;
2919: if (low_stack(lim,stack_lim(av2,1)))
2920: {
2921: GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
2922: gptr[4]=&uze; gptr[5]=&um1;
2923: if(DEBUGMEM>1) err(warnmem,"subresext, dr = %ld",dr);
2924: gerepilemany(av2,gptr,6);
2925: }
2926: }
2927: z = (GEN)v[2];
2928: if (dv > 4)
2929: {
2930: /* z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4)); */
2931: p1 = gpowgs(gdiv(z,h),dv-4);
2932: z = gmul(z,p1); uze = gmul(uze,p1);
2933: }
2934: if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); }
2935: p1 = gadd(z, gneg(gmul(uze,xprim)));
2936: if (!poldivis(p1,yprim,&vze)) err(bugparier,"subresext");
2937: /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
2938: p2 = gun;
2939: if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));
2940: if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));
2941: cu = cu? gdiv(p2,cu): p2;
2942: cv = cv? gdiv(p2,cv): p2;
2943:
2944: tetpil = avma;
2945: z = gmul(z,p2); uze = gmul(uze,cu); vze = gmul(vze,cv);
2946: {
2947: GEN *gptr[3]; gptr[0]=&z; gptr[1]=&uze; gptr[2]=&vze;
2948: gerepilemanysp(av,tetpil,gptr,3);
2949: }
2950: *U = uze; *V = vze; return z;
2951: }
2952:
2953: static GEN
2954: scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
2955: {
2956: *U=gzero; *V=ginv(y); return polun[varn(x)];
2957: }
2958:
2959: static GEN
2960: zero_bezout(GEN y, GEN *U, GEN *V)
2961: {
2962: GEN x=content(y);
2963: *U=gzero; *V = ginv(x); return gmul(y,*V);
2964: }
2965:
2966: /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
2967: GEN
2968: bezoutpol(GEN x, GEN y, GEN *U, GEN *V)
2969: {
2970: ulong av,av2,tetpil,lim;
2971: long degq,tx,ty,dx,dy,du,dv,dr;
2972: GEN q,z,g,h,r,p1,cu,cv,u,v,um1,uze,vze,xprim,yprim, *gptr[3];
2973:
2974: if (gcmp0(x)) return zero_bezout(y,U,V);
2975: if (gcmp0(y)) return zero_bezout(x,V,U);
2976: tx=typ(x); ty=typ(y); av=avma;
2977: if (is_scalar_t(tx) || is_scalar_t(ty))
2978: {
2979: if (tx==t_POL) return scalar_bezout(x,y,U,V);
2980: if (ty==t_POL) return scalar_bezout(y,x,V,U);
2981: *U = ginv(x); *V = gzero; return polun[0];
2982: }
2983: if (tx!=t_POL || ty!=t_POL) err(typeer,"bezoutpol");
2984: if (varn(x)!=varn(y))
2985: return (varn(x)<varn(y))? scalar_bezout(x,y,U,V)
2986: : scalar_bezout(y,x,V,U);
2987: dx = lgef(x); dy = lgef(y);
2988: if (dx < dy)
2989: {
2990: pswap(U,V); lswap(dx,dy); swap(x,y);
2991: }
2992: if (dy==3) return scalar_bezout(x,y,U,V);
2993:
2994: u = primitive_part(x, &cu); xprim = u;
2995: v = primitive_part(y, &cv); yprim = v;
2996: g = h = gun; av2 = avma; lim = stack_lim(av2,1);
2997: um1 = gun; uze = gzero;
2998: for(;;)
2999: {
3000: q = pseudodiv(u,v, &r); dr = lgef(r);
3001: if (dr == 2) break;
3002:
3003: du = lgef(u); dv = lgef(v); degq = du-dv;
3004: p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));
3005: um1 = uze; uze = p1;
3006: u = v; p1 = g; g = leading_term(u);
3007: switch(degq)
3008: {
3009: case 0: break;
3010: case 1:
3011: p1 = gmul(h,p1); h = g; break;
3012: default:
3013: p1 = gmul(gpowgs(h,degq), p1);
3014: h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
3015: }
3016: v = gdivexact(r,p1);
3017: uze= gdivexact(uze,p1);
3018: if (dr==3) break;
3019: if (low_stack(lim,stack_lim(av2,1)))
3020: {
3021: GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
3022: gptr[4]=&uze; gptr[5]=&um1;
3023: if(DEBUGMEM>1) err(warnmem,"bezoutpol, dr = %ld",dr);
3024: gerepilemany(av2,gptr,6);
3025: }
3026: }
3027: if (!poldivis(gsub(v,gmul(uze,xprim)),yprim, &vze))
3028: err(warner,"inexact computation in bezoutpol");
3029: if (cu) uze = gdiv(uze,cu);
3030: if (cv) vze = gdiv(vze,cv);
3031: p1 = ginv(content(v)); tetpil = avma;
3032:
3033: uze = gmul(uze,p1);
3034: vze = gmul(vze,p1);
3035: z = gmul(v,p1);
3036: gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&z;
3037: gerepilemanysp(av,tetpil,gptr,3);
3038: *U = uze; *V = vze; return z;
3039: }
3040:
3041: /*******************************************************************/
3042: /* */
3043: /* RESULTANT USING DUCOS VARIANT */
3044: /* */
3045: /*******************************************************************/
3046:
3047: static GEN
3048: reductum(GEN P)
3049: {
3050: if (signe(P)==0) return P;
3051: return normalizepol_i(dummycopy(P),lgef(P)-1);
3052: }
3053:
3054: static GEN
3055: Lazard(GEN x, GEN y, long n)
3056: {
3057: long a, b;
3058: GEN c;
3059:
3060: if (n<=1) return x;
3061: a=1; while (n >= (b=a+a)) a=b;
3062: c=x; n-=a;
3063: while (a>1)
3064: {
3065: a>>=1; c=gdivexact(gsqr(c),y);
3066: if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
3067: }
3068: return c;
3069: }
3070:
3071: static GEN
3072: Lazard2(GEN F, GEN x, GEN y, long n)
3073: {
3074: if (n<=1) return F;
3075: x = Lazard(x,y,n-1); return gdivexact(gmul(x,F),y);
3076: }
3077:
3078: extern GEN addshiftpol(GEN x, GEN y, long d);
3079: #define addshift(x,y) addshiftpol((x),(y),1)
3080:
3081: static GEN
3082: nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
3083: {
3084: GEN p0, q0, z0, H, A;
3085: long p, q, j, av, lim, v = varn(P);
3086:
3087: z0 = leading_term(Z);
3088: p = degpol(P); p0 = leading_term(P); P = reductum(P);
3089: q = degpol(Q); q0 = leading_term(Q); Q = reductum(Q);
3090:
3091: av = avma; lim = stack_lim(av,1);
3092: H = gneg(reductum(Z));
3093: A = gmul((GEN)P[q+2],H);
3094: for (j = q+1; j < p; j++)
3095: {
3096: H = (degpol(H) == q-1) ?
3097: addshift(reductum(H), gdivexact(gmul(gneg((GEN)H[q+1]),Q), q0)) :
3098: addshift(H, zeropol(v));
3099: A = gadd(A,gmul((GEN)P[j+2],H));
3100: if (low_stack(lim,stack_lim(av,1)))
3101: {
3102: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&H;
3103: if(DEBUGMEM>1) err(warnmem,"nextSousResultant j = %ld/%ld",j,p);
3104: gerepilemany(av,gptr,2);
3105: }
3106: }
3107: P = normalizepol_i(P, q+2);
3108: A = gdivexact(gadd(A,gmul(z0,P)), p0);
3109: A = (degpol(H) == q-1) ?
3110: gadd(gmul(q0,addshift(reductum(H),A)), gmul(gneg((GEN)H[q+1]),Q)) :
3111: gmul(q0, addshift(H,A));
3112: return gdivexact(A, ((p-q)&1)? s: gneg(s));
3113: }
3114:
3115: GEN
3116: resultantducos(GEN P, GEN Q)
3117: {
3118: ulong av = avma, av2,lim;
3119: long dP,dQ,delta;
3120: GEN cP,cQ,Z,s;
3121:
3122: if ((Z = init_resultant(P,Q))) return Z;
3123: dP = degpol(P);
3124: dQ = degpol(Q);
3125: P = primitive_part(P, &cP);
3126: Q = primitive_part(Q, &cQ);
3127: delta = dP - dQ;
3128: if (delta < 0)
3129: {
3130: Z = (dP & dQ & 1)? gneg(Q): Q;
3131: Q = P; P = Z; delta = -delta;
3132: }
3133: s = gun;
3134: if (degpol(Q) > 0)
3135: {
3136: av2 = avma; lim = stack_lim(av2,1);
3137: s = gpuigs(leading_term(Q),delta);
3138: Z = Q;
3139: Q = pseudorem(P, gneg(Q));
3140: P = Z;
3141: while(degpol(Q) > 0)
3142: {
3143: if (low_stack(lim,stack_lim(av,1)))
3144: {
3145: GEN *gptr[2]; gptr[0]=&P; gptr[1]=&Q;
3146: if(DEBUGMEM>1) err(warnmem,"resultantducos, degpol Q = %ld",degpol(Q));
3147: gerepilemany(av2,gptr,2); s = leading_term(P);
3148: }
3149: delta = degpol(P) - degpol(Q);
3150: Z = Lazard2(Q, leading_term(Q), s, delta);
3151: Q = nextSousResultant(P, Q, Z, s);
3152: P = Z;
3153: s = leading_term(P);
3154: }
3155: }
3156: if (!signe(Q)) { avma = av; return gzero; }
3157: if (!degpol(P)){ avma = av; return gun; }
3158: s = Lazard(leading_term(Q), s, degpol(P));
3159: if (cP) s = gmul(s, gpowgs(cP,dQ));
3160: if (cQ) s = gmul(s, gpowgs(cQ,dP)); else if (!cP) s = gcopy(s);
3161: return gerepileupto(av, s);
3162: }
3163:
3164: /*******************************************************************/
3165: /* */
3166: /* RESULTANT USING SYLVESTER MATRIX */
3167: /* */
3168: /*******************************************************************/
3169: static GEN
3170: _zeropol(void)
3171: {
3172: GEN x = cgetg(3,t_POL);
3173: x[1] = evallgef(3);
3174: x[2] = zero; return x;
3175: }
3176:
3177: static GEN
3178: sylvester_col(GEN x, long j, long d, long D)
3179: {
3180: GEN c = cgetg(d+1,t_COL);
3181: long i;
3182: for (i=1; i< j; i++) c[i]=zero;
3183: for ( ; i<=D; i++) c[i]=x[D-i+2];
3184: for ( ; i<=d; i++) c[i]=zero;
3185: return c;
3186: }
3187:
3188: GEN
3189: sylvestermatrix_i(GEN x, GEN y)
3190: {
3191: long j,d,dx,dy;
3192: GEN M;
3193:
3194: dx = degpol(x); if (dx < 0) { dx = 0; x = _zeropol(); }
3195: dy = degpol(y); if (dy < 0) { dy = 0; y = _zeropol(); }
3196: d = dx+dy; M = cgetg(d+1,t_MAT);
3197: for (j=1; j<=dy; j++) M[j] = (long)sylvester_col(x,j,d,j+dx);
3198: for (j=1; j<=dx; j++) M[j+dy] = (long)sylvester_col(y,j,d,j+dy);
3199: return M;
3200: }
3201:
3202: GEN
3203: sylvestermatrix(GEN x, GEN y)
3204: {
3205: long i,j,lx;
3206: if (typ(x)!=t_POL || typ(y)!=t_POL) err(typeer,"sylvestermatrix");
3207: if (varn(x) != varn(y))
3208: err(talker,"not the same variables in sylvestermatrix");
3209: x = sylvestermatrix_i(x,y); lx = lg(x);
3210: for (i=1; i<lx; i++)
3211: for (j=1; j<lx; j++) coeff(x,i,j) = lcopy(gcoeff(x,i,j));
3212: return x;
3213: }
3214:
3215: GEN
3216: resultant2(GEN x, GEN y)
3217: {
3218: long av;
3219: GEN r;
3220: if ((r = init_resultant(x,y))) return r;
3221: av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
3222: }
3223:
3224: static GEN
3225: fix_pol(GEN x, long v, long *mx)
3226: {
3227: long vx;
3228: GEN p1;
3229:
3230: if (typ(x) == t_POL)
3231: {
3232: vx = varn(x);
3233: if (vx)
3234: {
3235: if (v>=vx) return gsubst(x,v,polx[0]);
3236: p1 = cgetg(3,t_POL);
3237: p1[1] = evalvarn(0)|evalsigne(signe(x))|evallgef(3);
3238: p1[2] = (long)x; return p1;
3239: }
3240: if (v)
3241: {
3242: *mx = 1;
3243: return gsubst(gsubst(x,0,polx[MAXVARN]),v,polx[0]);
3244: }
3245: }
3246: return x;
3247: }
3248:
3249: /* resultant of x and y with respect to variable v, or with respect to their
3250: * main variable if v < 0. */
3251: GEN
3252: polresultant0(GEN x, GEN y, long v, long flag)
3253: {
3254: long av = avma, m = 0;
3255:
3256: if (v >= 0)
3257: {
3258: x = fix_pol(x,v, &m);
3259: y = fix_pol(y,v, &m);
3260: }
3261: switch(flag)
3262: {
3263: case 0: x=subresall(x,y,NULL); break;
3264: case 1: x=resultant2(x,y); break;
3265: case 2: x=resultantducos(x,y); break;
3266: default: err(flagerr,"polresultant");
3267: }
3268: if (m) x = gsubst(x,MAXVARN,polx[0]);
3269: return gerepileupto(av,x);
3270: }
3271:
3272: /*******************************************************************/
3273: /* */
3274: /* GCD USING SUBRESULTANT */
3275: /* */
3276: /*******************************************************************/
3277: extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
3278: extern GEN to_polmod(GEN x, GEN mod);
3279:
3280: GEN
3281: srgcd(GEN x, GEN y)
3282: {
3283: long av,av1,tetpil,tx=typ(x),ty=typ(y),dx,dy,vx,vmod,lim;
3284: GEN d,g,h,p1,p2,u,v,mod,cx,cy;
3285:
3286: if (!signe(y)) return gcopy(x);
3287: if (!signe(x)) return gcopy(y);
3288: if (is_scalar_t(tx) || is_scalar_t(ty)) return gun;
3289: if (tx!=t_POL || ty!=t_POL) err(typeer,"srgcd");
3290: vx=varn(x); if (vx!=varn(y)) return gun;
3291: if (ismonome(x)) return gcdmonome(x,y);
3292: if (ismonome(y)) return gcdmonome(y,x);
3293: av = avma;
3294: mod = NULL; vmod = -1;
3295: if (can_use_modular_gcd(x, &mod, &vmod) &&
3296: can_use_modular_gcd(y, &mod, &vmod))
3297: {
3298: if (mod)
3299: { /* (Q[Y]/mod)[X]*/
3300: x = primitive_part(lift(x), &cx);
3301: y = primitive_part(lift(y), &cy);
3302: if (cx)
3303: { if (cy) cx = ggcd(cx,cy); }
3304: else
3305: cx = cy;
3306: d = nfgcd(x,y,mod,NULL); if (cx) d = gmul(d,cx);
3307: return gerepileupto(av, gmul(d, to_polmod(gun,mod)));
3308: }
3309: if (vmod < 0) return modulargcd(x,y); /* Q[X] */
3310: }
3311:
3312: if (issimplefield(x) || issimplefield(y)) x = polgcdnun(x,y);
3313: else
3314: {
3315: dx=lgef(x); dy=lgef(y);
3316: if (dx<dy) { swap(x,y); lswap(dx,dy); }
3317: p1=content(x); p2=content(y); d=ggcd(p1,p2);
3318:
3319: tetpil=avma; d=gmul(d,polun[vx]);
3320: if (dy==3) return gerepile(av,tetpil,d);
3321:
3322: av1=avma; lim=stack_lim(av1,1);
3323: u=gdiv(x,p1); v=gdiv(y,p2); g=h=gun;
3324: for(;;)
3325: {
3326: GEN r = pseudorem(u,v);
3327: long degq, du, dv, dr=lgef(r);
3328:
3329: if (dr <= 3)
3330: {
3331: if (gcmp0(r)) break;
3332: avma=av1; return gerepile(av,tetpil,d);
3333: }
3334: if (DEBUGLEVEL > 9) fprintferr("srgcd: dr = %ld\n", dr);
3335: du=lgef(u); dv=lgef(v); degq=du-dv; u=v;
3336: switch(degq)
3337: {
3338: case 0:
3339: v = gdiv(r,g);
3340: g = leading_term(u);
3341: break;
3342: case 1:
3343: v = gdiv(r,gmul(h,g));
3344: h = g = leading_term(u);
3345: break;
3346: default:
3347: v = gdiv(r,gmul(gpuigs(h,degq),g));
3348: g = leading_term(u);
3349: h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));
3350: }
3351: if (low_stack(lim, stack_lim(av1,1)))
3352: {
3353: GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
3354: if(DEBUGMEM>1) err(warnmem,"srgcd");
3355: gerepilemany(av1,gptr,4);
3356: }
3357: }
3358: p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);
3359: x = gmul(d,v);
3360: }
3361:
3362: if (typ(x)!=t_POL) x = gmul(polun[vx],x);
3363: else
3364: {
3365: p1=leading_term(x); ty=typ(p1);
3366: if ((is_frac_t(ty) || is_intreal_t(ty)) && gsigne(p1)<0) x = gneg(x);
3367: }
3368: return gerepileupto(av,x);
3369: }
3370:
3371: extern GEN qf_disc(GEN x, GEN y, GEN z);
3372:
3373: GEN
3374: poldisc0(GEN x, long v)
3375: {
3376: long av,tx=typ(x),i;
3377: GEN z,p1,p2;
3378:
3379: switch(tx)
3380: {
3381: case t_POL:
3382: if (gcmp0(x)) return gzero;
3383: av = avma; i = 0;
3384: if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
3385: p1 = subres(x, derivpol(x));
3386: p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
3387: if (degpol(x) & 2) p1 = gneg(p1);
3388: if (i) p1 = gsubst(p1, MAXVARN, polx[0]);
3389: return gerepileupto(av,p1);
3390:
3391: case t_COMPLEX:
3392: return stoi(-4);
3393:
3394: case t_QUAD: case t_POLMOD:
3395: return poldisc0((GEN)x[1], v);
3396:
3397: case t_QFR: case t_QFI:
3398: av = avma; return gerepileuptoint(av, qf_disc(x,NULL,NULL));
3399:
3400: case t_VEC: case t_COL: case t_MAT:
3401: i=lg(x); z=cgetg(i,tx);
3402: for (i--; i; i--) z[i]=(long)poldisc0((GEN)x[i], v);
3403: return z;
3404: }
3405: err(typeer,"discsr");
3406: return NULL; /* not reached */
3407: }
3408:
3409: GEN
3410: discsr(GEN x)
3411: {
3412: return poldisc0(x, -1);
3413: }
3414:
3415: GEN
3416: reduceddiscsmith(GEN pol)
3417: {
3418: long av=avma,tetpil,i,j,n;
3419: GEN polp,alpha,p1,m;
3420:
3421: if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith");
3422: n=degpol(pol);
3423: if (n<=0) err(constpoler,"reduceddiscsmith");
3424: check_pol_int(pol,"poldiscreduced");
3425: if (!gcmp1((GEN)pol[n+2]))
3426: err(talker,"non-monic polynomial in poldiscreduced");
3427: polp = derivpol(pol); alpha = polx[varn(pol)];
3428: m=cgetg(n+1,t_MAT);
3429: for (j=1; j<=n; j++)
3430: {
3431: p1=cgetg(n+1,t_COL); m[j]=(long)p1;
3432: for (i=1; i<=n; i++) p1[i]=(long)truecoeff(polp,i-1);
3433: if (j<n) polp = gres(gmul(alpha,polp), pol);
3434: }
3435: tetpil=avma; return gerepile(av,tetpil,smith(m));
3436: }
3437:
3438: /***********************************************************************/
3439: /** **/
3440: /** STURM ALGORITHM **/
3441: /** (number of real roots of x in ]a,b]) **/
3442: /** **/
3443: /***********************************************************************/
3444:
3445: /* if a (resp. b) is NULL, set it to -oo (resp. +oo) */
3446: long
3447: sturmpart(GEN x, GEN a, GEN b)
3448: {
3449: long av = avma, lim = stack_lim(av,1), sl,sr,s,t,r1;
3450: GEN g,h,u,v;
3451:
3452: if (typ(x) != t_POL) err(typeer,"sturm");
3453: if (gcmp0(x)) err(zeropoler,"sturm");
3454: s=lgef(x); if (s==3) return 0;
3455:
3456: sl = gsigne(leading_term(x));
3457: if (s==4)
3458: {
3459: t = a? gsigne(poleval(x,a)): -sl;
3460: if (t == 0) { avma = av; return 0; }
3461: s = b? gsigne(poleval(x,b)): sl;
3462: avma = av; return (s == t)? 0: 1;
3463: }
3464: u=gdiv(x,content(x)); v=derivpol(x);
3465: v=gdiv(v,content(v));
3466: g=gun; h=gun;
3467: s = b? gsigne(poleval(u,b)): sl;
3468: t = a? gsigne(poleval(u,a)): ((lgef(u)&1)? sl: -sl);
3469: r1=0;
3470: sr = b? gsigne(poleval(v,b)): s;
3471: if (sr)
3472: {
3473: if (!s) s=sr;
3474: else if (sr!=s) { s= -s; r1--; }
3475: }
3476: sr = a? gsigne(poleval(v,a)): -t;
3477: if (sr)
3478: {
3479: if (!t) t=sr;
3480: else if (sr!=t) { t= -t; r1++; }
3481: }
3482: for(;;)
3483: {
3484: GEN p1, r = pseudorem(u,v);
3485: long du=lgef(u), dv=lgef(v), dr=lgef(r), degq=du-dv;
3486:
3487: if (dr<=2) err(talker,"not a squarefree polynomial in sturm");
3488: if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
3489: sl = gsigne((GEN) r[dr-1]);
3490: sr = b? gsigne(poleval(r,b)): sl;
3491: if (sr)
3492: {
3493: if (!s) s=sr;
3494: else if (sr!=s) { s= -s; r1--; }
3495: }
3496: sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
3497: if (sr)
3498: {
3499: if (!t) t=sr;
3500: else if (sr!=t) { t= -t; r1++; }
3501: }
3502: if (dr==3) { avma=av; return r1; }
3503:
3504: u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
3505: switch(degq)
3506: {
3507: case 0: break;
3508: case 1:
3509: p1 = gmul(h,p1); h = g; break;
3510: default:
3511: p1 = gmul(gpuigs(h,degq),p1);
3512: h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));
3513: }
3514: v = gdivexact(r,p1);
3515: if (low_stack(lim,stack_lim(av,1)))
3516: {
3517: GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;
3518: if(DEBUGMEM>1) err(warnmem,"polsturm, dr = %ld",dr);
3519: gerepilemany(av,gptr,4);
3520: }
3521: }
3522: }
3523:
3524: /*******************************************************************/
3525: /* */
3526: /* QUADRATIC POLYNOMIAL ASSOCIATED TO A DISCRIMINANT */
3527: /* */
3528: /*******************************************************************/
3529:
3530: GEN
3531: quadpoly0(GEN x, long v)
3532: {
3533: long res,l,tetpil,i,sx, tx = typ(x);
3534: GEN y,p1;
3535:
3536: if (is_matvec_t(tx))
3537: {
3538: l=lg(x); y=cgetg(l,tx);
3539: for (i=1; i<l; i++) y[i]=(long)quadpoly0((GEN)x[i],v);
3540: return y;
3541: }
3542: if (tx!=t_INT) err(arither1);
3543: if (v < 0) v = 0;
3544: sx=signe(x);
3545: if (!sx) err(talker,"zero discriminant in quadpoly");
3546: y=cgetg(5,t_POL);
3547: y[1]=evalsigne(1) | evalvarn(v) | evallgef(5); y[4]=un;
3548: res=mod4(x); if (sx<0 && res) res=4-res;
3549: if (res>1) err(funder2,"quadpoly");
3550:
3551: l=avma; p1=shifti(x,-2); setsigne(p1,-signe(p1));
3552: y[2] = (long) p1;
3553: if (!res) y[3] = zero;
3554: else
3555: {
3556: if (sx<0) { tetpil=avma; y[2] = lpile(l,tetpil,addsi(1,p1)); }
3557: y[3] = lnegi(gun);
3558: }
3559: return y;
3560: }
3561:
3562: GEN
3563: quadpoly(GEN x)
3564: {
3565: return quadpoly0(x,-1);
3566: }
3567:
3568: GEN
3569: quadgen(GEN x)
3570: {
3571: GEN y=cgetg(4,t_QUAD);
3572: y[1]=lquadpoly(x); y[2]=zero; y[3]=un; return y;
3573: }
3574:
3575: /*******************************************************************/
3576: /* */
3577: /* GENERIC (modular) INVERSE */
3578: /* */
3579: /*******************************************************************/
3580:
3581: GEN
3582: ginvmod(GEN x, GEN y)
3583: {
3584: long tx=typ(x);
3585:
3586: switch(typ(y))
3587: {
3588: case t_POL:
3589: if (tx==t_POL) return polinvmod(x,y);
3590: if (is_scalar_t(tx)) return ginv(x);
3591: break;
3592: case t_INT:
3593: if (tx==t_INT) return mpinvmod(x,y);
3594: if (tx==t_POL) return gzero;
3595: }
3596: err(typeer,"ginvmod");
3597: return NULL; /* not reached */
3598: }
3599:
3600: /***********************************************************************/
3601: /** **/
3602: /** NEWTON POLYGON **/
3603: /** **/
3604: /***********************************************************************/
3605:
3606: /* assume leading coeff of x is non-zero */
3607: GEN
3608: newtonpoly(GEN x, GEN p)
3609: {
3610: GEN y;
3611: long n,ind,a,b,c,u1,u2,r1,r2;
3612: long *vval, num[] = {evaltyp(t_INT) | m_evallg(3), 0, 0};
3613:
3614: if (typ(x)!=t_POL) err(typeer,"newtonpoly");
3615: n=degpol(x); if (n<=0) { y=cgetg(1,t_VEC); return y; }
3616: y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
3617: vval = (long *) gpmalloc(sizeof(long)*(n+1));
3618: for (a=0; a<=n; a++) vval[a] = ggval((GEN)x[a],p);
3619: for (a=0, ind=1; a<n; a++)
3620: {
3621: if (vval[a] != VERYBIGINT) break;
3622: y[ind++] = lstoi(VERYBIGINT);
3623: }
3624: for (b=a+1; b<=n; a=b, b=a+1)
3625: {
3626: while (vval[b] == VERYBIGINT) b++;
3627: u1=vval[a]-vval[b]; u2=b-a;
3628: for (c=b+1; c<=n; c++)
3629: {
3630: if (vval[c] == VERYBIGINT) continue;
3631: r1=vval[a]-vval[c]; r2=c-a;
3632: if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
3633: }
3634: while (ind<=b) { affsi(u1,num); y[ind++] = ldivgs(num,u2); }
3635: }
3636: free(vval); return y;
3637: }
3638:
3639: extern int cmp_pol(GEN x, GEN y);
3640: extern GEN ZY_ZXY_resultant(GEN A, GEN B0, long *lambda);
3641: GEN matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom);
3642: GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
3643:
3644: /* Factor polynomial a on the number field defined by polynomial t */
3645: GEN
3646: polfnf(GEN a, GEN t)
3647: {
3648: ulong av = avma;
3649: GEN x0,y,p1,p2,u,fa,n,unt,dent,alift;
3650: long lx,i,k,e;
3651: int sqfree, tmonic;
3652:
3653: if (typ(a)!=t_POL || typ(t)!=t_POL) err(typeer,"polfnf");
3654: if (gcmp0(a)) return gcopy(a);
3655: a = fix_relative_pol(t,a,0);
3656: alift = lift(a);
3657: p1 = content(alift); if (!gcmp1(p1)) { a = gdiv(a, p1); alift = lift(a); }
3658: p1 = content(t); if (!gcmp1(t)) t = gdiv(t, p1);
3659: tmonic = is_pm1(leading_term(t));
3660:
3661: dent = ZX_disc(t); unt = gmodulsg(1,t);
3662: u = nfgcd(alift,derivpol(alift), t, dent);
3663: sqfree = gcmp1(u);
3664: u = sqfree? alift: lift_intern(gdiv(a, gmul(unt,u)));
3665: k = 0; n = ZY_ZXY_resultant(t, u, &k);
3666: if (DEBUGLEVEL > 4) fprintferr("polfnf: choosing k = %ld\n",k);
3667:
3668: /* n guaranteed to be squarefree */
3669: fa = squff2(n,0); lx = lg(fa);
3670: y = cgetg(3,t_MAT);
3671: p1 = cgetg(lx,t_COL); y[1] = (long)p1;
3672: p2 = cgetg(lx,t_COL); y[2] = (long)p2;
3673: x0 = gadd(polx[varn(a)], gmulsg(-k, gmodulcp(polx[varn(t)], t)));
3674: for (i=lx-1; i>1; i--)
3675: {
3676: GEN b, F = lift_intern(poleval((GEN)fa[i], x0));
3677: if (!tmonic) F = gdiv(F, content(F));
3678: F = nfgcd(u, F, t, dent);
3679: if (typ(F) != t_POL || degpol(F) == 0)
3680: err(talker,"reducible modulus in factornf");
3681: F = gmul(unt, F);
3682: F = gdiv(F, leading_term(F));
3683: u = lift_intern(gdeuc(u, F));
3684: u = gdiv(u, content(u));
3685: if (sqfree) e = 1;
3686: else
3687: {
3688: e = 0; while (poldivis(a,F, &b)) { a = b; e++; }
3689: if (degpol(a) == degpol(u)) sqfree = 1;
3690: }
3691: p1[i] = (long)F;
3692: p2[i] = lstoi(e);
3693: }
3694: u = gmul(unt, u);
3695: u = gdiv(u, leading_term(u));
3696: p1[1] = (long)u;
3697: e = sqfree? 1: degpol(a)/degpol(u);
3698: p2[1] = lstoi(e);
3699: return gerepilecopy(av, sort_factor(y, cmp_pol));
3700: }
3701:
3702: extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);
3703: extern GEN FpM(GEN z, GEN p);
3704: extern GEN polpol_to_mat(GEN v, long n);
3705: extern GEN mat_to_polpol(GEN x, long v, long w);
3706:
3707: static GEN
3708: to_frac(GEN a, GEN b)
3709: {
3710: GEN f = cgetg(3, t_FRAC);
3711: f[1] = (long)a;
3712: f[2] = (long)b; return f;
3713: }
3714:
3715: static GEN
3716: lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN denom)
3717: {
3718: GEN a, b;
3719: if (signe(t) < 0) t = addii(t, mod); /* in case t is a centerlift */
3720: if (!ratlift(t, mod, &a, &b, amax, bmax)
3721: || (denom && !divise(denom,b))
3722: || !gcmp1(mppgcd(a,b))) return NULL;
3723: if (!is_pm1(b)) a = to_frac(a, b);
3724: return a;
3725: }
3726:
3727: /* compute rational lifting for all the components of M modulo mod.
3728: * If one components fails, return NULL.
3729: * See ratlift.
3730: * If denom is not NULL, check that the denominators divide denom
3731: *
3732: * FIXME: NOT stack clean ! a & b stay on the stack.
3733: * If we suppose mod and denom coprime, then a and b are coprime
3734: * so we can do a cgetg(t_FRAC).
3735: */
3736: GEN
3737: matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom)
3738: {
3739: ulong ltop = avma;
3740: GEN N, a;
3741: long l2, l3;
3742: long i, j;
3743: if (typ(M)!=t_MAT) err(typeer,"matratlift");
3744: l2 = lg(M); l3 = lg((GEN)M[1]);
3745: N = cgetg(l2, t_MAT);
3746: for (j = 1; j < l2; ++j)
3747: {
3748: N[j] = lgetg(l3, t_COL);
3749: for (i = 1; i < l3; ++i)
3750: {
3751: a = lift_to_frac(gcoeff(M,i,j), mod, amax,bmax,denom);
3752: if (!a) { avma = ltop; return NULL; }
3753: coeff(N, i, j) = (long)a;
3754: }
3755: }
3756: return N;
3757: }
3758:
3759: GEN
3760: polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)
3761: {
3762: ulong ltop = avma;
3763: GEN Q,a;
3764: long l2, j;
3765: if (typ(P)!=t_POL) err(typeer,"polratlift");
3766: l2 = lg(P);
3767: Q = cgetg(l2, t_POL); Q[1] = P[1];
3768: for (j = 2; j < l2; ++j)
3769: {
3770: a = lift_to_frac((GEN)P[j], mod, amax,bmax,denom);
3771: if (!a) { avma = ltop; return NULL; }
3772: Q[j] = (long)a;
3773: }
3774: return Q;
3775: }
3776:
3777: /* P,Q in Z[X,Y], nf in Z[Y] irreducible. compute GCD in Q[Y]/(nf)[X].
3778: *
3779: * We essentially follow M. Encarnacion "On a modular Algorithm for computing
3780: * GCDs of polynomials over number fields" (ISSAC'94).
3781: *
3782: * We procede as follows
3783: * 1:compute the gcd modulo primes discarding bad primes as they are detected.
3784: * 2:reconstruct the result via matratlift, stoping as soon as we get weird
3785: * denominators.
3786: * 3:if matratlift succeeds, try the full division.
3787: * Suppose we does not have sufficient accuracy to get the result right:
3788: * it is extremly rare that matratlift will succeed, and even if it does, the
3789: * polynomial we get has sensible coefficients, so the full division will
3790: * not be too costly.
3791: *
3792: * FIXME: Handle rational coefficient for P and Q.
3793: * If not NULL, den must a a multiple of the denominator of the gcd,
3794: * for example the discriminant of nf.
3795: *
3796: * NOTE: if nf is not irreducible, nfgcd may loop forever, especially if the
3797: * gcd divides nf !
3798: */
3799: GEN
3800: nfgcd(GEN P, GEN Q, GEN nf, GEN den)
3801: {
3802: ulong ltop = avma;
3803: GEN sol, mod = NULL;
3804: long x = varn(P);
3805: long y = varn(nf);
3806: long d = degpol(nf);
3807: GEN lP, lQ;
3808: if (!signe(P) || !signe(Q))
3809: return zeropol(x);
3810: /*Compute denominators*/
3811: if (!den) den = ZX_disc(nf);
3812: lP = leading_term(P);
3813: lQ = leading_term(Q);
3814: if ( !((typ(lP)==t_INT && is_pm1(lP)) || (typ(lQ)==t_INT && is_pm1(lQ))) )
3815: den = mulii(den, mppgcd(ZX_resultant(lP, nf), ZX_resultant(lQ, nf)));
3816: /*Modular GCD*/
3817: {
3818: ulong btop = avma;
3819: ulong st_lim = stack_lim(btop, 1);
3820: long p;
3821: long dM=0, dR;
3822: GEN M, dsol, dens;
3823: GEN R, ax, bo;
3824: byteptr primepointer;
3825: for (p = 27449, primepointer = diffptr + 3000; ; p += *(primepointer++))
3826: {
3827: if (!*primepointer) err(primer1);
3828: if (!smodis(den, p))
3829: continue;/*Discard primes dividing the disc(T) or leadingcoeff(PQ) */
3830: if (DEBUGLEVEL>=5) fprintferr("nfgcd: p=%d\n",p);
3831: if ((R = FpXQX_safegcd(P, Q, nf, stoi(p))) == NULL)
3832: continue;/*Discard primes when modular gcd does not exist*/
3833: dR = degpol(R);
3834: if (dR == 0) return scalarpol(gun, x);
3835: if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
3836:
3837: R = polpol_to_mat(R, d);
3838: /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
3839: if (!mod || dR < dM) { M = R; mod = stoi(p); dM = dR; continue; }
3840: if (low_stack(st_lim, stack_lim(btop, 1)))
3841: {
3842: GEN *bptr[2]; bptr[0]=&M; bptr[1]=&mod;
3843: if (DEBUGMEM>1) err(warnmem,"nfgcd");
3844: gerepilemany(btop, bptr, 2);
3845: }
3846:
3847: ax = gmulgs(mpinvmod(stoi(p), mod), p);
3848: M = gadd(R, gmul(ax, gsub(M, R)));
3849: mod = mulis(mod, p);
3850: M = lift(FpM(M, mod));
3851: /* I suspect it must be better to take amax > bmax*/
3852: bo = racine(shifti(mod, -1));
3853: if ((sol = matratlift(M, mod, bo, bo, den)) == NULL)
3854: continue;
3855: dens = denom(sol);
3856: sol = mat_to_polpol(sol,x,y);
3857: dsol = gmul(sol, gmodulcp(dens, nf));
3858: if (gdivise(P, dsol) && gdivise(Q, dsol))
3859: break;
3860: }
3861: }
3862: return gerepilecopy(ltop, sol);
3863: }
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-2.2 / src / basemath