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