Annotation of OpenXM_contrib/pari/src/basemath/polarit3.c, Revision 1.1
1.1 ! maekawa 1: /***********************************************************************/
! 2: /** **/
! 3: /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
! 4: /** (third part) **/
! 5: /** **/
! 6: /***********************************************************************/
! 7: /* $Id: polarit3.c,v 1.1.1.1 1999/09/16 13:47:36 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: /*******************************************************************/
! 11: /* */
! 12: /* KARATSUBA (for polynomials) */
! 13: /* */
! 14: /*******************************************************************/
! 15: #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
! 16:
! 17: #if 1 /* for tunings */
! 18: long SQR_LIMIT = 6;
! 19: long MUL_LIMIT = 10;
! 20:
! 21: void
! 22: setsqpol(long a) { SQR_LIMIT=a; }
! 23: void
! 24: setmulpol(long a) { MUL_LIMIT=a; }
! 25:
! 26: GEN
! 27: specpol(GEN x, long nx)
! 28: {
! 29: GEN z = cgetg(nx+2,t_POL);
! 30: long i;
! 31: for (i=0; i<nx; i++) z[i+2] = x[i];
! 32: z[1]=evalsigne(1)|evallgef(nx+2);
! 33: return z;
! 34: }
! 35: #else
! 36: # define SQR_LIMIT 6
! 37: # define MUL_LIMIT 10
! 38: #endif
! 39:
! 40: static GEN
! 41: addpol(GEN x, GEN y, long lx, long ly)
! 42: {
! 43: long i,lz;
! 44: GEN z;
! 45:
! 46: if (ly>lx) swapspec(x,y, lx,ly);
! 47: lz = lx+2; z = cgetg(lz,t_POL) + 2;
! 48: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 49: for ( ; i<lx; i++) z[i]=x[i];
! 50: z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 51: }
! 52:
! 53: static GEN
! 54: addpolcopy(GEN x, GEN y, long lx, long ly)
! 55: {
! 56: long i,lz;
! 57: GEN z;
! 58:
! 59: if (ly>lx) swapspec(x,y, lx,ly);
! 60: lz = lx+2; z = cgetg(lz,t_POL) + 2;
! 61: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 62: for ( ; i<lx; i++) z[i]=lcopy((GEN)x[i]);
! 63: z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 64: }
! 65:
! 66: #ifdef INLINE
! 67: INLINE
! 68: #endif
! 69: GEN
! 70: mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b)
! 71: {
! 72: GEN p1 = NULL;
! 73: long i,av = avma;
! 74: for (i=a; i<b; i++)
! 75: if (ynonzero[i])
! 76: {
! 77: GEN p2 = gmul((GEN)y[i],(GEN)x[-i]);
! 78: p1 = p1 ? gadd(p1, p2): p2;
! 79: }
! 80: return p1 ? gerepileupto(av, p1): gzero;
! 81: }
! 82:
! 83: static GEN
! 84: mulpol(GEN x, GEN y, long nx, long ny)
! 85: {
! 86: long i,lz,nz;
! 87: GEN z;
! 88: char *p1;
! 89:
! 90: if (!ny) return zeropol(0);
! 91: lz = nx+ny+1; nz = lz-2;
! 92: z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
! 93: p1 = gpmalloc(ny);
! 94: for (i=0; i<ny; i++)
! 95: {
! 96: p1[i] = !isexactzero((GEN)y[i]);
! 97: z[i] = (long)mulpol_limb(x+i,y,p1,0,i+1);
! 98: }
! 99: for ( ; i<nx; i++) z[i] = (long)mulpol_limb(x+i,y,p1,0,ny);
! 100: for ( ; i<nz; i++) z[i] = (long)mulpol_limb(x+i,y,p1,i-nx+1,ny);
! 101: free(p1); z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 102: }
! 103:
! 104: /* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
! 105: GEN
! 106: addshiftw(GEN x, GEN y, long d)
! 107: {
! 108: GEN xd,yd,zd = (GEN)avma;
! 109: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
! 110:
! 111: x += 2; y += 2; a = ny-d;
! 112: if (a <= 0)
! 113: {
! 114: lz = (a>nx)? ny+2: nx+d+2;
! 115: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
! 116: while (xd > x) *--zd = *--xd;
! 117: x = zd + a;
! 118: while (zd > x) *--zd = zero;
! 119: }
! 120: else
! 121: {
! 122: xd = new_chunk(d); yd = y+d;
! 123: x = addpol(x,yd, nx,a);
! 124: lz = (a>nx)? ny+2: lgef(x)+d;
! 125: x += 2; while (xd > x) *--zd = *--xd;
! 126: }
! 127: while (yd > y) *--zd = *--yd;
! 128: *--zd = evalsigne(1) | evallgef(lz);
! 129: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
! 130: }
! 131:
! 132: GEN
! 133: addshiftpol(GEN x, GEN y, long d)
! 134: {
! 135: long v = varn(x);
! 136: if (!signe(x)) return y;
! 137: x = addshiftw(x,y,d);
! 138: setvarn(x,v); return x;
! 139: }
! 140:
! 141: /* as above, producing a clean stack */
! 142: static GEN
! 143: addshiftwcopy(GEN x, GEN y, long d)
! 144: {
! 145: GEN xd,yd,zd = (GEN)avma;
! 146: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
! 147:
! 148: x += 2; y += 2; a = ny-d;
! 149: if (a <= 0)
! 150: {
! 151: lz = nx+d+2;
! 152: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
! 153: while (xd > x) *--zd = lcopy((GEN)*--xd);
! 154: x = zd + a;
! 155: while (zd > x) *--zd = zero;
! 156: }
! 157: else
! 158: {
! 159: xd = new_chunk(d); yd = y+d;
! 160: x = addpolcopy(x,yd, nx,a);
! 161: lz = (a>nx)? ny+2: lgef(x)+d;
! 162: x += 2; while (xd > x) *--zd = *--xd;
! 163: }
! 164: while (yd > y) *--zd = lcopy((GEN)*--yd);
! 165: *--zd = evalsigne(1) | evallgef(lz);
! 166: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
! 167: }
! 168:
! 169: /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
! 170: * b+2 were sent instead. na, nb = number of terms of a, b.
! 171: * Only c, c0, c1, c2 are genuine GEN.
! 172: */
! 173: GEN
! 174: quickmul(GEN a, GEN b, long na, long nb)
! 175: {
! 176: GEN a0,c,c0;
! 177: long av,n0,n0a,i;
! 178:
! 179: if (na < nb) swapspec(a,b, na,nb);
! 180: if (nb < MUL_LIMIT) return mulpol(a,b,na,nb);
! 181: i=(na>>1); n0=na-i; na=i;
! 182: av=avma; a0=a+n0; n0a=n0;
! 183: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
! 184:
! 185: if (nb > n0)
! 186: {
! 187: GEN b0,c1,c2;
! 188: long n0b;
! 189:
! 190: nb -= n0; b0 = b+n0; n0b = n0;
! 191: while (n0b && isexactzero((GEN)b[n0b-1])) n0b--;
! 192: c = quickmul(a,b,n0a,n0b);
! 193: c0 = quickmul(a0,b0, na,nb);
! 194:
! 195: c2 = addpol(a0,a, na,n0a);
! 196: c1 = addpol(b0,b, nb,n0b);
! 197:
! 198: c1 = quickmul(c1+2,c2+2, lgef(c1)-2,lgef(c2)-2);
! 199: c2 = gneg_i(gadd(c0,c));
! 200: c0 = addshiftw(c0, gadd(c1,c2), n0);
! 201: }
! 202: else
! 203: {
! 204: c = quickmul(a,b,n0a,nb);
! 205: c0 = quickmul(a0,b,na,nb);
! 206: }
! 207: c0 = addshiftwcopy(c0,c,n0);
! 208: return gerepileupto(av,c0);
! 209: }
! 210:
! 211: GEN
! 212: sqrpol(GEN x, long nx)
! 213: {
! 214: long av,i,j,l,lz,nz;
! 215: GEN p1,z;
! 216: char *p2;
! 217:
! 218: if (!nx) return zeropol(0);
! 219: lz = (nx << 1) + 1, nz = lz-2;
! 220: z = cgetg(lz,t_POL) + 2;
! 221: p2 = gpmalloc(nx);
! 222: for (i=0; i<nx; i++)
! 223: {
! 224: p2[i] = !isexactzero((GEN)x[i]);
! 225: p1=gzero; av=avma; l=(i+1)>>1;
! 226: for (j=0; j<l; j++)
! 227: if (p2[j] && p2[i-j])
! 228: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
! 229: p1 = gshift(p1,1);
! 230: if ((i&1) == 0 && p2[i>>1])
! 231: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
! 232: z[i] = lpileupto(av,p1);
! 233: }
! 234: for ( ; i<nz; i++)
! 235: {
! 236: p1=gzero; av=avma; l=(i+1)>>1;
! 237: for (j=i-nx+1; j<l; j++)
! 238: if (p2[j] && p2[i-j])
! 239: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
! 240: p1 = gshift(p1,1);
! 241: if ((i&1) == 0 && p2[i>>1])
! 242: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
! 243: z[i] = lpileupto(av,p1);
! 244: }
! 245: free(p2); z -= 2; z[1]=0; return normalizepol_i(z,lz);
! 246: }
! 247:
! 248: GEN
! 249: quicksqr(GEN a, long na)
! 250: {
! 251: GEN a0,c,c0,c1;
! 252: long av,n0,n0a,i;
! 253:
! 254: if (na<SQR_LIMIT) return sqrpol(a,na);
! 255: i=(na>>1); n0=na-i; na=i;
! 256: av=avma; a0=a+n0; n0a=n0;
! 257: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
! 258:
! 259: c = quicksqr(a,n0a);
! 260: c0 = quicksqr(a0,na);
! 261: c1 = gmul2n(quickmul(a0,a, na,n0a), 1);
! 262: c0 = addshiftw(c0,c1, n0);
! 263: c0 = addshiftwcopy(c0,c,n0);
! 264: return gerepileupto(av,c0);
! 265: }
! 266:
! 267: /* x,pol in Z[X], p in Z, n in N, compute lift(x^n mod (p, pol)) */
! 268: GEN
! 269: Fp_pow_mod_pol(GEN x, GEN n, GEN pol, GEN p)
! 270: {
! 271: long m,i,j,av=avma, lim=stack_lim(av,1), vx = varn(x);
! 272: GEN p1 = n+2, y = x, z;
! 273: if (!signe(n)) return polun[vx];
! 274: if (is_pm1(n)) return gcopy(x);
! 275: m = *p1;
! 276: j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
! 277: for (i=lgefint(n)-2;;)
! 278: {
! 279: for (; j; m<<=1,j--)
! 280: {
! 281: z = quicksqr(y+2, lgef(y)-2);
! 282: y = Fp_pol_red(z, p);
! 283: y = Fp_res(y,pol, p);
! 284: if (low_stack(lim, stack_lim(av,1)))
! 285: {
! 286: if(DEBUGMEM>1) err(warnmem,"[1]: Fp_pow_mod_pol");
! 287: y = gerepileupto(av, y);
! 288: }
! 289: if (m<0)
! 290: {
! 291: z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2);
! 292: y = Fp_pol_red(z, p);
! 293: y = Fp_res(y,pol, p);
! 294: }
! 295: if (low_stack(lim, stack_lim(av,1)))
! 296: {
! 297: if(DEBUGMEM>1) err(warnmem,"[2]: Fp_pow_mod_pol");
! 298: y = gerepileupto(av, y);
! 299: }
! 300: }
! 301: if (--i == 0) break;
! 302: m = *++p1, j = BITS_IN_LONG;
! 303: }
! 304: setvarn(y,vx); return gerepileupto(av,y);
! 305: }
! 306:
! 307: int ff_poltype(GEN *x, GEN *p, GEN *pol);
! 308:
! 309: /* z in Z[X], return z * Mod(1,p), normalized*/
! 310: GEN
! 311: Fp_pol(GEN z, GEN p)
! 312: {
! 313: long i,l = lgef(z);
! 314: GEN a,x = cgetg(l,t_POL);
! 315: if (isonstack(p)) p = icopy(p);
! 316: for (i=2; i<l; i++)
! 317: {
! 318: a = cgetg(3,t_INTMOD); x[i] = (long)a;
! 319: a[1] = (long)p;
! 320: a[2] = lmodii((GEN)z[i],p);
! 321: }
! 322: x[1] = z[1]; return normalizepol_i(x,l);
! 323: }
! 324:
! 325: /* z in Z^n, return z * Mod(1,p), normalized*/
! 326: GEN
! 327: Fp_vec(GEN z, GEN p)
! 328: {
! 329: long i,l = lg(z);
! 330: GEN a,x = cgetg(l,typ(z));
! 331: if (isonstack(p)) p = icopy(p);
! 332: for (i=1; i<l; i++)
! 333: {
! 334: a = cgetg(3,t_INTMOD); x[i] = (long)a;
! 335: a[1] = (long)p;
! 336: a[2] = lmodii((GEN)z[i],p);
! 337: }
! 338: return x;
! 339: }
! 340:
! 341: /* z in Z[X], return lift(z * Mod(1,p)), normalized*/
! 342: GEN
! 343: Fp_pol_red(GEN z, GEN p)
! 344: {
! 345: long i,l = lgef(z);
! 346: GEN x = cgetg(l,t_POL);
! 347: for (i=2; i<l; i++) x[i] = lmodii((GEN)z[i],p);
! 348: x[1] = z[1]; return normalizepol_i(x,l);
! 349: }
! 350:
! 351: /* z in Z^n, return lift(z * Mod(1,p)) */
! 352: GEN
! 353: Fp_vec_red(GEN z, GEN p)
! 354: {
! 355: long i,l = lg(z);
! 356: GEN x = cgetg(l,typ(z));
! 357: for (i=1; i<l; i++) x[i] = lmodii((GEN)z[i],p);
! 358: return x;
! 359: }
! 360:
! 361: /* no garbage collection, divide by leading coeff, mod p */
! 362: GEN
! 363: normalize_mod_p(GEN z, GEN p)
! 364: {
! 365: long l = lgef(z)-1;
! 366: GEN p1 = (GEN)z[l]; /* leading term */
! 367: if (gcmp1(p1)) return z;
! 368: z = gmul(z, mpinvmod(p1,p));
! 369: return Fp_pol_red(z, p);
! 370: }
! 371:
! 372: /* as above, p is guaranteed small, and coeffs of z are C longs in [0,p-1],
! 373: * coeffs are in z[0..l-1] (instead of z[2] for regular pols)
! 374: * Set varn(z) = 0
! 375: */
! 376: GEN
! 377: Fp_pol_small(GEN z, GEN p, long l)
! 378: {
! 379: long i;
! 380: GEN a,x = cgetg(l,t_POL);
! 381: if (isonstack(p)) p = icopy(p);
! 382: if (is_bigint(p)) err(talker, "not a small prime in Fp_pol_small");
! 383: z -= 2;
! 384: for (i=2; i<l; i++) {
! 385: a = cgetg(3,t_INTMOD); x[i] = (long)a;
! 386: a[1] = (long)p;
! 387: a[2] = lstoi(z[i]);
! 388: }
! 389: return normalizepol_i(x,l);
! 390: }
! 391:
! 392: /* assume z[i] % p has been done. But we may have z[i] < 0 */
! 393: GEN
! 394: small_to_pol(GEN z, long l, long p)
! 395: {
! 396: GEN x = cgetg(l,t_POL);
! 397: long i;
! 398: z -= 2; for (i=2; i<l; i++) x[i] = lstoi(z[i]<0? p+z[i]: z[i]);
! 399: return normalizepol_i(x,l);
! 400: }
! 401:
! 402: /* z in ?[X,Y] mod Q(Y) in Kronecker form ((subst(lift(z), x, y^(2deg(z)-1)))
! 403: * Recover the "real" z, normalized */
! 404: GEN
! 405: from_Kronecker(GEN z, GEN pol)
! 406: {
! 407: long i,j,lx,l = lgef(z), N = ((lgef(pol)-3)<<1) + 1;
! 408: GEN a,x, t = cgetg(N,t_POL);
! 409: t[1] = pol[1] & VARNBITS;
! 410: lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
! 411: if (isonstack(pol)) pol = gcopy(pol);
! 412: for (i=2; i<lx+2; i++)
! 413: {
! 414: a = cgetg(3,t_POLMOD); x[i] = (long)a;
! 415: a[1] = (long)pol;
! 416: for (j=2; j<N; j++) t[j] = z[j];
! 417: z += (N-2);
! 418: a[2] = lres(normalizepol_i(t,N), pol);
! 419: }
! 420: a = cgetg(3,t_POLMOD); x[i] = (long)a;
! 421: a[1] = (long)pol;
! 422: N = (l-2) % (N-2) + 2;
! 423: for (j=2; j<N; j++) t[j] = z[j];
! 424: a[2] = lres(normalizepol_i(t,N), pol);
! 425: return normalizepol_i(x, i+1);
! 426: }
! 427:
! 428: /*******************************************************************/
! 429: /* */
! 430: /* MODULAR GCD */
! 431: /* */
! 432: /*******************************************************************/
! 433: static GEN
! 434: maxnorm(GEN p)
! 435: {
! 436: long i,n=lgef(p)-3,ltop=avma,lbot;
! 437: GEN x, m = gzero;
! 438:
! 439: p += 2;
! 440: for (i=0; i<n; i++)
! 441: {
! 442: x = (GEN)p[i];
! 443: if (absi_cmp(x,m) > 0) m = absi(x);
! 444: }
! 445: m = divii(m, absi((GEN)p[n])); lbot = avma;
! 446: return gerepile(ltop,lbot,addis(m,1));
! 447: }
! 448:
! 449: /* return x[0 .. dx] mod p as C-long in a malloc'ed array */
! 450: static GEN
! 451: Fp_to_pol_long(GEN x, long dx, long p, long *d)
! 452: {
! 453: long i, m;
! 454: GEN a;
! 455:
! 456: for (i=dx; i>=0; i--)
! 457: {
! 458: m = smodis((GEN)x[i],p);
! 459: if (m) break;
! 460: }
! 461: if (i < 0) { *d = -1; return NULL; }
! 462: a = (GEN) gpmalloc((i+1)*sizeof(long));
! 463: *d = i; a[i] = m;
! 464: for (i--; i>=0; i--) a[i] = smodis((GEN)x[i],p);
! 465: return a;
! 466: }
! 467:
! 468: /* idem as Fp_poldivres but working only on C-long types
! 469: * ASSUME pr != ONLY_DIVIDES (TODO ???)
! 470: */
! 471: static long *
! 472: Fp_poldivres_long(long *x,long *y,long p,long dx, long dy, long *dc, GEN *pr)
! 473: {
! 474: long dz,i,j,p1,*z,*c,inv;
! 475:
! 476: if (!dy) { *dc=-1; return NULL; }
! 477: dz=dx-dy;
! 478: if (dz<0)
! 479: {
! 480: if (pr)
! 481: {
! 482: c=(long *) gpmalloc((dx+1)*sizeof(long));
! 483: for (i=0; i<=dx; i++) c[i]=x[i];
! 484: *dc = dx;
! 485: if (pr == ONLY_REM) return c;
! 486: *pr = c;
! 487: }
! 488: return NULL;
! 489: }
! 490: z=(long *) gpmalloc((dz+1)*sizeof(long));
! 491: inv = y[dy];
! 492: if (inv!=1)
! 493: {
! 494: long av = avma;
! 495: GEN res = mpinvmod(stoi(y[dy]),stoi(p));
! 496: inv = itos(res); avma = av;
! 497: }
! 498:
! 499: z[dz]=(inv*x[dx])%p;
! 500: for (i=dx-1; i>=dy; --i)
! 501: {
! 502: p1=x[i];
! 503: for (j=i-dy+1; j<=i && j<=dz; j++)
! 504: {
! 505: p1 -= z[j]*y[i-j];
! 506: if (p1 & (HIGHBIT>>1)) p1=p1%p;
! 507: }
! 508: z[i-dy]=((p1%p)*inv)%p;
! 509: }
! 510: if (!pr) return z;
! 511:
! 512: c=(long *) gpmalloc(dy*sizeof(long));
! 513: for (i=0; i<dy; i++)
! 514: {
! 515: p1=z[0]*y[i];
! 516: for (j=1; j<=i && j<=dz; j++)
! 517: {
! 518: p1 += z[j]*y[i-j];
! 519: if (p1 & (HIGHBIT>>1)) p1=p1%p;
! 520: }
! 521: c[i]=(x[i]-p1)%p;
! 522: }
! 523:
! 524: i=dy-1; while (i>=0 && c[i]==0) i--;
! 525: *dc=i;
! 526: if (pr == ONLY_REM) { free(z); return c; }
! 527: *pr = c; return z;
! 528: }
! 529:
! 530: /* x and y in Z[X] */
! 531: GEN
! 532: Fp_poldivres(GEN x, GEN y, GEN p, GEN *pr)
! 533: {
! 534: long vx,dx,dy,dz,i,j,av0,av,tetpil,sx,lrem;
! 535: GEN z,p1,rem,lead;
! 536:
! 537: if (!signe(y)) err(talker,"division by zero in Fp_poldivres");
! 538: vx=varn(x); dy=lgef(y)-3; dx=lgef(x)-3;
! 539: if (dx < dy)
! 540: {
! 541: if (pr)
! 542: {
! 543: av0 = avma; x = Fp_pol_red(x, p);
! 544: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
! 545: if (pr == ONLY_REM) return x;
! 546: *pr = x;
! 547: }
! 548: return zeropol(vx);
! 549: }
! 550: lead = leading_term(y);
! 551: if (!dy) /* y is constant */
! 552: {
! 553: if (pr && pr != ONLY_DIVIDES)
! 554: {
! 555: if (pr == ONLY_REM) return zeropol(vx);
! 556: *pr = zeropol(vx);
! 557: }
! 558: if (gcmp1(lead)) return gcopy(x);
! 559: av0 = avma; x = gmul(x, mpinvmod(lead,p)); tetpil = avma;
! 560: return gerepile(av0,tetpil,Fp_pol_red(x,p));
! 561: }
! 562: av0 = avma; dz = dx-dy;
! 563: if (2*expi(p)+6<BITS_IN_LONG)
! 564: { /* assume ab != 0 mod p */
! 565: long *a, *b, *zz, da,db,dr, pp = p[2];
! 566: a = Fp_to_pol_long(x+2, dx, pp, &da);
! 567: b = Fp_to_pol_long(y+2, dy, pp, &db);
! 568: zz = Fp_poldivres_long(a,b,pp,da,db, &dr, pr);
! 569: if (pr == ONLY_REM) dz = dr;
! 570: else if (pr && pr != ONLY_DIVIDES)
! 571: {
! 572: rem = small_to_pol(*pr, dr+3, pp);
! 573: free(*pr); setvarn(rem, vx); *pr = rem;
! 574: }
! 575: z = small_to_pol(zz, dz+3, pp);
! 576: free(zz); free(a); free(b); setvarn(z, vx); return z;
! 577: }
! 578: lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p));
! 579: avma = av0;
! 580: z=cgetg(dz+3,t_POL);
! 581: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 582: x += 2; y += 2; z += 2;
! 583:
! 584: p1 = (GEN)x[dx]; av = avma;
! 585: z[dz] = lead? lpileupto(av, modii(mulii(p1,lead), p)): licopy(p1);
! 586: for (i=dx-1; i>=dy; i--)
! 587: {
! 588: av=avma; p1=(GEN)x[i];
! 589: for (j=i-dy+1; j<=i && j<=dz; j++)
! 590: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 591: if (lead) p1 = mulii(p1,lead);
! 592: tetpil=avma; z[i-dy] = lpile(av,tetpil,modii(p1, p));
! 593: }
! 594: if (!pr) { if (lead) gunclone(lead); return z-2; }
! 595:
! 596: rem = (GEN)avma; av = (long)new_chunk(dx+3);
! 597: for (sx=0; ; i--)
! 598: {
! 599: p1 = (GEN)x[i];
! 600: for (j=0; j<=i && j<=dz; j++)
! 601: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 602: tetpil=avma; p1 = modii(p1,p); if (signe(p1)) { sx = 1; break; }
! 603: if (!i) break;
! 604: avma=av;
! 605: }
! 606: if (pr == ONLY_DIVIDES)
! 607: {
! 608: if (lead) gunclone(lead);
! 609: if (sx) { avma=av0; return NULL; }
! 610: avma = (long)rem; return z-2;
! 611: }
! 612: lrem=i+3; rem -= lrem;
! 613: rem[0]=evaltyp(t_POL) | evallg(lrem);
! 614: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
! 615: p1 = gerepile((long)rem,tetpil,p1);
! 616: rem += 2; rem[i]=(long)p1;
! 617: for (i--; i>=0; i--)
! 618: {
! 619: av=avma; p1 = (GEN)x[i];
! 620: for (j=0; j<=i && j<=dz; j++)
! 621: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 622: tetpil=avma; rem[i]=lpile(av,tetpil, modii(p1,p));
! 623: }
! 624: rem -= 2;
! 625: if (lead) gunclone(lead);
! 626: if (!sx) normalizepol_i(rem, lrem);
! 627: if (pr == ONLY_REM) return gerepileupto(av0,rem);
! 628: *pr = rem; return z-2;
! 629: }
! 630:
! 631: static GEN
! 632: Fp_pol_gcd_long(GEN x, GEN y, GEN p)
! 633: {
! 634: long *a,*b,*c,da,db,dc, pp = (long)p[2];
! 635: GEN z;
! 636:
! 637: a = Fp_to_pol_long(x+2, lgef(x)-3, pp, &da);
! 638: if (!a) return Fp_pol_red(y,p);
! 639: b = Fp_to_pol_long(y+2, lgef(y)-3, pp, &db);
! 640: while (db>=0)
! 641: {
! 642: c = Fp_poldivres_long(a,b, pp, da,db,&dc, ONLY_REM);
! 643: free(a); a=b; da=db; b=c; db=dc;
! 644: }
! 645: if (b) free(b);
! 646: z = small_to_pol(a, da+3, pp);
! 647: setvarn(z, varn(x));
! 648: free(a); return z;
! 649: }
! 650:
! 651: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)) */
! 652: GEN
! 653: Fp_pol_gcd(GEN x, GEN y, GEN p)
! 654: {
! 655: GEN a,b,c;
! 656: long av0,av;
! 657:
! 658: if (2*expi(p)+6<BITS_IN_LONG) return Fp_pol_gcd_long(x,y,p);
! 659: av0=avma;
! 660: a = Fp_pol_red(x, p); av = avma;
! 661: b = Fp_pol_red(y, p);
! 662: while (signe(b))
! 663: {
! 664: av = avma; c = Fp_res(a,b,p); a=b; b=c;
! 665: }
! 666: avma = av; return gerepileupto(av0, a);
! 667: }
! 668:
! 669: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
! 670: * ux + vy = gcd (mod p)
! 671: */
! 672: GEN
! 673: Fp_pol_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
! 674: {
! 675: GEN a,b,q,r,u,v,d,d1,v1;
! 676: long ltop,lbot;
! 677:
! 678: #if 0 /* TODO */
! 679: if (2*expi(p)+6<BITS_IN_LONG) return Fp_pol_extgcd_long(x,y,p);
! 680: #endif
! 681: ltop=avma;
! 682: a = Fp_pol_red(x, p);
! 683: b = Fp_pol_red(y, p);
! 684: d = a; d1 = b; v = gzero; v1 = gun;
! 685: while (signe(d1))
! 686: {
! 687: q = Fp_poldivres(d,d1,p, &r);
! 688: v = gadd(v, gneg_i(gmul(q,v1)));
! 689: v = Fp_pol_red(v,p);
! 690: u=v; v=v1; v1=u;
! 691: u=r; d=d1; d1=u;
! 692: }
! 693: u = gadd(d, gneg_i(gmul(b,v)));
! 694: u = Fp_pol_red(u, p);
! 695: lbot = avma;
! 696: u = Fp_deuc(u,a,p);
! 697: d = gcopy(d);
! 698: v = gcopy(v);
! 699: {
! 700: GEN *gptr[3]; gptr[0] = &d; gptr[1] = &u; gptr[2] = &v;
! 701: gerepilemanysp(ltop,lbot,gptr,3);
! 702: }
! 703: *ptu = u; *ptv = v; return d;
! 704: }
! 705:
! 706: GEN chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1);
! 707:
! 708: /* a and b in Q[X] */
! 709: GEN
! 710: modulargcd(GEN a, GEN b)
! 711: {
! 712: GEN D,A,B,Cp,p,q,H,g,limit,ma,mb,p1;
! 713: long av=avma,avlim=stack_lim(av,1), m,n,nA,nB,av2,lbot,i;
! 714: long prime[]={evaltyp(t_INT)|m_evallg(3),evalsigne(1)|evallgefint(3),0};
! 715: byteptr d = diffptr;
! 716:
! 717: if (typ(a)!=t_POL || typ(b)!=t_POL) err(notpoler,"modulargcd");
! 718: if (!signe(a)) return gcopy(b);
! 719: if (!signe(b)) return gcopy(a);
! 720: A = content(a);
! 721: B = content(b); D = ggcd(A,B);
! 722: A = gcmp1(A)? a: gdiv(a,A); nA=lgef(A)-3;
! 723: B = gcmp1(B)? b: gdiv(b,B); nB=lgef(B)-3;
! 724: g = mppgcd((GEN)A[nA+2], (GEN)B[nB+2]);
! 725: av2=avma; n=1+min(nA,nB);
! 726: ma=maxnorm(A); mb=maxnorm(B);
! 727: if (cmpii(ma,mb) > 0) limit=mb; else limit=ma;
! 728: limit = shifti(mulii(limit,g), n+1);
! 729:
! 730: /* initial p could be 1 << (BITS_IN_LONG/2-6), but diffptr is nice */
! 731: prime[2] = 1021; d += 172; /* p = prime(172) = precprime(1<<10) */
! 732: p = prime; H = NULL;
! 733: for(;;)
! 734: {
! 735: do
! 736: {
! 737: if (*d) p[2] += *d++;
! 738: else p = nextprime(addis(p,1)); /* never used */
! 739: }
! 740: while (!signe(resii(g,p)));
! 741: Cp = Fp_pol_gcd(A,B,p);
! 742: m = lgef(Cp)-3;
! 743: if (m==0) return gerepileupto(av,gmul(D,polun[varn(A)]));
! 744: if (gcmp1(g))
! 745: Cp = normalize_mod_p(Cp, p);
! 746: else
! 747: { /* very rare */
! 748: p1 = mulii(g, mpinvmod((GEN)Cp[m+2],p));
! 749: Cp = Fp_pol_red(gmul(p1,Cp), p);
! 750: }
! 751: if (m<n) { q=icopy(p); H=Cp; limit=shifti(limit,m-n); n=m; }
! 752: else
! 753: if (m==n && H)
! 754: {
! 755: GEN q2 = mulii(q,p);
! 756: for (i=2; i<=n+2; i++)
! 757: H[i]=(long) chinois_int_coprime((GEN)H[i],(GEN)Cp[i],q,p,q2);
! 758: q = q2;
! 759: if (cmpii(limit,q) <= 0)
! 760: {
! 761: GEN limit2=shifti(limit,-1);
! 762: for (i=2; i<=n+2; i++)
! 763: {
! 764: p1 = (GEN)H[i];
! 765: if (cmpii(p1,limit2) > 0) H[i]=lsubii(p1,q);
! 766: }
! 767: p1 = content(H); if (!gcmp1(p1)) H = gdiv(H,p1);
! 768: if (!signe(gres(A,H)) && !signe(gres(B,H)))
! 769: {
! 770: lbot=avma;
! 771: return gerepile(av,lbot,gmul(D,H));
! 772: }
! 773: H = NULL; /* failed */
! 774: }
! 775: }
! 776: if (low_stack(avlim, stack_lim(av,1)))
! 777: {
! 778: GEN *gptr[4]; gptr[0]=&H; gptr[1]=&p; gptr[2]=&q; gptr[3]=&limit;
! 779: if (DEBUGMEM>1) err(warnmem,"modulargcd");
! 780: gerepilemany(av2,gptr,4);
! 781: }
! 782: }
! 783: }
! 784:
! 785: /* returns a polynomial in variable v, whose coeffs correspond to the digits
! 786: * of m (in base p)
! 787: */
! 788: GEN
! 789: stopoly(long m, long p, long v)
! 790: {
! 791: GEN y = cgetg(BITS_IN_LONG + 2, t_POL);
! 792: long l=2;
! 793:
! 794: do { y[l++] = lstoi(m%p); m=m/p; } while (m);
! 795: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
! 796: return y;
! 797: }
! 798:
! 799: GEN
! 800: stopoly_gen(GEN m, GEN p, long v)
! 801: {
! 802: GEN y = cgetg(bit_accuracy(lgefint(m))+2, t_POL);
! 803: long l=2;
! 804:
! 805: do { y[l++] = lmodii(m,p); m=divii(m,p); } while (signe(m));
! 806: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
! 807: return y;
! 808: }
! 809:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to polarit3.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath