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