Annotation of OpenXM_contrib/pari-2.2/src/basemath/polarit3.c, Revision 1.1
1.1 ! noro 1: /* $Id: polarit3.c,v 1.82 2001/10/01 17:19:06 bill 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: /** (third part) **/
! 20: /** **/
! 21: /***********************************************************************/
! 22: #include "pari.h"
! 23:
! 24: #define lswap(x,y) {long _z=x; x=y; y=_z;}
! 25: #define pswap(x,y) {GEN *_z=x; x=y; y=_z;}
! 26: #define swap(x,y) {GEN _z=x; x=y; y=_z;}
! 27: #define swapspec(x,y, nx,ny) {swap(x,y); lswap(nx,ny);}
! 28: /* 2p^2 < 2^BITS_IN_LONG */
! 29: #ifdef LONG_IS_64BIT
! 30: # define u_OK_ULONG(p) ((ulong)p <= 3037000493UL)
! 31: #else
! 32: # define u_OK_ULONG(p) ((ulong)p <= 46337UL)
! 33: #endif
! 34: #define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
! 35:
! 36: /*******************************************************************/
! 37: /* */
! 38: /* KARATSUBA (for polynomials) */
! 39: /* */
! 40: /*******************************************************************/
! 41:
! 42: #if 1 /* for tunings */
! 43: long SQR_LIMIT = 6;
! 44: long MUL_LIMIT = 10;
! 45: long u_SQR_LIMIT = 6;
! 46: long u_MUL_LIMIT = 10;
! 47:
! 48: void
! 49: setsqpol(long a) { SQR_LIMIT=a; u_SQR_LIMIT=a; }
! 50: void
! 51: setmulpol(long a) { MUL_LIMIT=a; u_MUL_LIMIT=a; }
! 52:
! 53: GEN
! 54: u_specpol(GEN x, long nx)
! 55: {
! 56: GEN z = cgetg(nx+2,t_POL);
! 57: long i;
! 58: for (i=0; i<nx; i++) z[i+2] = lstoi(x[i]);
! 59: z[1]=evalsigne(1)|evallgef(nx+2);
! 60: return z;
! 61: }
! 62:
! 63: GEN
! 64: specpol(GEN x, long nx)
! 65: {
! 66: GEN z = cgetg(nx+2,t_POL);
! 67: long i;
! 68: for (i=0; i<nx; i++) z[i+2] = x[i];
! 69: z[1]=evalsigne(1)|evallgef(nx+2);
! 70: return z;
! 71: }
! 72: #else
! 73: # define SQR_LIMIT 6
! 74: # define MUL_LIMIT 10
! 75: #endif
! 76:
! 77: #ifndef HIGHWORD
! 78: # define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG)
! 79: #endif
! 80:
! 81: /* multiplication in Fp[X], p small */
! 82:
! 83: static GEN
! 84: u_normalizepol(GEN x, long lx)
! 85: {
! 86: long i;
! 87: for (i=lx-1; i>1; i--)
! 88: if (x[i]) break;
! 89: setlgef(x,i+1);
! 90: setsigne(x, (i > 1)); return x;
! 91: }
! 92:
! 93: GEN
! 94: u_FpX_deriv(GEN z, ulong p)
! 95: {
! 96: long i,l = lgef(z)-1;
! 97: GEN x;
! 98: if (l < 2) l = 2;
! 99: x = cgetg(l,t_VECSMALL); x[1] = z[1]; z++;
! 100: if (HIGHWORD(l | p))
! 101: for (i=2; i<l; i++) x[i] = mulssmod(i-1, z[i], p);
! 102: else
! 103: for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
! 104: return u_normalizepol(x,l);
! 105: }
! 106:
! 107: static GEN
! 108: u_FpX_gcd_i(GEN a, GEN b, ulong p)
! 109: {
! 110: GEN c;
! 111: if (lgef(b) > lgef(a)) swap(a, b);
! 112: while (signe(b))
! 113: {
! 114: c = u_FpX_rem(a,b,p);
! 115: a = b; b = c;
! 116: }
! 117: return a;
! 118: }
! 119:
! 120: GEN
! 121: u_FpX_gcd(GEN a, GEN b, ulong p)
! 122: {
! 123: ulong av = avma;
! 124: return gerepileupto(av, u_FpX_gcd_i(a,b,p));
! 125: }
! 126:
! 127: int
! 128: u_FpX_is_squarefree(GEN z, ulong p)
! 129: {
! 130: ulong av = avma;
! 131: GEN d = u_FpX_gcd_i(z, u_FpX_deriv(z,p) , p);
! 132: avma = av; return (degpol(d) == 0);
! 133: }
! 134:
! 135: static GEN
! 136: u_FpX_addspec(GEN x, GEN y, long p, long lx, long ly)
! 137: {
! 138: long i,lz;
! 139: GEN z;
! 140:
! 141: if (ly>lx) swapspec(x,y, lx,ly);
! 142: lz = lx+2; z = cgetg(lz,t_VECSMALL) + 2;
! 143: for (i=0; i<ly; i++) z[i] = addssmod(x[i],y[i],p);
! 144: for ( ; i<lx; i++) z[i] = x[i];
! 145: z -= 2; z[1]=0; return u_normalizepol(z, lz);
! 146: }
! 147:
! 148: #ifdef INLINE
! 149: INLINE
! 150: #endif
! 151: ulong
! 152: u_FpX_mullimb(GEN x, GEN y, ulong p, long a, long b)
! 153: {
! 154: ulong p1 = 0;
! 155: long i;
! 156: for (i=a; i<b; i++)
! 157: if (y[i])
! 158: {
! 159: p1 += y[i] * x[-i];
! 160: if (p1 & HIGHBIT) p1 %= p;
! 161: }
! 162: return p1 % p;
! 163: }
! 164:
! 165: ulong
! 166: u_FpX_mullimb_safe(GEN x, GEN y, ulong p, long a, long b)
! 167: {
! 168: ulong p1 = 0;
! 169: long i;
! 170: for (i=a; i<b; i++)
! 171: if (y[i])
! 172: p1 = addssmod(p1, mulssmod(y[i],x[-i],p), p);
! 173: return p1;
! 174: }
! 175:
! 176: /* assume nx >= ny > 0 */
! 177: static GEN
! 178: s_FpX_mulspec(GEN x, GEN y, ulong p, long nx, long ny)
! 179: {
! 180: long i,lz,nz;
! 181: GEN z;
! 182:
! 183: lz = nx+ny+1; nz = lz-2;
! 184: z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
! 185: if (u_OK_ULONG(p))
! 186: {
! 187: for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb(x+i,y,p,0,i+1);
! 188: for ( ; i<nx; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,0,ny);
! 189: for ( ; i<nz; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,i-nx+1,ny);
! 190: }
! 191: else
! 192: {
! 193: for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,i+1);
! 194: for ( ; i<nx; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,ny);
! 195: for ( ; i<nz; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,i-nx+1,ny);
! 196: }
! 197: z -= 2; z[1]=0; return u_normalizepol(z, lz);
! 198: }
! 199:
! 200: /* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
! 201: GEN
! 202: u_FpX_addshift(GEN x, GEN y, ulong p, long d)
! 203: {
! 204: GEN xd,yd,zd = (GEN)avma;
! 205: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
! 206:
! 207: x += 2; y += 2; a = ny-d;
! 208: if (a <= 0)
! 209: {
! 210: lz = (a>nx)? ny+2: nx+d+2;
! 211: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
! 212: while (xd > x) *--zd = *--xd;
! 213: x = zd + a;
! 214: while (zd > x) *--zd = 0;
! 215: }
! 216: else
! 217: {
! 218: xd = new_chunk(d); yd = y+d;
! 219: x = u_FpX_addspec(x,yd,p, nx,a);
! 220: lz = (a>nx)? ny+2: lgef(x)+d;
! 221: x += 2; while (xd > x) *--zd = *--xd;
! 222: }
! 223: while (yd > y) *--zd = *--yd;
! 224: *--zd = evalsigne(1) | evallgef(lz);
! 225: *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
! 226: }
! 227:
! 228: #define u_zeropol(malloc) u_allocpol(-1,malloc)
! 229: #define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2)
! 230: #define u_FpX_sqr(x, p) u_FpX_Ksqr(x+2,p, lgef(x)-2)
! 231: #define u_FpX_add(x,y, p) u_FpX_addspec(x+2,y+2,p, lgef(x)-2,lgef(y)-2)
! 232:
! 233: static GEN
! 234: u_allocpol(long d, int malloc)
! 235: {
! 236: GEN z = malloc? (GEN)gpmalloc((d+3) * sizeof(long)): new_chunk(d+3);
! 237: z[0] = evaltyp(t_VECSMALL) | evallg(d+3);
! 238: z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0);
! 239: return z;
! 240: }
! 241:
! 242: static GEN
! 243: u_scalarpol(ulong x, int malloc)
! 244: {
! 245: GEN z = u_allocpol(0,malloc);
! 246: z[2] = x; return z;
! 247: }
! 248:
! 249: static GEN
! 250: u_FpX_neg_i(GEN x, ulong p)
! 251: {
! 252: long i, l = lgef(x);
! 253: for (i=2; i<l; i++)
! 254: if (x[i]) x[i] = p - x[i];
! 255: return x;
! 256: }
! 257:
! 258: /* shift polynomial + gerepile */
! 259: static GEN
! 260: u_FpX_shiftip(long av, GEN x, long v)
! 261: {
! 262: long i, lx = lgef(x), ly;
! 263: GEN y;
! 264: if (v <= 0 || !signe(x)) return gerepileupto(av, x);
! 265: avma = av; ly = lx + v;
! 266: x += lx; y = new_chunk(ly) + ly;
! 267: for (i = 2; i<lx; i++) *--y = *--x;
! 268: for (i = 0; i< v; i++) *--y = 0;
! 269: *--y = evalsigne(1) | evallgef(ly);
! 270: *--y = evaltyp(t_VECSMALL) | evallg(ly); return y;
! 271: }
! 272:
! 273: /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
! 274: * b+2 were sent instead. na, nb = number of terms of a, b.
! 275: * Only c, c0, c1, c2 are genuine GEN.
! 276: */
! 277: GEN
! 278: u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb)
! 279: {
! 280: GEN a0,c,c0;
! 281: long av,n0,n0a,i, v = 0;
! 282:
! 283: while (na && !a[0]) { a++; na--; v++; }
! 284: while (nb && !b[0]) { b++; nb--; v++; }
! 285: if (na < nb) swapspec(a,b, na,nb);
! 286: if (!nb) return u_zeropol(0);
! 287:
! 288: av = avma;
! 289: if (nb < u_MUL_LIMIT)
! 290: return u_FpX_shiftip(av,s_FpX_mulspec(a,b,p,na,nb), v);
! 291: i=(na>>1); n0=na-i; na=i;
! 292: a0=a+n0; n0a=n0;
! 293: while (n0a && !a[n0a-1]) n0a--;
! 294:
! 295: if (nb > n0)
! 296: {
! 297: GEN b0,c1,c2;
! 298: long n0b;
! 299:
! 300: nb -= n0; b0 = b+n0; n0b = n0;
! 301: while (n0b && !b[n0b-1]) n0b--;
! 302: c = u_FpX_Kmul(a,b,p,n0a,n0b);
! 303: c0 = u_FpX_Kmul(a0,b0,p,na,nb);
! 304:
! 305: c2 = u_FpX_addspec(a0,a,p,na,n0a);
! 306: c1 = u_FpX_addspec(b0,b,p,nb,n0b);
! 307:
! 308: c1 = u_FpX_mul(c1,c2,p);
! 309: c2 = u_FpX_add(c0,c,p);
! 310:
! 311: c2 = u_FpX_neg_i(c2,p);
! 312: c2 = u_FpX_add(c1,c2,p);
! 313: c0 = u_FpX_addshift(c0,c2 ,p, n0);
! 314: }
! 315: else
! 316: {
! 317: c = u_FpX_Kmul(a,b,p,n0a,nb);
! 318: c0 = u_FpX_Kmul(a0,b,p,na,nb);
! 319: }
! 320: c0 = u_FpX_addshift(c0,c,p,n0);
! 321: return u_FpX_shiftip(av,c0, v);
! 322: }
! 323:
! 324: GEN
! 325: u_FpX_sqrpol(GEN x, ulong p, long nx)
! 326: {
! 327: long i,j,l,lz,nz;
! 328: ulong p1;
! 329: GEN z;
! 330:
! 331: if (!nx) return u_zeropol(0);
! 332: lz = (nx << 1) + 1, nz = lz-2;
! 333: z = cgetg(lz,t_VECSMALL) + 2;
! 334: for (i=0; i<nx; i++)
! 335: {
! 336: p1 = 0; l = (i+1)>>1;
! 337: for (j=0; j<l; j++)
! 338: {
! 339: p1 += x[j] * x[i-j];
! 340: if (p1 & HIGHBIT) p1 %= p;
! 341: }
! 342: p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
! 343: if ((i&1) == 0)
! 344: p1 += x[i>>1] * x[i>>1];
! 345: z[i] = p1 % p;
! 346: }
! 347: for ( ; i<nz; i++)
! 348: {
! 349: p1 = 0; l = (i+1)>>1;
! 350: for (j=i-nx+1; j<l; j++)
! 351: {
! 352: p1 += x[j] * x[i-j];
! 353: if (p1 & HIGHBIT) p1 %= p;
! 354: }
! 355: p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
! 356: if ((i&1) == 0)
! 357: p1 += x[i>>1] * x[i>>1];
! 358: z[i] = p1 % p;
! 359: }
! 360: z -= 2; z[1]=0; return u_normalizepol(z,lz);
! 361: }
! 362:
! 363: static GEN
! 364: u_Fp_gmul2_1(GEN x, ulong p)
! 365: {
! 366: long i, l = lgef(x);
! 367: GEN z = cgetg(l, t_VECSMALL);
! 368: for (i=2; i<l; i++)
! 369: {
! 370: ulong p1 = x[i]<<1;
! 371: if (p1 > p) p1 -= p;
! 372: z[i] = p1;
! 373: }
! 374: z[1] = x[1]; return z;
! 375: }
! 376:
! 377: GEN
! 378: u_FpX_Ksqr(GEN a, ulong p, long na)
! 379: {
! 380: GEN a0,c,c0,c1;
! 381: long av,n0,n0a,i, v = 0;
! 382:
! 383: while (na && !a[0]) { a++; na--; v += 2; }
! 384: av = avma;
! 385: if (na < u_SQR_LIMIT) return u_FpX_shiftip(av, u_FpX_sqrpol(a,p,na), v);
! 386: i=(na>>1); n0=na-i; na=i;
! 387: a0=a+n0; n0a=n0;
! 388: while (n0a && !a[n0a-1]) n0a--;
! 389:
! 390: c = u_FpX_Ksqr(a,p,n0a);
! 391: c0= u_FpX_Ksqr(a0,p,na);
! 392: if (p == 2) n0 *= 2;
! 393: else
! 394: {
! 395: c1 = u_Fp_gmul2_1(u_FpX_Kmul(a0,a,p, na,n0a), p);
! 396: c0 = u_FpX_addshift(c0,c1,p,n0);
! 397: }
! 398: c0 = u_FpX_addshift(c0,c,p,n0);
! 399: return u_FpX_shiftip(av,c0,v);
! 400: }
! 401:
! 402: /* generic multiplication */
! 403:
! 404: static GEN
! 405: addpol(GEN x, GEN y, long lx, long ly)
! 406: {
! 407: long i,lz;
! 408: GEN z;
! 409:
! 410: if (ly>lx) swapspec(x,y, lx,ly);
! 411: lz = lx+2; z = cgetg(lz,t_POL) + 2;
! 412: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 413: for ( ; i<lx; i++) z[i]=x[i];
! 414: z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 415: }
! 416:
! 417: static GEN
! 418: addpolcopy(GEN x, GEN y, long lx, long ly)
! 419: {
! 420: long i,lz;
! 421: GEN z;
! 422:
! 423: if (ly>lx) swapspec(x,y, lx,ly);
! 424: lz = lx+2; z = cgetg(lz,t_POL) + 2;
! 425: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 426: for ( ; i<lx; i++) z[i]=lcopy((GEN)x[i]);
! 427: z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 428: }
! 429:
! 430: #ifdef INLINE
! 431: INLINE
! 432: #endif
! 433: GEN
! 434: mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b)
! 435: {
! 436: GEN p1 = NULL;
! 437: long i,av = avma;
! 438: for (i=a; i<b; i++)
! 439: if (ynonzero[i])
! 440: {
! 441: GEN p2 = gmul((GEN)y[i],(GEN)x[-i]);
! 442: p1 = p1 ? gadd(p1, p2): p2;
! 443: }
! 444: return p1 ? gerepileupto(av, p1): gzero;
! 445: }
! 446:
! 447: /* assume nx >= ny > 0 */
! 448: static GEN
! 449: mulpol(GEN x, GEN y, long nx, long ny)
! 450: {
! 451: long i,lz,nz;
! 452: GEN z;
! 453: char *p1;
! 454:
! 455: lz = nx+ny+1; nz = lz-2;
! 456: z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
! 457: p1 = gpmalloc(ny);
! 458: for (i=0; i<ny; i++)
! 459: {
! 460: p1[i] = !isexactzero((GEN)y[i]);
! 461: z[i] = (long)mulpol_limb(x+i,y,p1,0,i+1);
! 462: }
! 463: for ( ; i<nx; i++) z[i] = (long)mulpol_limb(x+i,y,p1,0,ny);
! 464: for ( ; i<nz; i++) z[i] = (long)mulpol_limb(x+i,y,p1,i-nx+1,ny);
! 465: free(p1); z -= 2; z[1]=0; return normalizepol_i(z, lz);
! 466: }
! 467:
! 468: /* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
! 469: GEN
! 470: addshiftw(GEN x, GEN y, long d)
! 471: {
! 472: GEN xd,yd,zd = (GEN)avma;
! 473: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
! 474:
! 475: x += 2; y += 2; a = ny-d;
! 476: if (a <= 0)
! 477: {
! 478: lz = (a>nx)? ny+2: nx+d+2;
! 479: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
! 480: while (xd > x) *--zd = *--xd;
! 481: x = zd + a;
! 482: while (zd > x) *--zd = zero;
! 483: }
! 484: else
! 485: {
! 486: xd = new_chunk(d); yd = y+d;
! 487: x = addpol(x,yd, nx,a);
! 488: lz = (a>nx)? ny+2: lgef(x)+d;
! 489: x += 2; while (xd > x) *--zd = *--xd;
! 490: }
! 491: while (yd > y) *--zd = *--yd;
! 492: *--zd = evalsigne(1) | evallgef(lz);
! 493: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
! 494: }
! 495:
! 496: GEN
! 497: addshiftpol(GEN x, GEN y, long d)
! 498: {
! 499: long v = varn(x);
! 500: if (!signe(x)) return y;
! 501: x = addshiftw(x,y,d);
! 502: setvarn(x,v); return x;
! 503: }
! 504:
! 505: /* as above, producing a clean malloc */
! 506: static GEN
! 507: addshiftwcopy(GEN x, GEN y, long d)
! 508: {
! 509: GEN xd,yd,zd = (GEN)avma;
! 510: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
! 511:
! 512: x += 2; y += 2; a = ny-d;
! 513: if (a <= 0)
! 514: {
! 515: lz = nx+d+2;
! 516: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
! 517: while (xd > x) *--zd = lcopy((GEN)*--xd);
! 518: x = zd + a;
! 519: while (zd > x) *--zd = zero;
! 520: }
! 521: else
! 522: {
! 523: xd = new_chunk(d); yd = y+d;
! 524: x = addpolcopy(x,yd, nx,a);
! 525: lz = (a>nx)? ny+2: lgef(x)+d;
! 526: x += 2; while (xd > x) *--zd = *--xd;
! 527: }
! 528: while (yd > y) *--zd = lcopy((GEN)*--yd);
! 529: *--zd = evalsigne(1) | evallgef(lz);
! 530: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
! 531: }
! 532:
! 533: /* shift polynomial in place. assume v free cells have been left before x */
! 534: static GEN
! 535: shiftpol_ip(GEN x, long v)
! 536: {
! 537: long i, lx;
! 538: GEN y;
! 539: if (v <= 0 || !signe(x)) return x;
! 540: lx = lgef(x);
! 541: x += 2; y = x + v;
! 542: for (i = lx-3; i>=0; i--) y[i] = x[i];
! 543: for (i = 0 ; i< v; i++) x[i] = zero;
! 544: lx += v;
! 545: *--x = evalsigne(1) | evallgef(lx);
! 546: *--x = evaltyp(t_POL) | evallg(lx); return x;
! 547: }
! 548:
! 549: /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
! 550: * b+2 were sent instead. na, nb = number of terms of a, b.
! 551: * Only c, c0, c1, c2 are genuine GEN.
! 552: */
! 553: GEN
! 554: quickmul(GEN a, GEN b, long na, long nb)
! 555: {
! 556: GEN a0,c,c0;
! 557: long av,n0,n0a,i, v = 0;
! 558:
! 559: while (na && isexactzero((GEN)a[0])) { a++; na--; v++; }
! 560: while (nb && isexactzero((GEN)b[0])) { b++; nb--; v++; }
! 561: if (na < nb) swapspec(a,b, na,nb);
! 562: if (!nb) return zeropol(0);
! 563:
! 564: if (v) (void)cgetg(v,t_STR); /* v gerepile-safe cells for shiftpol_ip */
! 565: if (nb < MUL_LIMIT)
! 566: return shiftpol_ip(mulpol(a,b,na,nb), v);
! 567: i=(na>>1); n0=na-i; na=i;
! 568: av=avma; a0=a+n0; n0a=n0;
! 569: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
! 570:
! 571: if (nb > n0)
! 572: {
! 573: GEN b0,c1,c2;
! 574: long n0b;
! 575:
! 576: nb -= n0; b0 = b+n0; n0b = n0;
! 577: while (n0b && isexactzero((GEN)b[n0b-1])) n0b--;
! 578: c = quickmul(a,b,n0a,n0b);
! 579: c0 = quickmul(a0,b0, na,nb);
! 580:
! 581: c2 = addpol(a0,a, na,n0a);
! 582: c1 = addpol(b0,b, nb,n0b);
! 583:
! 584: c1 = quickmul(c1+2,c2+2, lgef(c1)-2,lgef(c2)-2);
! 585: c2 = gneg_i(gadd(c0,c));
! 586: c0 = addshiftw(c0, gadd(c1,c2), n0);
! 587: }
! 588: else
! 589: {
! 590: c = quickmul(a,b,n0a,nb);
! 591: c0 = quickmul(a0,b,na,nb);
! 592: }
! 593: c0 = addshiftwcopy(c0,c,n0);
! 594: return shiftpol_ip(gerepileupto(av,c0), v);
! 595: }
! 596:
! 597: GEN
! 598: sqrpol(GEN x, long nx)
! 599: {
! 600: long av,i,j,l,lz,nz;
! 601: GEN p1,z;
! 602: char *p2;
! 603:
! 604: if (!nx) return zeropol(0);
! 605: lz = (nx << 1) + 1, nz = lz-2;
! 606: z = cgetg(lz,t_POL) + 2;
! 607: p2 = gpmalloc(nx);
! 608: for (i=0; i<nx; i++)
! 609: {
! 610: p2[i] = !isexactzero((GEN)x[i]);
! 611: p1=gzero; av=avma; l=(i+1)>>1;
! 612: for (j=0; j<l; j++)
! 613: if (p2[j] && p2[i-j])
! 614: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
! 615: p1 = gshift(p1,1);
! 616: if ((i&1) == 0 && p2[i>>1])
! 617: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
! 618: z[i] = lpileupto(av,p1);
! 619: }
! 620: for ( ; i<nz; i++)
! 621: {
! 622: p1=gzero; av=avma; l=(i+1)>>1;
! 623: for (j=i-nx+1; j<l; j++)
! 624: if (p2[j] && p2[i-j])
! 625: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
! 626: p1 = gshift(p1,1);
! 627: if ((i&1) == 0 && p2[i>>1])
! 628: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
! 629: z[i] = lpileupto(av,p1);
! 630: }
! 631: free(p2); z -= 2; z[1]=0; return normalizepol_i(z,lz);
! 632: }
! 633:
! 634: GEN
! 635: quicksqr(GEN a, long na)
! 636: {
! 637: GEN a0,c,c0,c1;
! 638: long av,n0,n0a,i, v = 0;
! 639:
! 640: while (na && isexactzero((GEN)a[0])) { a++; na--; v += 2; }
! 641: if (v) (void)new_chunk(v);
! 642: if (na<SQR_LIMIT) return shiftpol_ip(sqrpol(a,na), v);
! 643: i=(na>>1); n0=na-i; na=i;
! 644: av=avma; a0=a+n0; n0a=n0;
! 645: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
! 646:
! 647: c = quicksqr(a,n0a);
! 648: c0 = quicksqr(a0,na);
! 649: c1 = gmul2n(quickmul(a0,a, na,n0a), 1);
! 650: c0 = addshiftw(c0,c1, n0);
! 651: c0 = addshiftwcopy(c0,c,n0);
! 652: return shiftpol_ip(gerepileupto(av,c0), v);
! 653: }
! 654: /*****************************************
! 655: * Arithmetic in Z/pZ[X] *
! 656: *****************************************/
! 657:
! 658: /*********************************************************************
! 659: This functions supposes polynomials already reduced.
! 660: There are clean and memory efficient.
! 661: **********************************************************************/
! 662:
! 663: GEN
! 664: FpX_center(GEN T,GEN mod)
! 665: {/*OK centermod exists, but is not so clean*/
! 666: ulong av;
! 667: long i, l=lg(T);
! 668: GEN P,mod2;
! 669: P=cgetg(l,t_POL);
! 670: P[1]=T[1];
! 671: av=avma;
! 672: mod2=gclone(shifti(mod,-1));/*clone*/
! 673: avma=av;
! 674: for(i=2;i<l;i++)
! 675: P[i]=cmpii((GEN)T[i],mod2)<0?licopy((GEN)T[i]):lsubii((GEN)T[i],mod);
! 676: gunclone(mod2);/*unclone*/
! 677: return P;
! 678: }
! 679:
! 680: GEN
! 681: FpX_neg(GEN x,GEN p)
! 682: {
! 683: long i,d=lgef(x);
! 684: GEN y;
! 685: y=cgetg(d,t_POL); y[1]=x[1];
! 686: for(i=2;i<d;i++)
! 687: if (signe(x[i])) y[i]=lsubii(p,(GEN)x[i]);
! 688: else y[i]=zero;
! 689: return y;
! 690: }
! 691: /**********************************************************************
! 692: Unclean functions, do not garbage collect.
! 693: This is a feature: The malloc is corrupted only by the call to FpX_red
! 694: so garbage collecting so often is not desirable.
! 695: FpX_red can sometime be avoided by passing NULL for p.
! 696: In this case the function is usually clean (see below for detail)
! 697: Added to help not using POLMOD of INTMOD which are deadly slow.
! 698: gerepileupto of the result is legible. Bill.
! 699: I don't like C++. I am wrong.
! 700: **********************************************************************/
! 701: /*
! 702: *If p is NULL no reduction is performed and the function is clean.
! 703: * for FpX_add,FpX_mul,FpX_sqr,FpX_Fp_mul
! 704: */
! 705: GEN
! 706: FpX_add(GEN x,GEN y,GEN p)
! 707: {
! 708: long lx,ly,i;
! 709: GEN z;
! 710: lx = lgef(x); ly = lgef(y); if (lx < ly) swapspec(x,y, lx,ly);
! 711: z = cgetg(lx,t_POL); z[1] = x[1];
! 712: for (i=2; i<ly; i++) z[i]=laddii((GEN)x[i],(GEN)y[i]);
! 713: for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]);
! 714: (void)normalizepol_i(z, lx);
! 715: if (lgef(z) == 2) { avma = (long)(z + lx); z = zeropol(varn(x)); }
! 716: if (p) z= FpX_red(z, p);
! 717: return z;
! 718: }
! 719: GEN
! 720: FpX_sub(GEN x,GEN y,GEN p)
! 721: {
! 722: long lx,ly,i,lz;
! 723: GEN z;
! 724: lx = lgef(x); ly = lgef(y);
! 725: lz=max(lx,ly);
! 726: z = cgetg(lz,t_POL);
! 727: if (lx >= ly)
! 728: {
! 729: z[1] = x[1];
! 730: for (i=2; i<ly; i++) z[i]=lsubii((GEN)x[i],(GEN)y[i]);
! 731: for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]);
! 732: (void)normalizepol_i(z, lz);
! 733: }
! 734: else
! 735: {
! 736: z[1] = y[1];
! 737: for (i=2; i<lx; i++) z[i]=lsubii((GEN)x[i],(GEN)y[i]);
! 738: for ( ; i<ly; i++) z[i]=lnegi((GEN)y[i]);
! 739: /*polynomial is always normalized*/
! 740: }
! 741: if (lgef(z) == 2) { avma = (long)(z + lz); z = zeropol(varn(x)); }
! 742: if (p) z= FpX_red(z, p);
! 743: return z;
! 744: }
! 745: GEN
! 746: FpX_mul(GEN x,GEN y,GEN p)
! 747: {
! 748: GEN z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2);
! 749: setvarn(z,varn(x));
! 750: if (!p) return z;
! 751: return FpX_red(z, p);
! 752: }
! 753: GEN
! 754: FpX_sqr(GEN x,GEN p)
! 755: {
! 756: GEN z = quicksqr(x+2, lgef(x)-2);
! 757: setvarn(z,varn(x));
! 758: if (!p) return z;
! 759: return FpX_red(z, p);
! 760: }
! 761: #define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),(0),NULL)
! 762:
! 763: /* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
! 764: static GEN
! 765: u_FpXQ_mul(GEN y,GEN x,GEN pol,ulong p)
! 766: {
! 767: GEN z = u_FpX_mul(y,x,p);
! 768: return u_FpX_rem(z,pol,p);
! 769: }
! 770: /* Square of y in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
! 771: static GEN
! 772: u_FpXQ_sqr(GEN y,GEN pol,ulong p)
! 773: {
! 774: GEN z = u_FpX_sqr(y,p);
! 775: return u_FpX_rem(z,pol,p);
! 776: }
! 777:
! 778: /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
! 779: * return lift(1 / (x mod (p,pol))) */
! 780: GEN
! 781: FpXQ_invsafe(GEN x, GEN T, GEN p)
! 782: {
! 783: GEN z, U, V;
! 784:
! 785: if (typ(x) != t_POL) return mpinvmod(x, p);
! 786: z = FpX_extgcd(x, T, p, &U, &V);
! 787: if (lgef(z) != 3) return NULL;
! 788: z = mpinvmod((GEN)z[2], p);
! 789: return FpX_Fp_mul(U, z, p);
! 790: }
! 791:
! 792: /* Product of y and x in Z/pZ[X]/(T)
! 793: * return lift(lift(Mod(x*y*Mod(1,p),T*Mod(1,p)))) */
! 794: GEN
! 795: FpXQ_mul(GEN y,GEN x,GEN T,GEN p)
! 796: {
! 797: GEN z;
! 798: if (typ(x) == t_INT)
! 799: {
! 800: if (typ(y) == t_INT) return modii(mulii(x,y), p);
! 801: return FpX_Fp_mul(y, x, p);
! 802: }
! 803: if (typ(y) == t_INT) return FpX_Fp_mul(x, y, p);
! 804: z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2); setvarn(z,varn(y));
! 805: z = FpX_red(z, p); return FpX_res(z,T, p);
! 806: }
! 807:
! 808: /* Square of y in Z/pZ[X]/(pol)
! 809: * return lift(lift(Mod(y^2*Mod(1,p),pol*Mod(1,p)))) */
! 810: GEN
! 811: FpXQ_sqr(GEN y,GEN pol,GEN p)
! 812: {
! 813: GEN z = quicksqr(y+2,lgef(y)-2); setvarn(z,varn(y));
! 814: z = FpX_red(z, p); return FpX_res(z,pol, p);
! 815: }
! 816: /*Modify y[2].
! 817: *No reduction if p is NULL
! 818: */
! 819: GEN
! 820: FpX_Fp_add(GEN y,GEN x,GEN p)
! 821: {
! 822: if (!signe(x)) return y;
! 823: if (!signe(y))
! 824: return scalarpol(x,varn(y));
! 825: y[2]=laddii((GEN)y[2],x);
! 826: if (p) y[2]=lmodii((GEN)y[2],p);
! 827: if (!signe(y[2]) && degpol(y) == 0) return zeropol(varn(y));
! 828: return y;
! 829: }
! 830: GEN
! 831: ZX_s_add(GEN y,long x)
! 832: {
! 833: if (!x) return y;
! 834: if (!signe(y))
! 835: return scalarpol(stoi(x),varn(y));
! 836: y[2] = laddis((GEN)y[2],x);
! 837: if (!signe(y[2]) && degpol(y) == 0) return zeropol(varn(y));
! 838: return y;
! 839: }
! 840: /* y is a polynomial in ZZ[X] and x an integer.
! 841: * If p is NULL, no reduction is perfomed and return x*y
! 842: *
! 843: * else the result is lift(y*x*Mod(1,p))
! 844: */
! 845: GEN
! 846: FpX_Fp_mul(GEN y,GEN x,GEN p)
! 847: {
! 848: GEN z;
! 849: int i;
! 850: if (!signe(x))
! 851: return zeropol(varn(y));
! 852: z=cgetg(lg(y),t_POL);
! 853: z[1]=y[1];
! 854: for(i=2;i<lgef(y);i++)
! 855: z[i]=lmulii((GEN)y[i],x);
! 856: if(!p) return z;
! 857: return FpX_red(z,p);
! 858: }
! 859: /*****************************************************************
! 860: * End of unclean functions. *
! 861: * *
! 862: *****************************************************************/
! 863: /*****************************************************************
! 864: Clean and with no reduced hypothesis. Beware that some operations
! 865: will be much slower with big unreduced coefficient
! 866: *****************************************************************/
! 867: /* Inverse of x in Z[X] / (p,T)
! 868: * return lift(lift(Mod(x*Mod(1,p), T*Mod(1,p))^-1)); */
! 869: GEN
! 870: FpXQ_inv(GEN x,GEN T,GEN p)
! 871: {
! 872: ulong av;
! 873: GEN U;
! 874:
! 875: if (!T) return mpinvmod(x,p);
! 876: av = avma;
! 877: U = FpXQ_invsafe(x, T, p);
! 878: if (!U) err(talker,"non invertible polynomial in FpXQ_inv");
! 879: return gerepileupto(av, U);
! 880: }
! 881:
! 882: GEN
! 883: FpXV_FpV_dotproduct(GEN V, GEN W, GEN p)
! 884: {
! 885: long ltop=avma;
! 886: long i;
! 887: GEN z = FpX_Fp_mul((GEN)V[1],(GEN)W[1],NULL);
! 888: for(i=2;i<lg(V);i++)
! 889: z=FpX_add(z,FpX_Fp_mul((GEN)V[i],(GEN)W[i],NULL),NULL);
! 890: return gerepileupto(ltop,FpX_red(z,p));
! 891: }
! 892:
! 893: /* generates the list of powers of x of degree 0,1,2,...,l*/
! 894: GEN
! 895: FpXQ_powers(GEN x, long l, GEN T, GEN p)
! 896: {
! 897: GEN V=cgetg(l+2,t_VEC);
! 898: long i;
! 899: V[1] = (long) scalarpol(gun,varn(T));
! 900: if (l==0) return V;
! 901: V[2] = lcopy(x);
! 902: if (l==1) return V;
! 903: V[3] = (long) FpXQ_sqr(x,T,p);
! 904: for(i=4;i<l+2;i++)
! 905: V[i] = (long) FpXQ_mul((GEN) V[i-1],x,T,p);
! 906: #if 0
! 907: TODO: Karim proposes to use squaring:
! 908: V[i] = (long) ((i&1)?FpXQ_sqr((GEN) V[(i+1)>>1],T,p)
! 909: :FpXQ_mul((GEN) V[i-1],x,T,p));
! 910: Please profile it.
! 911: #endif
! 912: return V;
! 913: }
! 914: #if 0
! 915: static long brent_kung_nbmul(long d, long n, long p)
! 916: {
! 917: return p+n*((d+p-1)/p);
! 918: }
! 919: TODO: This the the optimal parameter for the purpose of reducing
! 920: multiplications, but profiling need to be done to ensure it is not slower
! 921: than the other option in practice
! 922: /*Return optimal parameter l for the evaluation of n polynomials of degree d*/
! 923: long brent_kung_optpow(long d, long n)
! 924: {
! 925: double r;
! 926: long f,c,pr;
! 927: if (n>=d ) return d;
! 928: pr=n*d;
! 929: if (pr<=1) return 1;
! 930: r=d/sqrt(pr);
! 931: c=(long)ceil(d/ceil(r));
! 932: f=(long)floor(d/floor(r));
! 933: return (brent_kung_nbmul(d, n, c) <= brent_kung_nbmul(d, n, f))?c:f;
! 934: }
! 935: #endif
! 936: /*Return optimal parameter l for the evaluation of n polynomials of degree d*/
! 937: long brent_kung_optpow(long d, long n)
! 938: {
! 939: double r;
! 940: long l,pr;
! 941: if (n>=d ) return d;
! 942: pr=n*d;
! 943: if (pr<=1) return 1;
! 944: r=d/sqrt(pr);
! 945: l=(long)r;
! 946: return (d+l-1)/l;
! 947: }
! 948:
! 949: /*Close to FpXV_FpV_dotproduct*/
! 950:
! 951: static GEN
! 952: spec_compo_powers(GEN P, GEN V, long a, long n)
! 953: {
! 954: long i;
! 955: GEN z;
! 956: z = scalarpol((GEN)P[2+a],varn(P));
! 957: for(i=1;i<=n;i++)
! 958: z=FpX_add(z,FpX_Fp_mul((GEN)V[i+1],(GEN)P[2+a+i],NULL),NULL);
! 959: return z;
! 960: }
! 961: /*Try to implement algorithm in Brent & Kung (Fast algorithms for
! 962: *manipulating formal power series, JACM 25:581-595, 1978)
! 963:
! 964: V must be as output by FpXQ_powers.
! 965: For optimal performance, l (of FpXQ_powers) must be as output by
! 966: brent_kung_optpow
! 967: */
! 968:
! 969: GEN
! 970: FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p)
! 971: {
! 972: ulong ltop=avma;
! 973: long l=lg(V)-1;
! 974: GEN z,u;
! 975: long d=degpol(P),cnt=0;
! 976: if (d==-1) return zeropol(varn(T));
! 977: if (d<l)
! 978: {
! 979: z=spec_compo_powers(P,V,0,d);
! 980: return gerepileupto(ltop,FpX_red(z,p));
! 981: }
! 982: if (l<=1)
! 983: err(talker,"powers is only [] or [1] in FpX_FpXQV_compo");
! 984: z=spec_compo_powers(P,V,d-l+1,l-1);
! 985: d-=l;
! 986: while(d>=l-1)
! 987: {
! 988: u=spec_compo_powers(P,V,d-l+2,l-2);
! 989: z=FpX_add(u,FpXQ_mul(z,(GEN)V[l],T,p),NULL);
! 990: d-=l-1;
! 991: cnt++;
! 992: }
! 993: u=spec_compo_powers(P,V,0,d);
! 994: z=FpX_add(u,FpXQ_mul(z,(GEN)V[d+2],T,p),NULL);
! 995: cnt++;
! 996: if (DEBUGLEVEL>=8) fprintferr("FpX_FpXQV_compo: %d FpXQ_mul [%d]\n",cnt,l-1);
! 997: return gerepileupto(ltop,FpX_red(z,p));
! 998: }
! 999:
! 1000: /* T in Z[X] and x in Z/pZ[X]/(pol)
! 1001: * return lift(lift(subst(T,variable(T),Mod(x*Mod(1,p),pol*Mod(1,p)))));
! 1002: */
! 1003: GEN
! 1004: FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p)
! 1005: {
! 1006: ulong ltop=avma;
! 1007: GEN z;
! 1008: long d=degpol(T),rtd;
! 1009: if (!signe(T)) return zeropol(varn(T));
! 1010: rtd = (long) sqrt((double)d);
! 1011: z = FpX_FpXQV_compo(T,FpXQ_powers(x,rtd,pol,p),pol,p);
! 1012: return gerepileupto(ltop,z);
! 1013: }
! 1014:
! 1015: /* Evaluation in Fp
! 1016: * x in Z[X] and y in Z return x(y) mod p
! 1017: *
! 1018: * If p is very large (several longs) and x has small coefficients(<<p),
! 1019: * then Brent & Kung algorithm is faster. */
! 1020: GEN
! 1021: FpX_eval(GEN x,GEN y,GEN p)
! 1022: {
! 1023: ulong av;
! 1024: GEN p1,r,res;
! 1025: long i,j;
! 1026: i=lgef(x)-1;
! 1027: if (i<=2)
! 1028: return (i==2)? modii((GEN)x[2],p): gzero;
! 1029: res=cgetg(lgefint(p),t_INT);
! 1030: av=avma; p1=(GEN)x[i];
! 1031: /* specific attention to sparse polynomials (see poleval)*/
! 1032: /*You've guess it! It's a copy-paste(tm)*/
! 1033: for (i--; i>=2; i=j-1)
! 1034: {
! 1035: for (j=i; !signe((GEN)x[j]); j--)
! 1036: if (j==2)
! 1037: {
! 1038: if (i!=j) y = powmodulo(y,stoi(i-j+1),p);
! 1039: p1=mulii(p1,y);
! 1040: goto fppoleval;/*sorry break(2) no implemented*/
! 1041: }
! 1042: r = (i==j)? y: powmodulo(y,stoi(i-j+1),p);
! 1043: p1 = modii(addii(mulii(p1,r), (GEN)x[j]),p);
! 1044: }
! 1045: fppoleval:
! 1046: modiiz(p1,p,res);
! 1047: avma=av;
! 1048: return res;
! 1049: }
! 1050: /* Tz=Tx*Ty where Tx and Ty coprime
! 1051: * return lift(chinese(Mod(x*Mod(1,p),Tx*Mod(1,p)),Mod(y*Mod(1,p),Ty*Mod(1,p))))
! 1052: * if Tz is NULL it is computed
! 1053: * As we do not return it, and the caller will frequently need it,
! 1054: * it must compute it and pass it.
! 1055: */
! 1056: GEN
! 1057: FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p)
! 1058: {
! 1059: long av = avma;
! 1060: GEN ax,p1;
! 1061: ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p);
! 1062: p1=FpX_mul(ax, FpX_sub(y,x,p),p);
! 1063: p1 = FpX_add(x,p1,p);
! 1064: if (!Tz) Tz=FpX_mul(Tx,Ty,p);
! 1065: p1 = FpX_res(p1,Tz,p);
! 1066: return gerepileupto(av,p1);
! 1067: }
! 1068:
! 1069: /* assume n > 0 */
! 1070: GEN
! 1071: u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p)
! 1072: {
! 1073: ulong av = avma;
! 1074: GEN p1 = n+2, y = x;
! 1075: long m,i,j;
! 1076: m = *p1;
! 1077: j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j;
! 1078: for (i=lgefint(n)-2;;)
! 1079: {
! 1080: for (; j; m<<=1,j--)
! 1081: {
! 1082: y = u_FpXQ_sqr(y,pol,p);
! 1083: if (m<0)
! 1084: y = u_FpXQ_mul(y,x,pol,p);
! 1085: }
! 1086: if (--i == 0) break;
! 1087: m = *++p1, j = BITS_IN_LONG;
! 1088: }
! 1089: return gerepileupto(av, y);
! 1090: }
! 1091:
! 1092: /* x,pol in Z[X], p in Z, n in Z, compute lift(x^n mod (p, pol)) */
! 1093: GEN
! 1094: FpXQ_pow(GEN x, GEN n, GEN pol, GEN p)
! 1095: {
! 1096: ulong av, ltop = avma, lim=stack_lim(avma,1);
! 1097: long m,i,j, vx = varn(x);
! 1098: GEN p1 = n+2, y;
! 1099: if (!signe(n)) return polun[vx];
! 1100: if (signe(n) < 0)
! 1101: {
! 1102: x=FpXQ_inv(x,pol,p);
! 1103: if (is_pm1(n)) return x; /* n = -1*/
! 1104: }
! 1105: else
! 1106: if (is_pm1(n)) return gcopy(x); /* n = 1 */
! 1107: if (OK_ULONG(p))
! 1108: {
! 1109: ulong pp = p[2];
! 1110: pol = u_Fp_FpX(pol,0, pp);
! 1111: x = u_Fp_FpX(x,0, pp);
! 1112: y = u_FpXQ_pow(x, n, pol, pp);
! 1113: y = small_to_pol(y,vx);
! 1114: }
! 1115: else
! 1116: {
! 1117: av = avma;
! 1118: m = *p1; y = x;
! 1119: j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j;
! 1120: for (i=lgefint(n)-2;;)
! 1121: {
! 1122: for (; j; m<<=1,j--)
! 1123: {
! 1124: y = FpXQ_sqr(y,pol,p);
! 1125: if (low_stack(lim, stack_lim(av,1)))
! 1126: {
! 1127: if(DEBUGMEM>1) err(warnmem,"[1]: FpXQ_pow");
! 1128: y = gerepileupto(av, y);
! 1129: }
! 1130: if (m<0)
! 1131: y = FpXQ_mul(y,x,pol,p);
! 1132: if (low_stack(lim, stack_lim(av,1)))
! 1133: {
! 1134: if(DEBUGMEM>1) err(warnmem,"[2]: FpXQ_pow");
! 1135: y = gerepileupto(av, y);
! 1136: }
! 1137: }
! 1138: if (--i == 0) break;
! 1139: m = *++p1, j = BITS_IN_LONG;
! 1140: }
! 1141: }
! 1142: return gerepileupto(ltop,y);
! 1143: }
! 1144:
! 1145: /*******************************************************************/
! 1146: /* */
! 1147: /* Fp[X][Y] */
! 1148: /* */
! 1149: /*******************************************************************/
! 1150: /*Polynomials whose coefficients are either polynomials or integers*/
! 1151:
! 1152: GEN
! 1153: FpXX_red(GEN z, GEN p)
! 1154: {
! 1155: GEN res;
! 1156: int i;
! 1157: res=cgetg(lgef(z),t_POL);
! 1158: res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z));
! 1159: for(i=2;i<lgef(res);i++)
! 1160: if (typ(z[i])!=t_INT)
! 1161: {
! 1162: ulong av=avma;
! 1163: res[i]=(long)FpX_red((GEN)z[i],p);
! 1164: if (lgef(res[i])<=3)
! 1165: {
! 1166: if (lgef(res[i])==2) {avma=av;res[i]=zero;}
! 1167: else res[i]=lpilecopy(av,gmael(res,i,2));
! 1168: }
! 1169: }
! 1170: else
! 1171: res[i]=lmodii((GEN)z[i],p);
! 1172: res=normalizepol_i(res,lgef(res));
! 1173: return res;
! 1174: }
! 1175:
! 1176: /*******************************************************************/
! 1177: /* */
! 1178: /* (Fp[X]/(Q))[Y] */
! 1179: /* */
! 1180: /*******************************************************************/
! 1181:
! 1182: extern GEN to_Kronecker(GEN P, GEN Q);
! 1183: GEN /*Somewhat copy-pasted...*/
! 1184: /*Not malloc nor warn-clean.*/
! 1185: FpXQX_from_Kronecker(GEN z, GEN pol, GEN p)
! 1186: {
! 1187: long i,j,lx,l = lgef(z), N = (degpol(pol)<<1) + 1;
! 1188: GEN x, t = cgetg(N,t_POL);
! 1189: t[1] = pol[1] & VARNBITS;
! 1190: lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
! 1191: if (isonstack(pol)) pol = gcopy(pol);
! 1192: for (i=2; i<lx+2; i++)
! 1193: {
! 1194: for (j=2; j<N; j++) t[j] = z[j];
! 1195: z += (N-2);
! 1196: x[i] = (long)FpX_res(normalizepol_i(t,N), pol, p);
! 1197: }
! 1198: N = (l-2) % (N-2) + 2;
! 1199: for (j=2; j<N; j++) t[j] = z[j];
! 1200: x[i] = (long)FpX_res(normalizepol_i(t,N), pol, p);
! 1201: return normalizepol_i(x, i+1);
! 1202: }
! 1203: /*Unused/untested*/
! 1204: GEN
! 1205: FpXQX_red(GEN z, GEN T, GEN p)
! 1206: {
! 1207: GEN res;
! 1208: int i;
! 1209: res=cgetg(lgef(z),t_POL);
! 1210: res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z));
! 1211: for(i=2;i<lgef(res);i++)
! 1212: if (typ(z[i])!=t_INT)
! 1213: {
! 1214: if (T)
! 1215: res[i]=(long)FpX_res((GEN)z[i],T,p);
! 1216: else
! 1217: res[i]=(long)FpX_red((GEN)z[i],p);
! 1218: }
! 1219: else
! 1220: res[i]=lmodii((GEN)z[i],p);
! 1221: res=normalizepol_i(res,lgef(res));
! 1222: return res;
! 1223: }
! 1224: /* if T = NULL, call FpX_mul */
! 1225: GEN
! 1226: FpXQX_mul(GEN x, GEN y, GEN T, GEN p)
! 1227: {
! 1228: ulong ltop;
! 1229: GEN z,kx,ky;
! 1230: long vx;
! 1231: if (!T) return FpX_mul(x,y,p);
! 1232: ltop = avma;
! 1233: vx = min(varn(x),varn(y));
! 1234: kx= to_Kronecker(x,T);
! 1235: ky= to_Kronecker(y,T);
! 1236: z = quickmul(ky+2, kx+2, lgef(ky)-2, lgef(kx)-2);
! 1237: z = FpX_red(z,p);
! 1238: z = FpXQX_from_Kronecker(z,T,p);
! 1239: setvarn(z,vx);/*quickmul and Fq_from_Kronecker are not varn-clean*/
! 1240: return gerepileupto(ltop,z);
! 1241: }
! 1242: GEN/*Unused/untested*/
! 1243: FpXQX_sqr(GEN x, GEN T, GEN p)
! 1244: {
! 1245: ulong ltop=avma;
! 1246: GEN z,kx;
! 1247: long vx=varn(x);
! 1248: kx= to_Kronecker(x,T);
! 1249: z = quicksqr(kx+2, lgef(kx)-2);
! 1250: z = FpX_red(z,p);
! 1251: z = FpXQX_from_Kronecker(z,T,p);
! 1252: setvarn(z,vx);/*quickmul and Fq_from_Kronecker are nor varn-clean*/
! 1253: return gerepileupto(ltop,z);
! 1254: }
! 1255: GEN
! 1256: FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p)
! 1257: {
! 1258: int i, lP = lgef(P);
! 1259: GEN res = cgetg(lP,t_POL);
! 1260: res[1] = evalsigne(1) | evalvarn(varn(P)) | evallgef(lP);
! 1261: for(i=2; i<lP; i++)
! 1262: res[i] = (long)FpXQ_mul(U,(GEN)P[i], T,p);
! 1263: return normalizepol_i(res,lgef(res));
! 1264: }
! 1265:
! 1266: /* a X^degpol, assume degpol >= 0 */
! 1267: static GEN
! 1268: monomial(GEN a, int degpol, int v)
! 1269: {
! 1270: long i, lP = degpol+3;
! 1271: GEN P = cgetg(lP, t_POL);
! 1272: P[1] = evalsigne(1) | evalvarn(v) | evallgef(lP);
! 1273: lP--; P[lP] = (long)a;
! 1274: for (i=2; i<lP; i++) P[i] = zero;
! 1275: return P;
! 1276: }
! 1277:
! 1278: /* return NULL if Euclidean remainder sequence fails (==> T reducible mod p)
! 1279: * last non-zero remainder otherwise */
! 1280: GEN
! 1281: FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p)
! 1282: {
! 1283: ulong ltop = avma, btop, st_lim;
! 1284: long dg, vx = varn(P);
! 1285: GEN U, q;
! 1286: P = FpXX_red(P, p);
! 1287: Q = FpXX_red(Q, p);
! 1288: if (!signe(P)) return gerepileupto(ltop, Q);
! 1289: if (!signe(Q)) { avma = (ulong)P; return P; }
! 1290: T = FpX_red(T, p);
! 1291:
! 1292: btop = avma; st_lim = stack_lim(btop, 1);
! 1293: dg = lgef(P)-lgef(Q);
! 1294: if (dg < 0) { swap(P, Q); dg = -dg; }
! 1295: for(;;)
! 1296: {
! 1297: U = FpXQ_invsafe(leading_term(Q), T, p);
! 1298: if (!U) { avma = ltop; return NULL; }
! 1299: do /* set P := P % Q */
! 1300: {
! 1301: q = FpXQ_mul(U, gneg(leading_term(P)), T, p);
! 1302: P = gadd(P, FpXQX_mul(monomial(q, dg, vx), Q, T, p));
! 1303: P = FpXQX_red(P, T, p); /* wasteful, but negligible */
! 1304: dg = lgef(P)-lgef(Q);
! 1305: } while (dg >= 0);
! 1306: if (!signe(P)) break;
! 1307:
! 1308: if (low_stack(st_lim, stack_lim(btop, 1)))
! 1309: {
! 1310: GEN *bptr[2]; bptr[0]=&P; bptr[1]=&Q;
! 1311: if (DEBUGLEVEL>1) err(warnmem,"FpXQX_safegcd");
! 1312: gerepilemany(btop, bptr, 2);
! 1313: }
! 1314: swap(P, Q); dg = -dg;
! 1315: }
! 1316: Q = FpXQX_FpXQ_mul(Q, U, T, p); /* normalize GCD */
! 1317: return gerepileupto(ltop, Q);
! 1318: }
! 1319:
! 1320: /*******************************************************************/
! 1321: /* */
! 1322: /* Fq[X] */
! 1323: /* */
! 1324: /*******************************************************************/
! 1325:
! 1326: /*Well nothing, for now. So we reuse FpXQX*/
! 1327: #define FqX_mul FpXQX_mul
! 1328: #define FqX_sqr FpXQX_sqr
! 1329: #define FqX_red FpXQX_red
! 1330: static GEN
! 1331: Fq_neg(GEN x, GEN T, GEN p)/*T is not used but it is for consistency*/
! 1332: {
! 1333: switch(typ(x)==t_POL)
! 1334: {
! 1335: case 0: return signe(x)?subii(p,x):gzero;
! 1336: case 1: return FpX_neg(x,p);
! 1337: }
! 1338: return NULL;
! 1339: }
! 1340: static GEN modulo,Tmodulo;
! 1341: static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodulo,modulo);}
! 1342: GEN
! 1343: FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v)
! 1344: {
! 1345: ulong ltop=avma;
! 1346: long k;
! 1347: GEN W=cgetg(lg(V),t_VEC);
! 1348: for(k=1;k<lg(V);k++)
! 1349: W[k]=(long)deg1pol(gun,Fq_neg((GEN)V[k],Tp,p),v);
! 1350: modulo=p;Tmodulo=Tp;
! 1351: W=divide_conquer_prod(W,&fgmul);
! 1352: return gerepileupto(ltop,W);
! 1353: }
! 1354:
! 1355: /*******************************************************************/
! 1356: /* */
! 1357: /* n-th ROOT in Fq */
! 1358: /* */
! 1359: /*******************************************************************/
! 1360: /*NO clean malloc*/
! 1361: static GEN fflgen(GEN l, long e, GEN r, GEN T ,GEN p, GEN *zeta)
! 1362: {
! 1363: ulong av1;
! 1364: GEN z,m,m1;
! 1365: long x=varn(T),k,u,v,pp,i;
! 1366: if (is_bigint(p))
! 1367: pp=VERYBIGINT;
! 1368: else
! 1369: pp=itos(p);
! 1370: z=(lgef(T)==4)?polun[x]:polx[x];
! 1371:
! 1372: av1 = avma;
! 1373: for (k=1; ; k++)
! 1374: {
! 1375: u=k;v=0;
! 1376: while (u%pp==0){u/=pp;v++;}
! 1377: if(!v)
! 1378: z=gadd(z,gun);
! 1379: else
! 1380: {
! 1381: z=gadd(z,gpowgs(polx[x],v));
! 1382: if (DEBUGLEVEL>=6)
! 1383: fprintferr("FF l-Gen:next %Z",z);
! 1384: }
! 1385: m1 = m = FpXQ_pow(z,r,T,p);
! 1386: for (i=1; i<e; i++)
! 1387: if (gcmp1(m=FpXQ_pow(m,l,T,p))) break;
! 1388: if (i==e) break;
! 1389: avma = av1;
! 1390: }
! 1391: *zeta=m;
! 1392: return m1;
! 1393: }
! 1394: /* resoud x^l=a mod (p,T)
! 1395: * l doit etre premier
! 1396: * q=p^degpol(T)-1
! 1397: * q=(l^e)*r
! 1398: * e>=1
! 1399: * pgcd(r,l)=1
! 1400: * m=y^(q/l)
! 1401: * y n'est pas une puissance l-ieme
! 1402: * m!=1
! 1403: * ouf!
! 1404: */
! 1405: GEN
! 1406: ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, GEN r, GEN y, GEN m)
! 1407: {
! 1408: ulong av = avma,lim;
! 1409: long i,k;
! 1410: GEN p1,p2,u1,u2,v,w,z;
! 1411:
! 1412: bezout(r,l,&u1,&u2);
! 1413: v=FpXQ_pow(a,u2,T,p);
! 1414: w=FpXQ_pow(a,modii(mulii(negi(u1),r),q),T,p);
! 1415: lim = stack_lim(av,1);
! 1416: while (!gcmp1(w))
! 1417: {
! 1418: /* if p is not prime, next loop will not end */
! 1419: k=0;
! 1420: p1=w;
! 1421: do
! 1422: {
! 1423: z=p1;
! 1424: p1=FpXQ_pow(p1,l,T,p);
! 1425: k++;
! 1426: }while(!gcmp1(p1));
! 1427: if (k==e) { avma=av; return NULL; }
! 1428: p2 = FpXQ_mul(z,m,T,p);
! 1429: for(i=1; !gcmp1(p2); i++) p2 = FpXQ_mul(p2,m,T,p);/*should be a baby step
! 1430: giant step instead*/
! 1431: p1= FpXQ_pow(y,modii(mulsi(i,gpowgs(l,e-k-1)),q),T,p);
! 1432: m = FpXQ_pow(m,stoi(i),T,p);
! 1433: e = k;
! 1434: v = FpXQ_mul(p1,v,T,p);
! 1435: y = FpXQ_pow(p1,l,T,p);
! 1436: w = FpXQ_mul(y,w,T,p);
! 1437: if (low_stack(lim, stack_lim(av,1)))
! 1438: {
! 1439: GEN *gptr[4];
! 1440: if(DEBUGMEM>1) err(warnmem,"ffsqrtlmod");
! 1441: gptr[0]=&y; gptr[1]=&v; gptr[2]=&w; gptr[3]=&m;
! 1442: gerepilemany(av,gptr,4);
! 1443: }
! 1444: }
! 1445: return gerepilecopy(av,v);
! 1446: }
! 1447: /* n is an integer, a is in Fp[X]/(T), p is prime, T is irreducible mod p
! 1448:
! 1449: return a solution of
! 1450:
! 1451: x^n=a mod p
! 1452:
! 1453: 1)If there is no solution return NULL and if zetan is not NULL set zetan to gzero.
! 1454:
! 1455: 2) If there is solution there are exactly m=gcd(p-1,n) of them.
! 1456:
! 1457: If zetan is not NULL, zetan is set to a primitive mth root of unity so that
! 1458: the set of solutions is {x*zetan^k;k=0 to m-1}
! 1459:
! 1460: If a=0 ,return 0 and if zetan is not NULL zetan is set to gun
! 1461: */
! 1462: GEN ffsqrtnmod(GEN a, GEN n, GEN T, GEN p, GEN *zetan)
! 1463: {
! 1464: ulong ltop=avma,lbot=0,av1,lim;
! 1465: long i,j,e;
! 1466: GEN m,u1,u2;
! 1467: GEN q,r,zeta,y,l,z;
! 1468: GEN *gptr[2];
! 1469: if (typ(a) != t_POL || typ(n) != t_INT || typ(T) != t_POL || typ(p)!=t_INT)
! 1470: err(typeer,"ffsqrtnmod");
! 1471: if (lgef(T)==3)
! 1472: err(constpoler,"ffsqrtnmod");
! 1473: if(!signe(n))
! 1474: err(talker,"1/0 exponent in ffsqrtnmod");
! 1475: if(gcmp1(n)) {if (zetan) *zetan=gun;return gcopy(a);}
! 1476: if(gcmp0(a)) {if (zetan) *zetan=gun;return gzero;}
! 1477: q=addsi(-1,gpowgs(p,degpol(T)));
! 1478: m=bezout(n,q,&u1,&u2);
! 1479: if (gcmp(m,n))
! 1480: {
! 1481: GEN b=modii(u1,q);
! 1482: lbot=avma;
! 1483: a=FpXQ_pow(a,b,T,p);
! 1484: }
! 1485: if (zetan) z=polun[varn(T)];
! 1486: lim=stack_lim(ltop,1);
! 1487: if (!gcmp1(m))
! 1488: {
! 1489: m=decomp(m);
! 1490: av1=avma;
! 1491: for (i = lg(m[1])-1; i; i--)
! 1492: {
! 1493: l=gcoeff(m,i,1); j=itos(gcoeff(m,i,2));
! 1494: e=pvaluation(q,l,&r);
! 1495: y=fflgen(l,e,r,T,p,&zeta);
! 1496: if (zetan) z=FpXQ_mul(z,FpXQ_pow(y,gpowgs(l,e-j),T,p),T,p);
! 1497: do
! 1498: {
! 1499: lbot=avma;
! 1500: a=ffsqrtlmod(a,l,T,p,q,e,r,y,zeta);
! 1501: if (!a){avma=ltop;return NULL;}
! 1502: j--;
! 1503: }while (j);
! 1504: if (low_stack(lim, stack_lim(ltop,1)))
! 1505: /* n can have lots of prime factors*/
! 1506: {
! 1507: if(DEBUGMEM>1) err(warnmem,"ffsqrtnmod");
! 1508: if (zetan)
! 1509: {
! 1510: z=gcopy(z);
! 1511: gptr[0]=&a;gptr[1]=&z;
! 1512: gerepilemanysp(av1,lbot,gptr,2);
! 1513: }
! 1514: else
! 1515: a=gerepileupto(av1,a);
! 1516: lbot=av1;
! 1517: }
! 1518: }
! 1519: }
! 1520: if (zetan)
! 1521: {
! 1522: *zetan=gcopy(z);
! 1523: gptr[0]=&a;gptr[1]=zetan;
! 1524: gerepilemanysp(ltop,lbot,gptr,2);
! 1525: }
! 1526: else
! 1527: a=gerepileupto(ltop,a);
! 1528: return a;
! 1529: }
! 1530: /*******************************************************************/
! 1531: /* Isomorphisms between finite fields */
! 1532: /* */
! 1533: /*******************************************************************/
! 1534: GEN
! 1535: matrixpow(long n, long m, GEN y, GEN P,GEN l)
! 1536: {
! 1537: return vecpol_to_mat(FpXQ_powers(y,m-1,P,l),n);
! 1538: }
! 1539: /* compute the reciprocical isomorphism of S mod T,p, i.e. V such that
! 1540: V(S)=X mod T,p*/
! 1541: GEN
! 1542: Fp_inv_isom(GEN S,GEN T, GEN p)
! 1543: {
! 1544: ulong ltop = avma, lbot;
! 1545: GEN M, V;
! 1546: int n, i;
! 1547: long x;
! 1548: x = varn(T);
! 1549: n = degpol(T);
! 1550: M = matrixpow(n,n,S,T,p);
! 1551: V = cgetg(n + 1, t_COL);
! 1552: for (i = 1; i <= n; i++)
! 1553: V[i] = zero;
! 1554: V[2] = un;
! 1555: V = FpM_invimage(M,V,p);
! 1556: lbot = avma;
! 1557: V = gtopolyrev(V, x);
! 1558: return gerepile(ltop, lbot, V);
! 1559: }
! 1560: #if 0
! 1561: /*Old, trivial, implementation*/
! 1562: GEN
! 1563: intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU)
! 1564: {
! 1565: ulong ltop=avma;
! 1566: long vp=varn(P), np=degpol(P);
! 1567: GEN A;
! 1568: A=FqM_ker(gaddmat(gneg(polx[MAXVARN]), *MA),lU,l);
! 1569: if (lg(A)!=2)
! 1570: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
! 1571: ,l,polx[vp],P);
! 1572: A=gmul((GEN)A[1],gmodulcp(gmodulcp(gun,l),U));
! 1573: return gerepileupto(ltop,gtopolyrev(A,vp));
! 1574: }
! 1575: #endif
! 1576: /* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in
! 1577: * FFp[X]/T. Compute P(M)
! 1578: * not stack clean
! 1579: */
! 1580: static GEN
! 1581: polfrobenius(GEN M, GEN p, long r, long v)
! 1582: {
! 1583: GEN V,W;
! 1584: long i;
! 1585: V = cgetg(r+2,t_VEC);
! 1586: V[1] = (long) polx[v];
! 1587: if (r == 0) return V;
! 1588: V[2] = (long) gtopolyrev((GEN)M[2],v);
! 1589: W = (GEN) M[2];
! 1590: for (i = 3; i <= r+1; ++i)
! 1591: {
! 1592: W = FpM_FpV_mul(M,W,p);
! 1593: V[i] = (long) gtopolyrev(W,v);
! 1594: }
! 1595: return V;
! 1596: }
! 1597:
! 1598: static GEN
! 1599: matpolfrobenius(GEN V, GEN P, GEN T, GEN p)
! 1600: {
! 1601: long i;
! 1602: long l=degpol(T);
! 1603: long v=varn(T);
! 1604: GEN M,W;
! 1605: GEN PV=gtovec(P);
! 1606: PV=cgetg(degpol(P)+2,t_VEC);
! 1607: for(i=1;i<lg(PV);i++)
! 1608: PV[i]=P[1+i];
! 1609: M=cgetg(l+1,t_VEC);
! 1610: M[1]=(long)scalarpol(poleval(P,gun),v);
! 1611: M[2]=(long)FpXV_FpV_dotproduct(V,PV,p);
! 1612: W=cgetg(lg(V),t_VEC);
! 1613: for(i=1;i<lg(W);i++)
! 1614: W[i]=V[i];
! 1615: for(i=3;i<=l;i++)
! 1616: {
! 1617: long j;
! 1618: for(j=1;j<lg(W);j++)
! 1619: W[j]=(long)FpXQ_mul((GEN)W[j],(GEN)V[j],T,p);
! 1620: M[i]=(long)FpXV_FpV_dotproduct(W,PV,p);
! 1621: }
! 1622: return vecpol_to_mat(M,l);
! 1623: }
! 1624:
! 1625: /* Essentially we want to compute
! 1626: * FqM_ker(MA-polx[MAXVARN],lU,l)
! 1627: * To avoid use of matrix in Fq we procede as follows:
! 1628: * We compute FpM_ker(U(MA)) and then we recover
! 1629: * the eigen value by Galois action, see formula.
! 1630: */
! 1631: static GEN
! 1632: intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU)
! 1633: {
! 1634: ulong ltop=avma;
! 1635: long vp=varn(P);
! 1636: long vu=varn(lU), r=degpol(lU);
! 1637: long i;
! 1638: GEN A, R, M, ib0, V;
! 1639: if (DEBUGLEVEL>=4) timer2();
! 1640: V=polfrobenius(MA,l,r,varn(U));
! 1641: if (DEBUGLEVEL>=4) msgtimer("pol[frobenius]");
! 1642: M=matpolfrobenius(V,lU,P,l);
! 1643: if (DEBUGLEVEL>=4) msgtimer("matrix cyclo");
! 1644: A=FpM_ker(M,l);
! 1645: if (DEBUGLEVEL>=4) msgtimer("kernel");
! 1646: A=gerepileupto(ltop,A);
! 1647: if (lg(A)!=r+1)
! 1648: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
! 1649: ,l,polx[vp],P);
! 1650: /*The formula is
! 1651: * a_{r-1}=-\phi(a_0)/b_0
! 1652: * a_{i-1}=\phi(a_i)+b_ia_{r-1} i=r-1 to 1
! 1653: * Where a_0=A[1] and b_i=U[i+2]
! 1654: */
! 1655: ib0=negi(mpinvmod((GEN)lU[2],l));
! 1656: R=cgetg(r+1,t_MAT);
! 1657: R[1]=A[1];
! 1658: R[r]=(long)FpM_FpV_mul(MA,gmul((GEN)A[1],ib0),l);
! 1659: for(i=r-1;i>1;i--)
! 1660: R[i]=(long)FpV_red(gadd(FpM_FpV_mul(MA,(GEN)R[i+1],l),
! 1661: gmul((GEN)lU[i+2],(GEN)R[r])),l);
! 1662: R=gtrans_i(R);
! 1663: for(i=1;i<lg(R);i++)
! 1664: R[i]=(long)gtopolyrev((GEN)R[i],vu);
! 1665: A=gtopolyrev(R,vp);
! 1666: A=gmul(A,gmodulcp(gmodulcp(gun,l),U));
! 1667: return gerepileupto(ltop,A);
! 1668: }
! 1669:
! 1670: /* n must divide both the degree of P and Q. Compute SP and SQ such
! 1671: that the subfield of FF_l[X]/(P) generated by SP and the subfield of
! 1672: FF_l[X]/(Q) generated by SQ are isomorphic of degree n. P and Q do
! 1673: not need to be of the same variable. if MA (resp. MB) is not NULL,
! 1674: must be the matrix of the frobenius map in FF_l[X]/(P) (resp.
! 1675: FF_l[X]/(Q) ). */
! 1676: /* Note on the implementation choice:
! 1677: * We assume the prime p is very large
! 1678: * so we handle Frobenius as matrices.
! 1679: */
! 1680: void
! 1681: Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN *SQ, GEN MA, GEN MB)
! 1682: {
! 1683: ulong ltop=avma,lbot;
! 1684: long vp,vq,np,nq,e,pg;
! 1685: GEN q;
! 1686: GEN A,B,Ap,Bp;
! 1687: GEN *gptr[2];
! 1688: vp=varn(P);vq=varn(Q);
! 1689: np=degpol(P);nq=degpol(Q);
! 1690: if (np<=0 || nq<=0 || n<=0 || np%n!=0 || nq%n!=0)
! 1691: err(talker,"bad degrees in Fp_intersect: %d,%d,%d",n,degpol(P),degpol(Q));
! 1692: e=pvaluation(stoi(n),l,&q);
! 1693: pg=itos(q);
! 1694: avma=ltop;
! 1695: if (DEBUGLEVEL>=4) timer2();
! 1696: if(!MA) MA=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l);
! 1697: if(!MB) MB=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l);
! 1698: if (DEBUGLEVEL>=4) msgtimer("matrixpow");
! 1699: A=Ap=zeropol(vp);
! 1700: B=Bp=zeropol(vq);
! 1701: if (pg>1)
! 1702: {
! 1703: if (gcmp0(modis(addis(l,-1),pg)))
! 1704: /*We do not need to use relative extension in this setting, so
! 1705: we don't. (Well,now that we don't in the other case also, it is more
! 1706: dubious to treat cases apart...)*/
! 1707: {
! 1708: GEN L,An,Bn,ipg,z;
! 1709: z=rootmod(cyclo(pg,-1),l);
! 1710: if (lg(z)<2) err(talker,"%Z is not a prime in Fp_intersect",l);
! 1711: z=negi(lift((GEN)z[1]));
! 1712: ipg=stoi(pg);
! 1713: if (DEBUGLEVEL>=4) timer2();
! 1714: A=FpM_ker(gaddmat(z, MA),l);
! 1715: if (lg(A)!=2)
! 1716: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
! 1717: ,l,polx[vp],P);
! 1718: A=gtopolyrev((GEN)A[1],vp);
! 1719: B=FpM_ker(gaddmat(z, MB),l);
! 1720: if (lg(B)!=2)
! 1721: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
! 1722: ,l,polx[vq],Q);
! 1723: B=gtopolyrev((GEN)B[1],vq);
! 1724: if (DEBUGLEVEL>=4) msgtimer("FpM_ker");
! 1725: An=(GEN) FpXQ_pow(A,ipg,P,l)[2];
! 1726: Bn=(GEN) FpXQ_pow(B,ipg,Q,l)[2];
! 1727: z=modii(mulii(An,mpinvmod(Bn,l)),l);
! 1728: L=mpsqrtnmod(z,ipg,l,NULL);
! 1729: if ( !L )
! 1730: err(talker,"Polynomials not irreducible in Fp_intersect");
! 1731: if (DEBUGLEVEL>=4) msgtimer("mpsqrtnmod");
! 1732: B=FpX_Fp_mul(B,L,l);
! 1733: }
! 1734: else
! 1735: {
! 1736: GEN L,An,Bn,ipg,U,lU,z;
! 1737: U=gmael(factmod(cyclo(pg,MAXVARN),l),1,1);
! 1738: lU=lift(U);
! 1739: ipg=stoi(pg);
! 1740: A=intersect_ker(P, MA, l, U, lU);
! 1741: B=intersect_ker(Q, MB, l, U, lU);
! 1742: /*Somewhat ugly, but it is a proof that POLYMOD are useful and
! 1743: powerful.*/
! 1744: if (DEBUGLEVEL>=4) timer2();
! 1745: An=lift(lift((GEN)lift(gpowgs(gmodulcp(A,P),pg))[2]));
! 1746: Bn=lift(lift((GEN)lift(gpowgs(gmodulcp(B,Q),pg))[2]));
! 1747: if (DEBUGLEVEL>=4) msgtimer("pows [P,Q]");
! 1748: z=FpXQ_inv(Bn,lU,l);
! 1749: z=FpXQ_mul(An,z,lU,l);
! 1750: L=ffsqrtnmod(z,ipg,lU,l,NULL);
! 1751: if (DEBUGLEVEL>=4) msgtimer("ffsqrtn");
! 1752: if ( !L )
! 1753: err(talker,"Polynomials not irreducible in Fp_intersect");
! 1754: B=gsubst(lift(lift(gmul(B,L))),MAXVARN,gzero);
! 1755: A=gsubst(lift(lift(A)),MAXVARN,gzero);
! 1756: }
! 1757: }
! 1758: if (e!=0)
! 1759: {
! 1760: GEN VP,VQ,moinsun,Ay,By,lmun;
! 1761: int i,j;
! 1762: moinsun=stoi(-1);
! 1763: lmun=addis(l,-1);
! 1764: MA=gaddmat(moinsun,MA);
! 1765: MB=gaddmat(moinsun,MB);
! 1766: Ay=polun[vp];
! 1767: By=polun[vq];
! 1768: VP=cgetg(np+1,t_COL);
! 1769: VP[1]=un;
! 1770: for(i=2;i<=np;i++) VP[i]=zero;
! 1771: if (np==nq) VQ=VP;/*save memory*/
! 1772: else
! 1773: {
! 1774: VQ=cgetg(nq+1,t_COL);
! 1775: VQ[1]=un;
! 1776: for(i=2;i<=nq;i++) VQ[i]=zero;
! 1777: }
! 1778: for(j=0;j<e;j++)
! 1779: {
! 1780: if (j)
! 1781: {
! 1782: Ay=FpXQ_mul(Ay,FpXQ_pow(Ap,lmun,P,l),P,l);
! 1783: for(i=1;i<lgef(Ay)-1;i++) VP[i]=Ay[i+1];
! 1784: for(;i<=np;i++) VP[i]=zero;
! 1785: }
! 1786: Ap=FpM_invimage(MA,VP,l);
! 1787: Ap=gtopolyrev(Ap,vp);
! 1788: if (j)
! 1789: {
! 1790: By=FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l);
! 1791: for(i=1;i<lgef(By)-1;i++) VQ[i]=By[i+1];
! 1792: for(;i<=nq;i++) VQ[i]=zero;
! 1793: }
! 1794: Bp=FpM_invimage(MB,VQ,l);
! 1795: Bp=gtopolyrev(Bp,vq);
! 1796: if (DEBUGLEVEL>=4) msgtimer("FpM_invimage");
! 1797: }
! 1798: }/*FpX_add is not clean, so we must do it *before* lbot=avma*/
! 1799: A=FpX_add(A,Ap,NULL);
! 1800: B=FpX_add(B,Bp,NULL);
! 1801: lbot=avma;
! 1802: *SP=FpX_red(A,l);
! 1803: *SQ=FpX_red(B,l);
! 1804: gptr[0]=SP;gptr[1]=SQ;
! 1805: gerepilemanysp(ltop,lbot,gptr,2);
! 1806: }
! 1807: /* Let l be a prime number, P, Q in ZZ[X]. P and Q are both
! 1808: * irreducible modulo l and degree(P) divides degree(Q). Output a
! 1809: * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
! 1810: * that Q | P(R) mod l. If P and Q have the same degree, it is of course an
! 1811: * isomorphism. */
! 1812: GEN Fp_isom(GEN P,GEN Q,GEN l)
! 1813: {
! 1814: ulong ltop=avma;
! 1815: GEN SP,SQ,R;
! 1816: Fp_intersect(degpol(P),P,Q,l,&SP,&SQ,NULL,NULL);
! 1817: R=Fp_inv_isom(SP,P,l);
! 1818: R=FpX_FpXQ_compo(R,SQ,Q,l);
! 1819: return gerepileupto(ltop,R);
! 1820: }
! 1821: GEN
! 1822: Fp_factorgalois(GEN P,GEN l, long d, long w)
! 1823: {
! 1824: ulong ltop=avma;
! 1825: GEN R,V,ld,Tl;
! 1826: long n,k,m;
! 1827: long v;
! 1828: v=varn(P);
! 1829: n=degpol(P);
! 1830: m=n/d;
! 1831: ld=gpowgs(l,d);
! 1832: Tl=gcopy(P); setvarn(Tl,w);
! 1833: V=cgetg(m+1,t_VEC);
! 1834: V[1]=lpolx[w];
! 1835: for(k=2;k<=m;k++)
! 1836: V[k]=(long)FpXQ_pow((GEN)V[k-1],ld,Tl,l);
! 1837: R=FqV_roots_to_pol(V,Tl,l,v);
! 1838: return gerepileupto(ltop,R);
! 1839: }
! 1840: extern GEN mat_to_polpol(GEN x, long v,long w);
! 1841: extern GEN polpol_to_mat(GEN v, long n);
! 1842: /* P,Q irreducible over F_l. Factor P over FF_l[X] / Q [factors are ordered as
! 1843: * a Frobenius cycle] */
! 1844: GEN
! 1845: Fp_factor_irred(GEN P,GEN l, GEN Q)
! 1846: {
! 1847: ulong ltop=avma,av;
! 1848: GEN SP,SQ,MP,MQ,M,MF,E,V,IR,res;
! 1849: long np=degpol(P),nq=degpol(Q);
! 1850: long i,d=cgcd(np,nq);
! 1851: long vp=varn(P),vq=varn(Q);
! 1852: if (d==1)
! 1853: {
! 1854: res=cgetg(2,t_COL);
! 1855: res[1]=lcopy(P);
! 1856: return res;
! 1857: }
! 1858: MF=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l);
! 1859: Fp_intersect(d,P,Q,l,&SP,&SQ,NULL,MF);
! 1860: av=avma;
! 1861: E=Fp_factorgalois(P,l,d,vq);
! 1862: E=polpol_to_mat(E,np);
! 1863: MP = matrixpow(np,d,SP,P,l);
! 1864: IR = (GEN)FpM_sindexrank(MP,l)[1];
! 1865: E = rowextract_p(E, IR);
! 1866: M = rowextract_p(MP,IR);
! 1867: M = FpM_inv(M,l);
! 1868: MQ = matrixpow(nq,d,SQ,Q,l);
! 1869: M = FpM_mul(MQ,M,l);
! 1870: M = FpM_mul(M,E,l);
! 1871: M = gerepileupto(av,M);
! 1872: V = cgetg(d+1,t_VEC);
! 1873: V[1]=(long)M;
! 1874: for(i=2;i<=d;i++)
! 1875: V[i]=(long)FpM_mul(MF,(GEN)V[i-1],l);
! 1876: res=cgetg(d+1,t_COL);
! 1877: for(i=1;i<=d;i++)
! 1878: res[i]=(long)mat_to_polpol((GEN)V[i],vp,vq);
! 1879: return gerepilecopy(ltop,res);
! 1880: }
! 1881: GEN Fp_factor_rel0(GEN P,GEN l, GEN Q)
! 1882: {
! 1883: ulong ltop=avma,tetpil;
! 1884: GEN V,ex,F,y,R;
! 1885: long n,nbfact=0,nmax=lgef(P)-2;
! 1886: long i;
! 1887: F=factmod0(P,l);
! 1888: n=lg((GEN)F[1]);
! 1889: V=cgetg(nmax,t_VEC);
! 1890: ex=cgetg(nmax,t_VECSMALL);
! 1891: for(i=1;i<n;i++)
! 1892: {
! 1893: int j,r;
! 1894: R=Fp_factor_irred(gmael(F,1,i),l,Q);
! 1895: r=lg(R);
! 1896: for (j=1;j<r;j++)
! 1897: {
! 1898: V[++nbfact]=R[j];
! 1899: ex[nbfact]=mael(F,2,i);
! 1900: }
! 1901: }
! 1902: setlg(V,nbfact+1);
! 1903: setlg(ex,nbfact+1);
! 1904: tetpil=avma; y=cgetg(3,t_VEC);
! 1905: y[1]=lcopy(V);
! 1906: y[2]=lcopy(ex);
! 1907: (void)sort_factor(y,cmp_pol);
! 1908: return gerepile(ltop,tetpil,y);
! 1909: }
! 1910: GEN Fp_factor_rel(GEN P, GEN l, GEN Q)
! 1911: {
! 1912: long tetpil,av=avma;
! 1913: long nbfact;
! 1914: long j;
! 1915: GEN y,u,v;
! 1916: GEN z=Fp_factor_rel0(P,l,Q),t = (GEN) z[1],ex = (GEN) z[2];
! 1917: nbfact=lg(t);
! 1918: tetpil=avma; y=cgetg(3,t_MAT);
! 1919: u=cgetg(nbfact,t_COL); y[1]=(long)u;
! 1920: v=cgetg(nbfact,t_COL); y[2]=(long)v;
! 1921: for (j=1; j<nbfact; j++)
! 1922: {
! 1923: u[j] = lcopy((GEN)t[j]);
! 1924: v[j] = lstoi(ex[j]);
! 1925: }
! 1926: return gerepile(av,tetpil,y);
! 1927: }
! 1928:
! 1929: /*******************************************************************/
! 1930: extern int ff_poltype(GEN *x, GEN *p, GEN *pol);
! 1931:
! 1932: static GEN
! 1933: to_intmod(GEN x, GEN p)
! 1934: {
! 1935: GEN a = cgetg(3,t_INTMOD);
! 1936: a[1] = (long)p;
! 1937: a[2] = lmodii(x,p); return a;
! 1938: }
! 1939:
! 1940: /* z in Z[X], return z * Mod(1,p), normalized*/
! 1941: GEN
! 1942: FpX(GEN z, GEN p)
! 1943: {
! 1944: long i,l = lgef(z);
! 1945: GEN x = cgetg(l,t_POL);
! 1946: if (isonstack(p)) p = icopy(p);
! 1947: for (i=2; i<l; i++) x[i] = (long)to_intmod((GEN)z[i], p);
! 1948: x[1] = z[1]; return normalizepol_i(x,l);
! 1949: }
! 1950:
! 1951: /* z in Z^n, return z * Mod(1,p), normalized*/
! 1952: GEN
! 1953: FpV(GEN z, GEN p)
! 1954: {
! 1955: long i,l = lg(z);
! 1956: GEN x = cgetg(l,typ(z));
! 1957: if (isonstack(p)) p = icopy(p);
! 1958: for (i=1; i<l; i++) x[i] = (long)to_intmod((GEN)z[i], p);
! 1959: return x;
! 1960: }
! 1961: /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
! 1962: GEN
! 1963: FpM(GEN z, GEN p)
! 1964: {
! 1965: long i,j,l = lg(z), m = lg((GEN)z[1]);
! 1966: GEN x,y,zi;
! 1967: if (isonstack(p)) p = icopy(p);
! 1968: x = cgetg(l,t_MAT);
! 1969: for (i=1; i<l; i++)
! 1970: {
! 1971: x[i] = lgetg(m,t_COL);
! 1972: y = (GEN)x[i];
! 1973: zi= (GEN)z[i];
! 1974: for (j=1; j<m; j++) y[j] = (long)to_intmod((GEN)zi[j], p);
! 1975: }
! 1976: return x;
! 1977: }
! 1978: /* z in Z[X] or Z, return lift(z * Mod(1,p)), normalized*/
! 1979: GEN
! 1980: FpX_red(GEN z, GEN p)
! 1981: {
! 1982: long i,l;
! 1983: GEN x;
! 1984: if (typ(z) == t_INT) return modii(z,p);
! 1985: l = lgef(z);
! 1986: x = cgetg(l,t_POL);
! 1987: for (i=2; i<l; i++) x[i] = lmodii((GEN)z[i],p);
! 1988: x[1] = z[1]; return normalizepol_i(x,l);
! 1989: }
! 1990:
! 1991: GEN
! 1992: FpXV_red(GEN z, GEN p)
! 1993: {
! 1994: long i,l = lg(z);
! 1995: GEN x = cgetg(l, t_VEC);
! 1996: for (i=1; i<l; i++) x[i] = (long)FpX_red((GEN)z[i], p);
! 1997: return x;
! 1998: }
! 1999:
! 2000: /* z in Z^n, return lift(z * Mod(1,p)) */
! 2001: GEN
! 2002: FpV_red(GEN z, GEN p)
! 2003: {
! 2004: long i,l = lg(z);
! 2005: GEN x = cgetg(l,typ(z));
! 2006: for (i=1; i<l; i++) x[i] = lmodii((GEN)z[i],p);
! 2007: return x;
! 2008: }
! 2009:
! 2010: /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
! 2011: GEN
! 2012: FpM_red(GEN z, GEN p)
! 2013: {
! 2014: long i,j,l = lg(z), m = lg((GEN)z[1]);
! 2015: GEN x,y;
! 2016: x = cgetg(l,t_MAT);
! 2017: for (i=1; i<l; i++)
! 2018: {
! 2019: x[i]=lgetg(m,t_MAT);y=(GEN)x[i];
! 2020: for(j=1; j<m ; j++)
! 2021: y[j] = lmodii(gmael(z,i,j),p);
! 2022: }
! 2023: return x;
! 2024: }
! 2025:
! 2026: /* no garbage collection, divide by leading coeff, mod p */
! 2027: GEN
! 2028: FpX_normalize(GEN z, GEN p)
! 2029: {
! 2030: GEN p1 = leading_term(z);
! 2031: if (gcmp1(p1)) return z;
! 2032: return FpX_Fp_mul(z, mpinvmod(p1,p), p);
! 2033: }
! 2034:
! 2035: GEN
! 2036: FpXQX_normalize(GEN z, GEN T, GEN p)
! 2037: {
! 2038: GEN p1 = leading_term(z);
! 2039: if (gcmp1(p1)) return z;
! 2040: if (!T) return FpX_normalize(z,p);
! 2041: return FpXQX_FpXQ_mul(z, FpXQ_inv(p1,T,p), T,p);
! 2042: }
! 2043:
! 2044: GEN
! 2045: small_to_mat(GEN z)
! 2046: {
! 2047: long i, l = lg(z);
! 2048: GEN x = cgetg(l,t_MAT);
! 2049: for (i=1; i<l; i++) { x[i] = (long)gtovec((GEN)z[i]); settyp(x[i], t_COL); }
! 2050: return x;
! 2051: }
! 2052:
! 2053: GEN
! 2054: small_to_pol_i(GEN z, long l)
! 2055: {
! 2056: long i;
! 2057: GEN x = cgetg(l,t_POL);
! 2058: for (i=2; i<l; i++) x[i] = lstoi(z[i]);
! 2059: x[1] = z[1]; return x;
! 2060: }
! 2061:
! 2062: /* assume z[i] % p has been done */
! 2063: GEN
! 2064: small_to_pol(GEN z, long v)
! 2065: {
! 2066: GEN x = small_to_pol_i(z, lgef(z));
! 2067: setvarn(x,v); return x;
! 2068: }
! 2069:
! 2070: GEN
! 2071: pol_to_small(GEN x)
! 2072: {
! 2073: long i, lx = lgef(x);
! 2074: GEN a = u_allocpol(lx-3, 0);
! 2075: for (i=2; i<lx; i++) a[i] = itos((GEN)x[i]);
! 2076: return a;
! 2077: }
! 2078:
! 2079: /* z in Z[X,Y] representing an elt in F_p[X,Y] mod pol(Y) i.e F_q[X])
! 2080: * in Kronecker form. Recover the "real" z, normalized
! 2081: * If p = NULL, use generic functions and the coeff. ring implied by the
! 2082: * coefficients of z */
! 2083: GEN
! 2084: FqX_from_Kronecker(GEN z, GEN pol, GEN p)
! 2085: {
! 2086: long i,j,lx,l = lgef(z), N = (degpol(pol)<<1) + 1;
! 2087: GEN a,x, t = cgetg(N,t_POL);
! 2088: t[1] = pol[1] & VARNBITS;
! 2089: lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
! 2090: if (isonstack(pol)) pol = gcopy(pol);
! 2091: for (i=2; i<lx+2; i++)
! 2092: {
! 2093: a = cgetg(3,t_POLMOD); x[i] = (long)a;
! 2094: a[1] = (long)pol;
! 2095: for (j=2; j<N; j++) t[j] = z[j];
! 2096: z += (N-2);
! 2097: a[2] = (long)FpX_res(normalizepol_i(t,N), pol,p);
! 2098: }
! 2099: a = cgetg(3,t_POLMOD); x[i] = (long)a;
! 2100: a[1] = (long)pol;
! 2101: N = (l-2) % (N-2) + 2;
! 2102: for (j=2; j<N; j++) t[j] = z[j];
! 2103: a[2] = (long)FpX_res(normalizepol_i(t,N), pol,p);
! 2104: return normalizepol_i(x, i+1);
! 2105: }
! 2106:
! 2107: /* z in ?[X,Y] mod Q(Y) in Kronecker form ((subst(lift(z), x, y^(2deg(z)-1)))
! 2108: * Recover the "real" z, normalized */
! 2109: GEN
! 2110: from_Kronecker(GEN z, GEN pol)
! 2111: {
! 2112: return FqX_from_Kronecker(z,pol,NULL);
! 2113: }
! 2114:
! 2115: /*******************************************************************/
! 2116: /* */
! 2117: /* MODULAR GCD */
! 2118: /* */
! 2119: /*******************************************************************/
! 2120: ulong xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s);
! 2121: /* 1 / Mod(x,p) , or 0 if inverse doesn't exist */
! 2122: ulong
! 2123: u_invmod(ulong x, ulong p)
! 2124: {
! 2125: long s;
! 2126: ulong xv, xv1, g = xgcduu(p, x, 1, &xv, &xv1, &s);
! 2127: if (g != 1UL) return 0UL;
! 2128: xv = xv1 % p; if (s < 0) xv = p - xv;
! 2129: return xv;
! 2130: }
! 2131:
! 2132: #if 0
! 2133: static ulong
! 2134: umodratu(GEN a, ulong p)
! 2135: {
! 2136: if (typ(a) == t_INT)
! 2137: return umodiu(a,p);
! 2138: else { /* assume a a t_FRAC */
! 2139: ulong num = umodiu((GEN)a[1],p);
! 2140: ulong den = umodiu((GEN)a[2],p);
! 2141: return (ulong)mulssmod(num, u_invmod(den,p), p);
! 2142: }
! 2143: }
! 2144: #endif
! 2145:
! 2146: /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL/VECSMALL/INT */
! 2147: GEN
! 2148: u_Fp_FpX(GEN x, int malloc, ulong p)
! 2149: {
! 2150: long i, lx;
! 2151: GEN a;
! 2152:
! 2153: switch (typ(x))
! 2154: {
! 2155: case t_VECSMALL: return x;
! 2156: case t_INT: a = u_allocpol(0, malloc);
! 2157: a[2] = umodiu(x, p); return a;
! 2158: }
! 2159: lx = lgef(x); a = u_allocpol(lx-3, malloc);
! 2160: for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p);
! 2161: return u_normalizepol(a,lx);
! 2162: }
! 2163:
! 2164: GEN
! 2165: u_Fp_FpM(GEN x, ulong p)
! 2166: {
! 2167: long i,j,m,n = lg(x);
! 2168: GEN y = cgetg(n,t_MAT);
! 2169: if (n == 1) return y;
! 2170: m = lg(x[1]);
! 2171: for (j=1; j<n; j++)
! 2172: {
! 2173: y[j] = (long)cgetg(m,t_VECSMALL);
! 2174: for (i=1; i<m; i++) coeff(y,i,j) = (long)umodiu(gcoeff(x,i,j), p);
! 2175: }
! 2176: return y;
! 2177: }
! 2178:
! 2179: static GEN
! 2180: u_FpX_Fp_mul(GEN y, ulong x,ulong p, int malloc)
! 2181: {
! 2182: GEN z;
! 2183: int i, l;
! 2184: if (!x) return u_zeropol(malloc);
! 2185: l = lgef(y); z = u_allocpol(l-3, malloc);
! 2186: if (HIGHWORD(x | p))
! 2187: for(i=2; i<l; i++) z[i] = mulssmod(y[i], x, p);
! 2188: else
! 2189: for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
! 2190: return z;
! 2191: }
! 2192:
! 2193: GEN
! 2194: u_FpX_normalize(GEN z, ulong p)
! 2195: {
! 2196: long l = lgef(z)-1;
! 2197: ulong p1 = z[l]; /* leading term */
! 2198: if (p1 == 1) return z;
! 2199: return u_FpX_Fp_mul(z, u_invmod(p1,p), p, 0);
! 2200: }
! 2201:
! 2202: static GEN
! 2203: u_copy(GEN x, int malloc)
! 2204: {
! 2205: long i, l = lgef(x);
! 2206: GEN z = u_allocpol(l-3, malloc);
! 2207: for (i=2; i<l; i++) z[i] = x[i];
! 2208: return z;
! 2209: }
! 2210:
! 2211: /* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES */
! 2212: GEN
! 2213: u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *pr)
! 2214: {
! 2215: GEN z,q,c;
! 2216: long dx,dy,dz,i,j;
! 2217: ulong p1,inv;
! 2218:
! 2219: dy = degpol(y);
! 2220: if (!dy)
! 2221: {
! 2222: if (pr)
! 2223: {
! 2224: if (pr == ONLY_REM) return u_zeropol(malloc);
! 2225: *pr = u_zeropol(malloc);
! 2226: }
! 2227: if (y[2] == 1UL) return u_copy(x,malloc);
! 2228: return u_FpX_Fp_mul(x, u_invmod(y[2], p), p, malloc);
! 2229: }
! 2230: dx = degpol(x);
! 2231: dz = dx-dy;
! 2232: if (dz < 0)
! 2233: {
! 2234: if (pr)
! 2235: {
! 2236: c = u_copy(x, malloc);
! 2237: if (pr == ONLY_REM) return c;
! 2238: *pr = c;
! 2239: }
! 2240: return u_zeropol(malloc);
! 2241: }
! 2242: x += 2;
! 2243: y += 2;
! 2244: z = u_allocpol(dz, malloc || (pr == ONLY_REM)) + 2;
! 2245: inv = y[dy];
! 2246: if (inv != 1UL) inv = u_invmod(inv,p);
! 2247:
! 2248: if (u_OK_ULONG(p))
! 2249: {
! 2250: z[dz] = (inv*x[dx]) % p;
! 2251: for (i=dx-1; i>=dy; --i)
! 2252: {
! 2253: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
! 2254: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2255: {
! 2256: p1 += z[j]*y[i-j];
! 2257: if (p1 & HIGHBIT) p1 = p1 % p;
! 2258: }
! 2259: p1 %= p;
! 2260: z[i-dy] = p1? ((p - p1)*inv) % p: 0;
! 2261: }
! 2262: }
! 2263: else
! 2264: {
! 2265: z[dz] = mulssmod(inv, x[dx], p);
! 2266: for (i=dx-1; i>=dy; --i)
! 2267: {
! 2268: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
! 2269: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2270: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
! 2271: z[i-dy] = p1? mulssmod(p - p1, inv, p): 0;
! 2272: }
! 2273: }
! 2274: q = u_normalizepol(z-2, dz+3);
! 2275: if (!pr) return q;
! 2276:
! 2277: c = u_allocpol(dy,malloc) + 2;
! 2278: if (u_OK_ULONG(p))
! 2279: {
! 2280: for (i=0; i<dy; i++)
! 2281: {
! 2282: p1 = z[0]*y[i];
! 2283: for (j=1; j<=i && j<=dz; j++)
! 2284: {
! 2285: p1 += z[j]*y[i-j];
! 2286: if (p1 & HIGHBIT) p1 = p1 % p;
! 2287: }
! 2288: c[i] = subssmod(x[i], p1%p, p);
! 2289: }
! 2290: }
! 2291: else
! 2292: {
! 2293: for (i=0; i<dy; i++)
! 2294: {
! 2295: p1 = mulssmod(z[0],y[i],p);
! 2296: for (j=1; j<=i && j<=dz; j++)
! 2297: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
! 2298: c[i] = subssmod(x[i], p1, p);
! 2299: }
! 2300: }
! 2301: i=dy-1; while (i>=0 && !c[i]) i--;
! 2302: c = u_normalizepol(c-2, i+3);
! 2303: if (pr == ONLY_REM) { free(q); return c; }
! 2304: *pr = c; return q;
! 2305: }
! 2306:
! 2307: /* x and y in Z[X]. Possibly x in Z */
! 2308: GEN
! 2309: FpX_divres(GEN x, GEN y, GEN p, GEN *pr)
! 2310: {
! 2311: long vx,dx,dy,dz,i,j,av0,av,tetpil,sx,lrem;
! 2312: GEN z,p1,rem,lead;
! 2313:
! 2314: if (!p) return poldivres(x,y,pr);
! 2315: if (!signe(y)) err(talker,"division by zero in FpX_divres");
! 2316: vx=varn(x); dy=degpol(y); dx=(typ(x)==t_INT)? 0: degpol(x);
! 2317: if (dx < dy)
! 2318: {
! 2319: if (pr)
! 2320: {
! 2321: av0 = avma; x = FpX_red(x, p);
! 2322: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
! 2323: if (pr == ONLY_REM) return x;
! 2324: *pr = x;
! 2325: }
! 2326: return zeropol(vx);
! 2327: }
! 2328: lead = leading_term(y);
! 2329: if (!dy) /* y is constant */
! 2330: {
! 2331: if (pr && pr != ONLY_DIVIDES)
! 2332: {
! 2333: if (pr == ONLY_REM) return zeropol(vx);
! 2334: *pr = zeropol(vx);
! 2335: }
! 2336: if (gcmp1(lead)) return gcopy(x);
! 2337: av0 = avma; x = gmul(x, mpinvmod(lead,p)); tetpil = avma;
! 2338: return gerepile(av0,tetpil,FpX_red(x,p));
! 2339: }
! 2340: av0 = avma; dz = dx-dy;
! 2341: if (OK_ULONG(p))
! 2342: { /* assume ab != 0 mod p */
! 2343: ulong pp = (ulong)p[2];
! 2344: GEN a = u_Fp_FpX(x,1, pp);
! 2345: GEN b = u_Fp_FpX(y,1, pp);
! 2346: GEN zz= u_FpX_divrem(a,b,pp,1, pr);
! 2347: if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
! 2348: {
! 2349: rem = small_to_pol(*pr,vx);
! 2350: free(*pr); *pr = rem;
! 2351: }
! 2352: z = small_to_pol(zz,vx);
! 2353: free(zz); free(a); free(b); return z;
! 2354: }
! 2355: lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p));
! 2356: avma = av0;
! 2357: z=cgetg(dz+3,t_POL);
! 2358: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 2359: x += 2; y += 2; z += 2;
! 2360:
! 2361: p1 = (GEN)x[dx]; av = avma;
! 2362: z[dz] = lead? lpileupto(av, modii(mulii(p1,lead), p)): licopy(p1);
! 2363: for (i=dx-1; i>=dy; i--)
! 2364: {
! 2365: av=avma; p1=(GEN)x[i];
! 2366: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2367: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 2368: if (lead) p1 = mulii(p1,lead);
! 2369: tetpil=avma; z[i-dy] = lpile(av,tetpil,modii(p1, p));
! 2370: }
! 2371: if (!pr) { if (lead) gunclone(lead); return z-2; }
! 2372:
! 2373: rem = (GEN)avma; av = (long)new_chunk(dx+3);
! 2374: for (sx=0; ; i--)
! 2375: {
! 2376: p1 = (GEN)x[i];
! 2377: for (j=0; j<=i && j<=dz; j++)
! 2378: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 2379: tetpil=avma; p1 = modii(p1,p); if (signe(p1)) { sx = 1; break; }
! 2380: if (!i) break;
! 2381: avma=av;
! 2382: }
! 2383: if (pr == ONLY_DIVIDES)
! 2384: {
! 2385: if (lead) gunclone(lead);
! 2386: if (sx) { avma=av0; return NULL; }
! 2387: avma = (long)rem; return z-2;
! 2388: }
! 2389: lrem=i+3; rem -= lrem;
! 2390: rem[0]=evaltyp(t_POL) | evallg(lrem);
! 2391: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
! 2392: p1 = gerepile((long)rem,tetpil,p1);
! 2393: rem += 2; rem[i]=(long)p1;
! 2394: for (i--; i>=0; i--)
! 2395: {
! 2396: av=avma; p1 = (GEN)x[i];
! 2397: for (j=0; j<=i && j<=dz; j++)
! 2398: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
! 2399: tetpil=avma; rem[i]=lpile(av,tetpil, modii(p1,p));
! 2400: }
! 2401: rem -= 2;
! 2402: if (lead) gunclone(lead);
! 2403: if (!sx) normalizepol_i(rem, lrem);
! 2404: if (pr == ONLY_REM) return gerepileupto(av0,rem);
! 2405: *pr = rem; return z-2;
! 2406: }
! 2407:
! 2408: /* x and y in Z[Y][X]. Assume T irreducible mod p */
! 2409: GEN
! 2410: FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr)
! 2411: {
! 2412: long vx,dx,dy,dz,i,j,av0,av,tetpil,sx,lrem;
! 2413: GEN z,p1,rem,lead;
! 2414:
! 2415: if (!p) return poldivres(x,y,pr);
! 2416: if (!T) return FpX_divres(x,y,p,pr);
! 2417: if (!signe(y)) err(talker,"division by zero in FpX_divres");
! 2418: vx=varn(x); dy=degpol(y); dx=degpol(x);
! 2419: if (dx < dy)
! 2420: {
! 2421: if (pr)
! 2422: {
! 2423: av0 = avma; x = FpXQX_red(x, T, p);
! 2424: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
! 2425: if (pr == ONLY_REM) return x;
! 2426: *pr = x;
! 2427: }
! 2428: return zeropol(vx);
! 2429: }
! 2430: lead = leading_term(y);
! 2431: if (!dy) /* y is constant */
! 2432: {
! 2433: if (pr && pr != ONLY_DIVIDES)
! 2434: {
! 2435: if (pr == ONLY_REM) return zeropol(vx);
! 2436: *pr = zeropol(vx);
! 2437: }
! 2438: if (gcmp1(lead)) return gcopy(x);
! 2439: av0 = avma; x = gmul(x, FpXQ_inv(lead,T,p)); tetpil = avma;
! 2440: return gerepile(av0,tetpil,FpXQX_red(x,T,p));
! 2441: }
! 2442: av0 = avma; dz = dx-dy;
! 2443: #if 0 /* FIXME: to be done */
! 2444: if (OK_ULONG(p))
! 2445: { /* assume ab != 0 mod p */
! 2446: ulong pp = (ulong)p[2];
! 2447: GEN a = u_Fp_FpX(x,1, pp);
! 2448: GEN b = u_Fp_FpX(y,1, pp);
! 2449: GEN zz= u_FpX_divrem(a,b,pp,1, pr);
! 2450: if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
! 2451: {
! 2452: rem = small_to_pol(*pr,vx);
! 2453: free(*pr); *pr = rem;
! 2454: }
! 2455: z = small_to_pol(zz,vx);
! 2456: free(zz); free(a); free(b); return z;
! 2457: }
! 2458: #endif
! 2459: lead = gcmp1(lead)? NULL: gclone(FpXQ_inv(lead,T,p));
! 2460: avma = av0;
! 2461: z=cgetg(dz+3,t_POL);
! 2462: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 2463: x += 2; y += 2; z += 2;
! 2464:
! 2465: p1 = (GEN)x[dx]; av = avma;
! 2466: z[dz] = lead? lpileupto(av, FpX_res(gmul(p1,lead), T, p)): lcopy(p1);
! 2467: for (i=dx-1; i>=dy; i--)
! 2468: {
! 2469: av=avma; p1=(GEN)x[i];
! 2470: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2471: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2472: if (lead) p1 = gmul(FpX_res(p1, T, p), lead);
! 2473: tetpil=avma; z[i-dy] = lpile(av,tetpil,FpX_res(p1, T, p));
! 2474: }
! 2475: if (!pr) { if (lead) gunclone(lead); return z-2; }
! 2476:
! 2477: rem = (GEN)avma; av = (long)new_chunk(dx+3);
! 2478: for (sx=0; ; i--)
! 2479: {
! 2480: p1 = (GEN)x[i];
! 2481: for (j=0; j<=i && j<=dz; j++)
! 2482: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2483: tetpil=avma; p1 = FpX_res(p1, T, p); if (signe(p1)) { sx = 1; break; }
! 2484: if (!i) break;
! 2485: avma=av;
! 2486: }
! 2487: if (pr == ONLY_DIVIDES)
! 2488: {
! 2489: if (lead) gunclone(lead);
! 2490: if (sx) { avma=av0; return NULL; }
! 2491: avma = (long)rem; return z-2;
! 2492: }
! 2493: lrem=i+3; rem -= lrem;
! 2494: rem[0]=evaltyp(t_POL) | evallg(lrem);
! 2495: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
! 2496: p1 = gerepile((long)rem,tetpil,p1);
! 2497: rem += 2; rem[i]=(long)p1;
! 2498: for (i--; i>=0; i--)
! 2499: {
! 2500: av=avma; p1 = (GEN)x[i];
! 2501: for (j=0; j<=i && j<=dz; j++)
! 2502: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2503: tetpil=avma; rem[i]=lpile(av,tetpil, FpX_res(p1, T, p));
! 2504: }
! 2505: rem -= 2;
! 2506: if (lead) gunclone(lead);
! 2507: if (!sx) normalizepol_i(rem, lrem);
! 2508: if (pr == ONLY_REM) return gerepileupto(av0,rem);
! 2509: *pr = rem; return z-2;
! 2510: }
! 2511:
! 2512: static GEN
! 2513: FpX_gcd_long(GEN x, GEN y, GEN p)
! 2514: {
! 2515: ulong pp = (ulong)p[2], av = avma;
! 2516: GEN a,b;
! 2517:
! 2518: (void)new_chunk((lgef(x) + lgef(y)) << 2); /* scratch space */
! 2519: a = u_Fp_FpX(x,0, pp);
! 2520: if (!signe(a)) { avma = av; return FpX_red(y,p); }
! 2521: b = u_Fp_FpX(y,0, pp);
! 2522: a = u_FpX_gcd_i(a,b, pp);
! 2523: avma = av; return small_to_pol(a, varn(x));
! 2524: }
! 2525:
! 2526: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)) */
! 2527: GEN
! 2528: FpX_gcd(GEN x, GEN y, GEN p)
! 2529: {
! 2530: GEN a,b,c;
! 2531: long av0,av;
! 2532:
! 2533: if (OK_ULONG(p)) return FpX_gcd_long(x,y,p);
! 2534: av0=avma;
! 2535: a = FpX_red(x, p); av = avma;
! 2536: b = FpX_red(y, p);
! 2537: while (signe(b))
! 2538: {
! 2539: av = avma; c = FpX_res(a,b,p); a=b; b=c;
! 2540: }
! 2541: avma = av; return gerepileupto(av0, a);
! 2542: }
! 2543:
! 2544: GEN
! 2545: u_FpX_sub(GEN x, GEN y, ulong p)
! 2546: {
! 2547: long i,lz,lx = lgef(x), ly = lgef(y);
! 2548: GEN z;
! 2549:
! 2550: if (ly <= lx)
! 2551: {
! 2552: lz = lx; z = cgetg(lz,t_VECSMALL);
! 2553: for (i=2; i<ly; i++) z[i] = subssmod(x[i],y[i],p);
! 2554: for ( ; i<lx; i++) z[i] = x[i];
! 2555: }
! 2556: else
! 2557: {
! 2558: lz = ly; z = cgetg(lz,t_VECSMALL);
! 2559: for (i=2; i<lx; i++) z[i] = subssmod(x[i],y[i],p);
! 2560: for ( ; i<ly; i++) z[i] = y[i]? p - y[i]: y[i];
! 2561: }
! 2562: z[1]=0; return u_normalizepol(z, lz);
! 2563: }
! 2564:
! 2565: /* list of u_FpX in gptr, return then as GEN */
! 2566: static void
! 2567: u_gerepilemany(long av, GEN* gptr[], long n, long v)
! 2568: {
! 2569: GEN *l = (GEN*)gpmalloc(n*sizeof(GEN));
! 2570: long i;
! 2571:
! 2572: /* copy objects */
! 2573: for (i=0; i<n; i++) l[i] = gclone(*(gptr[i]));
! 2574: avma = av;
! 2575: /* copy them back, kill clones */
! 2576: for (--i; i>=0; i--)
! 2577: {
! 2578: *(gptr[i]) = small_to_pol(l[i],v);
! 2579: gunclone(l[i]);
! 2580: }
! 2581: free(l);
! 2582: }
! 2583:
! 2584: static GEN
! 2585: u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
! 2586: {
! 2587: GEN q,z,u,v, x = a, y = b;
! 2588:
! 2589: u = u_zeropol(0);
! 2590: v= u_scalarpol(1, 0); /* v = 1 */
! 2591: while (signe(y))
! 2592: {
! 2593: q = u_FpX_divrem(x,y,p, 0,&z);
! 2594: x = y; y = z; /* (x,y) = (y, x - q y) */
! 2595: z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
! 2596: u = v; v = z; /* (u,v) = (v, u - q v) */
! 2597: }
! 2598: z = u_FpX_sub(x, u_FpX_mul(b,u,p), p);
! 2599: z = u_FpX_div(z,a,p);
! 2600: *ptu = z;
! 2601: *ptv = u; return x;
! 2602: }
! 2603:
! 2604: static GEN
! 2605: FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
! 2606: {
! 2607: ulong pp = (ulong)p[2];
! 2608: long av = avma;
! 2609: GEN a, b, d;
! 2610:
! 2611: a = u_Fp_FpX(x,0, pp);
! 2612: b = u_Fp_FpX(y,0, pp);
! 2613: d = u_FpX_extgcd(a,b, pp, ptu,ptv);
! 2614: {
! 2615: GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv;
! 2616: u_gerepilemany(av, gptr, 3, varn(x));
! 2617: }
! 2618: return d;
! 2619: }
! 2620:
! 2621: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
! 2622: * ux + vy = gcd (mod p) */
! 2623: GEN
! 2624: FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
! 2625: {
! 2626: GEN a,b,q,r,u,v,d,d1,v1;
! 2627: long ltop,lbot;
! 2628:
! 2629: if (OK_ULONG(p)) return FpX_extgcd_long(x,y,p,ptu,ptv);
! 2630: ltop=avma;
! 2631: a = FpX_red(x, p);
! 2632: b = FpX_red(y, p);
! 2633: d = a; d1 = b; v = gzero; v1 = gun;
! 2634: while (signe(d1))
! 2635: {
! 2636: q = FpX_divres(d,d1,p, &r);
! 2637: v = gadd(v, gneg_i(gmul(q,v1)));
! 2638: v = FpX_red(v,p);
! 2639: u=v; v=v1; v1=u;
! 2640: u=r; d=d1; d1=u;
! 2641: }
! 2642: u = gadd(d, gneg_i(gmul(b,v)));
! 2643: u = FpX_red(u, p);
! 2644: lbot = avma;
! 2645: u = FpX_div(u,a,p);
! 2646: d = gcopy(d);
! 2647: v = gcopy(v);
! 2648: {
! 2649: GEN *gptr[3]; gptr[0] = &d; gptr[1] = &u; gptr[2] = &v;
! 2650: gerepilemanysp(ltop,lbot,gptr,3);
! 2651: }
! 2652: *ptu = u; *ptv = v; return d;
! 2653: }
! 2654:
! 2655: /* x and y in Z[Y][X], return lift(gcd(x mod T,p, y mod T,p)). Set u and v st
! 2656: * ux + vy = gcd (mod T,p) */
! 2657: GEN
! 2658: FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
! 2659: {
! 2660: GEN a,b,q,r,u,v,d,d1,v1;
! 2661: long ltop,lbot;
! 2662:
! 2663: #if 0 /* FIXME To be done...*/
! 2664: if (OK_ULONG(p)) return FpXQX_extgcd_long(x,y,T,p,ptu,ptv);
! 2665: #endif
! 2666: if (!T) return FpX_extgcd(x,y,p,ptu,ptv);
! 2667: ltop=avma;
! 2668: a = FpXQX_red(x, T, p);
! 2669: b = FpXQX_red(y, T, p);
! 2670: d = a; d1 = b; v = gzero; v1 = gun;
! 2671: while (signe(d1))
! 2672: {
! 2673: q = FpXQX_divres(d,d1,T,p, &r);
! 2674: v = gadd(v, gneg_i(gmul(q,v1)));
! 2675: v = FpXQX_red(v,T,p);
! 2676: u=v; v=v1; v1=u;
! 2677: u=r; d=d1; d1=u;
! 2678: }
! 2679: u = gadd(d, gneg_i(gmul(b,v)));
! 2680: u = FpXQX_red(u,T, p);
! 2681: lbot = avma;
! 2682: u = FpXQX_divres(u,a,T,p,NULL);
! 2683: d = gcopy(d);
! 2684: v = gcopy(v);
! 2685: {
! 2686: GEN *gptr[3]; gptr[0] = &d; gptr[1] = &u; gptr[2] = &v;
! 2687: gerepilemanysp(ltop,lbot,gptr,3);
! 2688: }
! 2689: *ptu = u; *ptv = v; return d;
! 2690: }
! 2691:
! 2692: extern GEN caractducos(GEN p, GEN x, int v);
! 2693:
! 2694: GEN
! 2695: FpXQ_charpoly(GEN x, GEN T, GEN p)
! 2696: {
! 2697: ulong ltop=avma;
! 2698: GEN R=lift(caractducos(FpX(T,p),FpX(x,p),varn(T)));
! 2699: return gerepileupto(ltop,R);
! 2700: }
! 2701:
! 2702: GEN
! 2703: FpXQ_minpoly(GEN x, GEN T, GEN p)
! 2704: {
! 2705: ulong ltop=avma;
! 2706: GEN R=FpXQ_charpoly(x, T, p);
! 2707: GEN G=FpX_gcd(R,derivpol(R),p);
! 2708: G=FpX_normalize(G,p);
! 2709: G=FpX_div(R,G,p);
! 2710: return gerepileupto(ltop,G);
! 2711: }
! 2712:
! 2713: /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
! 2714: static GEN
! 2715: u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
! 2716: {
! 2717: ulong av = avma, d, amod = umodiu(a,p);
! 2718: GEN ax;
! 2719:
! 2720: if (b == amod) return NULL;
! 2721: d = (b > amod)? b - amod: p - (amod - b); /* (b - a) mod p */
! 2722: (void)new_chunk(lgefint(pq)<<1); /* HACK */
! 2723: ax = mului(mulssmod(d,qinv,p), q); /* d mod p, 0 mod q */
! 2724: avma = av; return addii(a, ax); /* in ]-q, pq[ assuming a in -]-q,q[ */
! 2725: }
! 2726:
! 2727: /* centerlift(u mod p) */
! 2728: GEN
! 2729: u_center(ulong u, ulong p, ulong ps2)
! 2730: {
! 2731: return stoi((long) (u > ps2)? u - p: u);
! 2732: }
! 2733:
! 2734: GEN
! 2735: ZX_init_CRT(GEN Hp, ulong p, long v)
! 2736: {
! 2737: long i, l = lgef(Hp), lim = (long)(p>>1);
! 2738: GEN H = cgetg(l, t_POL);
! 2739: H[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
! 2740: for (i=2; i<l; i++)
! 2741: H[i] = (long)u_center(Hp[i], p, lim);
! 2742: return H;
! 2743: }
! 2744:
! 2745: /* assume lg(Hp) > 1 */
! 2746: GEN
! 2747: ZM_init_CRT(GEN Hp, ulong p)
! 2748: {
! 2749: long i,j, m = lg(Hp[1]), l = lg(Hp), lim = (long)(p>>1);
! 2750: GEN c,cp,H = cgetg(l, t_MAT);
! 2751: for (j=1; j<l; j++)
! 2752: {
! 2753: cp = (GEN)Hp[j];
! 2754: c = cgetg(m, t_COL);
! 2755: H[j] = (long)c;
! 2756: for (i=1; i<l; i++) c[i] = (long)u_center(cp[i],p, lim);
! 2757: }
! 2758: return H;
! 2759: }
! 2760:
! 2761: int
! 2762: Z_incremental_CRT(GEN *H, ulong Hp, GEN q, GEN qp, ulong p)
! 2763: {
! 2764: GEN h, lim = shifti(qp,-1);
! 2765: ulong qinv = u_invmod(umodiu(q,p), p);
! 2766: int stable = 1;
! 2767: h = u_chrem_coprime(*H,Hp,q,p,qinv,qp);
! 2768: if (h)
! 2769: {
! 2770: if (cmpii(h,lim) > 0) h = subii(h,qp);
! 2771: *H = h; stable = 0;
! 2772: }
! 2773: return stable;
! 2774: }
! 2775:
! 2776: int
! 2777: ZX_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p)
! 2778: {
! 2779: GEN h, lim = shifti(qp,-1);
! 2780: ulong qinv = u_invmod(umodiu(q,p), p);
! 2781: long i, l = lgef(H), lp = lgef(Hp);
! 2782: int stable = 1;
! 2783: for (i=2; i<lp; i++)
! 2784: {
! 2785: h = u_chrem_coprime((GEN)H[i],Hp[i],q,p,qinv,qp);
! 2786: if (h)
! 2787: {
! 2788: if (cmpii(h,lim) > 0) h = subii(h,qp);
! 2789: H[i] = (long)h; stable = 0;
! 2790: }
! 2791: }
! 2792: for ( ; i<l; i++)
! 2793: {
! 2794: h = u_chrem_coprime((GEN)H[i], 0,q,p,qinv,qp);
! 2795: if (h)
! 2796: {
! 2797: if (cmpii(h,lim) > 0) h = subii(h,qp);
! 2798: H[i] = (long)h; stable = 0;
! 2799: }
! 2800: }
! 2801: return stable;
! 2802: }
! 2803:
! 2804: int
! 2805: ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p)
! 2806: {
! 2807: GEN h, lim = shifti(qp,-1);
! 2808: ulong qinv = u_invmod(umodiu(q,p), p);
! 2809: long i,j, l = lg(H), m = lg(H[1]);
! 2810: int stable = 1;
! 2811: for (j=1; j<l; j++)
! 2812: for (i=1; i<m; i++)
! 2813: {
! 2814: h = u_chrem_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
! 2815: if (h)
! 2816: {
! 2817: if (cmpii(h,lim) > 0) h = subii(h,qp);
! 2818: coeff(H,i,j) = (long)h; stable = 0;
! 2819: }
! 2820: }
! 2821: return stable;
! 2822: }
! 2823:
! 2824: /* returns a polynomial in variable v, whose coeffs correspond to the digits
! 2825: * of m (in base p) */
! 2826: GEN
! 2827: stopoly(long m, long p, long v)
! 2828: {
! 2829: GEN y = cgetg(BITS_IN_LONG + 2, t_POL);
! 2830: long l=2;
! 2831:
! 2832: do { y[l++] = lstoi(m%p); m=m/p; } while (m);
! 2833: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
! 2834: return y;
! 2835: }
! 2836:
! 2837: GEN
! 2838: stopoly_gen(GEN m, GEN p, long v)
! 2839: {
! 2840: GEN y = cgetg(bit_accuracy(lgefint(m))+2, t_POL);
! 2841: long l=2;
! 2842:
! 2843: do { y[l++] = lmodii(m,p); m=divii(m,p); } while (signe(m));
! 2844: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
! 2845: return y;
! 2846: }
! 2847:
! 2848: /* modular power */
! 2849: ulong
! 2850: powuumod(ulong x, ulong n0, ulong p)
! 2851: {
! 2852: ulong y, z, n;
! 2853: if (n0 <= 2)
! 2854: { /* frequent special cases */
! 2855: if (n0 == 2) return mulssmod(x,x,p);
! 2856: if (n0 == 1) return x;
! 2857: if (n0 == 0) return 1;
! 2858: }
! 2859: y = 1; z = x; n = n0;
! 2860: for(;;)
! 2861: {
! 2862: if (n&1) y = mulssmod(y,z,p);
! 2863: n>>=1; if (!n) return y;
! 2864: z = mulssmod(z,z,p);
! 2865: }
! 2866: }
! 2867:
! 2868: /* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0 */
! 2869: GEN
! 2870: u_FpX_rem(GEN x, GEN y, ulong p)
! 2871: {
! 2872: GEN z, c;
! 2873: long dx,dy,dz,i,j;
! 2874: ulong p1,inv;
! 2875:
! 2876: dy = degpol(y); if (!dy) return u_zeropol(0);
! 2877: dx = degpol(x);
! 2878: dz = dx-dy; if (dz < 0) return u_copy(x, 0);
! 2879: x += 2;
! 2880: y += 2;
! 2881: z = u_allocpol(dz, 1) + 2;
! 2882: inv = y[dy];
! 2883: if (inv != 1UL) inv = u_invmod(inv,p);
! 2884:
! 2885: c = u_allocpol(dy,0) + 2;
! 2886: if (u_OK_ULONG(p))
! 2887: {
! 2888: z[dz] = (inv*x[dx]) % p;
! 2889: for (i=dx-1; i>=dy; --i)
! 2890: {
! 2891: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
! 2892: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2893: {
! 2894: p1 += z[j]*y[i-j];
! 2895: if (p1 & HIGHBIT) p1 = p1 % p;
! 2896: }
! 2897: p1 %= p;
! 2898: z[i-dy] = p1? ((p - p1)*inv) % p: 0;
! 2899: }
! 2900: for (i=0; i<dy; i++)
! 2901: {
! 2902: p1 = z[0]*y[i];
! 2903: for (j=1; j<=i && j<=dz; j++)
! 2904: {
! 2905: p1 += z[j]*y[i-j];
! 2906: if (p1 & HIGHBIT) p1 = p1 % p;
! 2907: }
! 2908: c[i] = subssmod(x[i], p1%p, p);
! 2909: }
! 2910: }
! 2911: else
! 2912: {
! 2913: z[dz] = mulssmod(inv, x[dx], p);
! 2914: for (i=dx-1; i>=dy; --i)
! 2915: {
! 2916: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
! 2917: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2918: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
! 2919: z[i-dy] = p1? mulssmod(p - p1, inv, p): 0;
! 2920: }
! 2921: for (i=0; i<dy; i++)
! 2922: {
! 2923: p1 = mulssmod(z[0],y[i],p);
! 2924: for (j=1; j<=i && j<=dz; j++)
! 2925: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
! 2926: c[i] = subssmod(x[i], p1, p);
! 2927: }
! 2928: }
! 2929: i = dy-1; while (i>=0 && !c[i]) i--;
! 2930: free(z-2); return u_normalizepol(c-2, i+3);
! 2931: }
! 2932:
! 2933: ulong
! 2934: u_FpX_resultant(GEN a, GEN b, ulong p)
! 2935: {
! 2936: long da,db,dc,cnt;
! 2937: ulong lb,av, res = 1UL;
! 2938: GEN c;
! 2939:
! 2940: if (!signe(a) || !signe(b)) return 0;
! 2941: da = degpol(a);
! 2942: db = degpol(b);
! 2943: if (db > da)
! 2944: {
! 2945: swapspec(a,b, da,db);
! 2946: if (da & db & 1) res = p-res;
! 2947: }
! 2948: if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
! 2949: cnt = 0; av = avma;
! 2950: while (db)
! 2951: {
! 2952: lb = b[db+2];
! 2953: c = u_FpX_rem(a,b, p);
! 2954: a = b; b = c; dc = degpol(c);
! 2955: if (dc < 0) { avma = av; return 0; }
! 2956:
! 2957: if (da & db & 1) res = p - res;
! 2958: if (lb != 1) res = mulssmod(res, powuumod(lb, da - dc, p), p);
! 2959: if (++cnt == 4) { cnt = 0; avma = av; }
! 2960: da = db; /* = degpol(a) */
! 2961: db = dc; /* = degpol(b) */
! 2962: }
! 2963: avma = av; return mulssmod(res, powuumod(b[2], da, p), p);
! 2964: }
! 2965:
! 2966: /* If resultant is 0, *ptU and *ptU are not set */
! 2967: ulong
! 2968: u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
! 2969: {
! 2970: GEN z,q,u,v, x = a, y = b;
! 2971: ulong lb, av = avma, res = 1UL;
! 2972: long dx,dy,dz;
! 2973:
! 2974: if (!signe(x) || !signe(y)) return 0;
! 2975: dx = degpol(x);
! 2976: dy = degpol(y);
! 2977: if (dy > dx)
! 2978: {
! 2979: swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
! 2980: a = x; b = y;
! 2981: if (dx & dy & 1) res = p-res;
! 2982: }
! 2983: u = u_zeropol(0);
! 2984: v = u_scalarpol(1, 0); /* v = 1 */
! 2985: while (dy)
! 2986: { /* b u = x (a), b v = y (a) */
! 2987: lb = y[dy+2];
! 2988: q = u_FpX_divrem(x,y, p, 0, &z);
! 2989: x = y; y = z; /* (x,y) = (y, x - q y) */
! 2990: dz = degpol(z); if (dz < 0) { avma = av; return 0; }
! 2991: z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
! 2992: u = v; v = z; /* (u,v) = (v, u - q v) */
! 2993:
! 2994: if (dx & dy & 1) res = p - res;
! 2995: if (lb != 1) res = mulssmod(res, powuumod(lb, dx-dz, p), p);
! 2996: dx = dy; /* = degpol(x) */
! 2997: dy = dz; /* = degpol(y) */
! 2998: }
! 2999: res = mulssmod(res, powuumod(y[2], dx, p), p);
! 3000: lb = mulssmod(res, u_invmod(y[2],p), p);
! 3001: v = gerepileupto(av, u_FpX_Fp_mul(v, lb, p, 0));
! 3002: av = avma;
! 3003: u = u_FpX_sub(u_scalarpol(res,0), u_FpX_mul(b,v,p), p);
! 3004: u = gerepileupto(av, u_FpX_div(u,a,p)); /* = (res - b v) / a */
! 3005: *ptU = u;
! 3006: *ptV = v; return res;
! 3007: }
! 3008:
! 3009: /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with "generic"
! 3010: * degree sequence as given by dglist, set *Ci and return resultant(a,b) */
! 3011: ulong
! 3012: u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
! 3013: {
! 3014: long da,db,dc,cnt,ind;
! 3015: ulong lb,av,cx = 1, res = 1UL;
! 3016: GEN c;
! 3017:
! 3018: if (C0) { *C0 = 1; *C1 = 0; }
! 3019: if (!signe(a) || !signe(b)) return 0;
! 3020: da = degpol(a);
! 3021: db = degpol(b);
! 3022: if (db > da)
! 3023: {
! 3024: swapspec(a,b, da,db);
! 3025: if (da & db & 1) res = p-res;
! 3026: }
! 3027: /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
! 3028: if (!da) return 1;
! 3029: cnt = ind = 0; av = avma;
! 3030: while (db)
! 3031: {
! 3032: lb = b[db+2];
! 3033: c = u_FpX_rem(a,b, p);
! 3034: a = b; b = c; dc = degpol(c);
! 3035: if (dc < 0) { avma = av; return 0; }
! 3036:
! 3037: ind++;
! 3038: if (C0)
! 3039: { /* check that Euclidean remainder sequence doesn't degenerate */
! 3040: if (dc != dglist[ind]) { avma = av; return 0; }
! 3041: /* update resultant */
! 3042: if (da & db & 1) res = p-res;
! 3043: if (lb != 1)
! 3044: {
! 3045: ulong t = powuumod(lb, da - dc, p);
! 3046: res = mulssmod(res, t, p);
! 3047: if (dc) cx = mulssmod(cx, t, p);
! 3048: }
! 3049: }
! 3050: else
! 3051: {
! 3052: if (dc > dglist[ind]) dglist[ind] = dc;
! 3053: }
! 3054: if (++cnt == 4) { cnt = 0; avma = av; }
! 3055: da = db; /* = degpol(a) */
! 3056: db = dc; /* = degpol(b) */
! 3057: }
! 3058: if (!C0)
! 3059: {
! 3060: if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
! 3061: return 0;
! 3062: }
! 3063:
! 3064: if (da == 1) /* last non-constant polynomial has degree 1 */
! 3065: {
! 3066: *C0 = mulssmod(cx, a[2], p);
! 3067: *C1 = mulssmod(cx, a[3], p);
! 3068: lb = b[2];
! 3069: } else lb = powuumod(b[2], da, p);
! 3070: avma = av; return mulssmod(res, lb, p);
! 3071: }
! 3072:
! 3073: static ulong global_pp;
! 3074: static GEN
! 3075: _u_FpX_mul(GEN a, GEN b)
! 3076: {
! 3077: return u_FpX_mul(a,b, global_pp);
! 3078: }
! 3079:
! 3080: /* compute prod (x - a[i]) */
! 3081: GEN
! 3082: u_FpV_roots_to_pol(GEN a, ulong p)
! 3083: {
! 3084: long i,k,lx = lg(a);
! 3085: GEN p1,p2;
! 3086: if (lx == 1) return polun[0];
! 3087: p1 = cgetg(lx, t_VEC); global_pp = p;
! 3088: for (k=1,i=1; i<lx-1; i+=2)
! 3089: {
! 3090: p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2;
! 3091: p2[2] = mulssmod(a[i], a[i+1], p);
! 3092: p2[3] = (p<<1) - (a[i] + a[i+1]);
! 3093: p2[4] = 1; p2[1] = evallgef(5);
! 3094: }
! 3095: if (i < lx)
! 3096: {
! 3097: p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
! 3098: p2[1] = evallgef(4);
! 3099: p2[2] = p - a[i];
! 3100: p2[3] = 1;
! 3101: }
! 3102: setlg(p1, k); return divide_conquer_prod(p1, _u_FpX_mul);
! 3103: }
! 3104:
! 3105:
! 3106: /* u P(X) + v P(-X) */
! 3107: static GEN
! 3108: pol_comp(GEN P, GEN u, GEN v)
! 3109: {
! 3110: long i, l = lgef(P);
! 3111: GEN y = cgetg(l, t_POL);
! 3112: for (i=2; i<l; i++)
! 3113: {
! 3114: GEN t = (GEN)P[i];
! 3115: y[i] = gcmp0(t)? zero:
! 3116: (i&1)? lmul(t, gsub(u,v)) /* odd degree */
! 3117: : lmul(t, gadd(u,v));/* even degree */
! 3118: }
! 3119: y[1] = P[1]; return normalizepol_i(y,l);
! 3120: }
! 3121:
! 3122: static GEN
! 3123: u_pol_comp(GEN P, ulong u, ulong v, ulong p)
! 3124: {
! 3125: long i, l = lgef(P);
! 3126: GEN y = u_allocpol(l-3, 0);
! 3127: for (i=2; i<l; i++)
! 3128: {
! 3129: ulong t = P[i];
! 3130: y[i] = (t == 0)? 0:
! 3131: (i&1)? mulssmod(t, u + (p - v), p)
! 3132: : mulssmod(t, u + v, p);
! 3133: }
! 3134: return u_normalizepol(y,l);
! 3135: }
! 3136:
! 3137: GEN roots_to_pol(GEN a, long v);
! 3138:
! 3139: GEN
! 3140: polint_triv(GEN xa, GEN ya)
! 3141: {
! 3142: GEN P = NULL, Q = roots_to_pol(xa,0);
! 3143: long i, n = lg(xa), av = avma, lim = stack_lim(av,2);
! 3144: for (i=1; i<n; i++)
! 3145: {
! 3146: GEN T,dP;
! 3147: if (gcmp0((GEN)ya[i])) continue;
! 3148: T = gdeuc(Q, gsub(polx[0], (GEN)xa[i]));
! 3149: if (i < n-1 && absi_equal((GEN)xa[i], (GEN)xa[i+1]))
! 3150: { /* x_i = -x_{i+1} */
! 3151: T = gdiv(T, poleval(T, (GEN)xa[i]));
! 3152: dP = pol_comp(T, (GEN)ya[i], (GEN)ya[i+1]);
! 3153: i++;
! 3154: }
! 3155: else
! 3156: {
! 3157: dP = gmul((GEN)ya[i], T);
! 3158: dP = gdiv(dP, poleval(T,(GEN)xa[i]));
! 3159: }
! 3160: P = P? gadd(P, dP): dP;
! 3161: if (low_stack(lim,stack_lim(av,2)))
! 3162: {
! 3163: if (DEBUGMEM>1) err(warnmem,"polint_triv2 (i = %ld)",i);
! 3164: P = gerepileupto(av, P);
! 3165: }
! 3166: }
! 3167: return P? P: zeropol(0);
! 3168: }
! 3169:
! 3170: ulong
! 3171: u_FpX_eval(GEN x, ulong y, ulong p)
! 3172: {
! 3173: ulong p1,r;
! 3174: long i,j;
! 3175: i=lgef(x)-1;
! 3176: if (i<=2)
! 3177: return (i==2)? x[2]: 0;
! 3178: p1 = x[i];
! 3179: /* specific attention to sparse polynomials (see poleval)*/
! 3180: if (u_OK_ULONG(p))
! 3181: {
! 3182: for (i--; i>=2; i=j-1)
! 3183: {
! 3184: for (j=i; !x[j]; j--)
! 3185: if (j==2)
! 3186: {
! 3187: if (i != j) y = powuumod(y, i-j+1, p);
! 3188: return (p1 * y) % p;
! 3189: }
! 3190: r = (i==j)? y: powuumod(y, i-j+1, p);
! 3191: p1 = ((p1*r) + x[j]) % p;
! 3192: }
! 3193: }
! 3194: else
! 3195: {
! 3196: for (i--; i>=2; i=j-1)
! 3197: {
! 3198: for (j=i; !x[j]; j--)
! 3199: if (j==2)
! 3200: {
! 3201: if (i != j) y = powuumod(y, i-j+1, p);
! 3202: return mulssmod(p1, y, p);
! 3203: }
! 3204: r = (i==j)? y: powuumod(y, i-j+1, p);
! 3205: p1 = addssmod(x[j], mulssmod(p1,r,p), p);
! 3206: }
! 3207: }
! 3208: return p1;
! 3209: }
! 3210:
! 3211: static GEN
! 3212: u_FpX_div_by_X_x(GEN a, ulong x, ulong p)
! 3213: {
! 3214: long l = lgef(a), i;
! 3215: GEN z = u_allocpol(l-4, 0), a0, z0;
! 3216: a0 = a + l-1;
! 3217: z0 = z + l-2; *z0 = *a0--;
! 3218: if (u_OK_ULONG(p))
! 3219: {
! 3220: for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
! 3221: {
! 3222: ulong t = (*a0-- + x * *z0--) % p;
! 3223: *z0 = t;
! 3224: }
! 3225: }
! 3226: else
! 3227: {
! 3228: for (i=l-3; i>1; i--)
! 3229: {
! 3230: ulong t = addssmod(*a0--, mulssmod(x, *z0--, p), p);
! 3231: *z0 = t;
! 3232: }
! 3233: }
! 3234: return z;
! 3235: }
! 3236:
! 3237: /* xa, ya = t_VECSMALL */
! 3238: GEN
! 3239: u_FpV_polint(GEN xa, GEN ya, ulong p)
! 3240: {
! 3241: GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p);
! 3242: long i, n = lg(xa);
! 3243: ulong av, inv;
! 3244: av = avma; (void)new_chunk(n+3); /* storage space for P */
! 3245: for (i=1; i<n; i++)
! 3246: {
! 3247: if (!ya[i]) continue;
! 3248: T = u_FpX_div_by_X_x(Q, xa[i], p);
! 3249: inv = u_invmod(u_FpX_eval(T,xa[i], p), p);
! 3250: if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p)
! 3251: {
! 3252: dP = u_pol_comp(T, mulssmod(ya[i],inv,p), mulssmod(ya[i+1],inv,p), p);
! 3253: i++; /* x_i = -x_{i+1} */
! 3254: }
! 3255: else
! 3256: dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p, 0);
! 3257: if (P) avma = av;
! 3258: P = P? u_FpX_add(P, dP, p): dP;
! 3259: }
! 3260: return P? P: u_zeropol(0);
! 3261: }
! 3262:
! 3263: static void
! 3264: u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p)
! 3265: {
! 3266: GEN T,Q = u_FpV_roots_to_pol(xa, p);
! 3267: GEN dP = NULL, P = NULL;
! 3268: GEN dP0 = NULL, P0= NULL;
! 3269: GEN dP1 = NULL, P1= NULL;
! 3270: long i, n = lg(xa);
! 3271: ulong inv;
! 3272: for (i=1; i<n; i++)
! 3273: {
! 3274: T = u_FpX_div_by_X_x(Q, xa[i], p);
! 3275: inv = u_invmod(u_FpX_eval(T,xa[i], p), p);
! 3276:
! 3277: #if 0
! 3278: if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p)
! 3279: {
! 3280: if (ya[i])
! 3281: dP = u_pol_comp(T, mulssmod(ya[i],inv,p), mulssmod(ya[i+1],inv,p), p);
! 3282: if (C0[i])
! 3283: dP0= u_pol_comp(T, mulssmod(C0[i],inv,p), mulssmod(C0[i+1],inv,p), p);
! 3284: if (C1[i])
! 3285: dP1= u_pol_comp(T, mulssmod(C1[i],inv,p), mulssmod(C1[i+1],inv,p), p);
! 3286: i++; /* x_i = -x_{i+1} */
! 3287: }
! 3288: else
! 3289: #endif
! 3290: {
! 3291: if (ya[i])
! 3292: {
! 3293: dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p, 0);
! 3294: P = P ? u_FpX_add(P , dP , p): dP;
! 3295: }
! 3296: if (C0[i])
! 3297: {
! 3298: dP0= u_FpX_Fp_mul(T, mulssmod(C0[i],inv,p), p, 0);
! 3299: P0= P0? u_FpX_add(P0, dP0, p): dP0;
! 3300: }
! 3301: if (C1[i])
! 3302: {
! 3303: dP1= u_FpX_Fp_mul(T, mulssmod(C1[i],inv,p), p, 0);
! 3304: P1= P1? u_FpX_add(P1, dP1, p): dP1;
! 3305: }
! 3306: }
! 3307: }
! 3308: ya[1] = (long) (P ? P : u_zeropol(0));
! 3309: C0[1] = (long) (P0? P0: u_zeropol(0));
! 3310: C1[1] = (long) (P1? P1: u_zeropol(0));
! 3311: }
! 3312:
! 3313: /* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
! 3314: * Return 0 in case of degree drop. */
! 3315: static GEN
! 3316: vec_u_FpX_eval(GEN b, ulong x, ulong p)
! 3317: {
! 3318: GEN z;
! 3319: long i, lb = lgef(b);
! 3320: ulong leadz = u_FpX_eval((GEN)b[lb-1], x, p);
! 3321: if (!leadz) return u_zeropol(0);
! 3322:
! 3323: z = u_allocpol(lb-3, 0);
! 3324: for (i=2; i<lb-1; i++)
! 3325: z[i] = u_FpX_eval((GEN)b[i], x, p);
! 3326: z[i] = leadz; return z;
! 3327: }
! 3328:
! 3329: /* as above, but don't care about degree drop */
! 3330: static GEN
! 3331: vec_u_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop)
! 3332: {
! 3333: GEN z;
! 3334: long i, lb = lgef(b);
! 3335: z = u_allocpol(lb-3, 0);
! 3336: for (i=2; i<lb; i++)
! 3337: z[i] = u_FpX_eval((GEN)b[i], x, p);
! 3338: z = u_normalizepol(z, lb);
! 3339: *drop = lb - lgef(z);
! 3340: return z;
! 3341: }
! 3342:
! 3343: /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
! 3344: * bound = N_2(A)^degpol B N_2(B)^degpol(A), where N_2(A) = sqrt(sum (N_1(Ai))^2)
! 3345: * Return B such that bound < 2^B */
! 3346: ulong
! 3347: ZY_ZXY_ResBound(GEN A, GEN B)
! 3348: {
! 3349: ulong av = avma;
! 3350: GEN a = gzero, b = gzero, run = realun(DEFAULTPREC);
! 3351: long i , lA = lgef(A), lB = lgef(B);
! 3352: for (i=2; i<lA; i++) a = addii(a, sqri((GEN)A[i]));
! 3353: for (i=2; i<lB; i++)
! 3354: {
! 3355: GEN t = (GEN)B[i];
! 3356: if (typ(t) == t_POL) t = gnorml1(t, 0);
! 3357: b = addii(b, sqri(t));
! 3358: }
! 3359: b = gmul(gpowgs(mulir(a,run), degpol(B)), gpowgs(mulir(b,run), degpol(A)));
! 3360: avma = av; return 1 + (gexpo(b)>>1);
! 3361: }
! 3362:
! 3363: /* 0, 1, -1, 2, -2, ... */
! 3364: #define next_lambda(a) (a>0 ? -a : 1-a)
! 3365:
! 3366: /* If lambda = NULL, assume A in Z[Y], B in Q[Y][X], return Res_Y(A,B)
! 3367: * Otherwise, find a small lambda (start from *lambda, use the sequence above)
! 3368: * such that R(X) = Res_Y(A, B(X + lambda Y)) is squarefree, reset *lambda to
! 3369: * the chosen value and return R
! 3370: *
! 3371: * If LERS is non-NULL, set it to the last non-constant polynomial in the PRS */
! 3372: GEN
! 3373: ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LERS)
! 3374: {
! 3375: int checksqfree = lambda? 1: 0, delete = 0, first = 1, stable;
! 3376: ulong av = avma, av2, lim, bound;
! 3377: long i,n, lb, dres = degpol(A)*degpol(B0), nmax = (dres+1)>>1;
! 3378: long vX = varn(B0), vY = varn(A); /* assume vX < vY */
! 3379: GEN x,y,dglist,cB,B,q,a,b,ev,H,H0,H1,Hp,H0p,H1p,C0,C1;
! 3380: byteptr d = diffptr + 3000;
! 3381: ulong p = 27449; /* p = prime(3000) */
! 3382:
! 3383: dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
! 3384: if (LERS)
! 3385: {
! 3386: if (!lambda) err(talker,"ZY_ZXY_resultant_all: LERS needs lambda");
! 3387: C0 = cgetg(dres+2, t_VECSMALL);
! 3388: C1 = cgetg(dres+2, t_VECSMALL);
! 3389: dglist = cgetg(dres+1, t_VECSMALL);
! 3390: }
! 3391: x = cgetg(dres+2, t_VECSMALL);
! 3392: y = cgetg(dres+2, t_VECSMALL);
! 3393: if (vY == MAXVARN)
! 3394: {
! 3395: vY = fetch_var(); delete = 1;
! 3396: B0 = gsubst(B0, MAXVARN, polx[vY]);
! 3397: A = dummycopy(A); setvarn(A, vY);
! 3398: }
! 3399: cB = content(B0);
! 3400: if (typ(cB) == t_POL) cB = content(cB);
! 3401: if (gcmp1(cB)) cB = NULL; else B0 = gdiv(B0, cB);
! 3402:
! 3403: if (lambda)
! 3404: B = poleval(B0, gadd(polx[MAXVARN], gmulsg(*lambda, polx[vY])));
! 3405: else
! 3406: B = poleval(B0, polx[MAXVARN]);
! 3407: av2 = avma; lim = stack_lim(av,2);
! 3408:
! 3409: INIT:
! 3410: if (first) first = 0;
! 3411: else
! 3412: { /* lambda != NULL */
! 3413: avma = av2; *lambda = next_lambda(*lambda);
! 3414: if (DEBUGLEVEL>4) fprintferr("Starting with lambda = %ld\n",*lambda);
! 3415: B = poleval(B0, gadd(polx[MAXVARN], gmulsg(*lambda, polx[vY])));
! 3416: av2 = avma;
! 3417: }
! 3418:
! 3419: if (degpol(A) <= 3)
! 3420: { /* sub-resultant faster for small degrees */
! 3421: if (LERS)
! 3422: {
! 3423: H = subresall(A,B,&q);
! 3424: if (typ(q) != t_POL || degpol(q)!=1 || !ZX_is_squarefree(H)) goto INIT;
! 3425: H0 = (GEN)q[2]; if (typ(H0) == t_POL) setvarn(H0,vX);
! 3426: H1 = (GEN)q[3]; if (typ(H1) == t_POL) setvarn(H1,vX);
! 3427: }
! 3428: else
! 3429: {
! 3430: H = subres(A,B);
! 3431: if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
! 3432: }
! 3433: goto END;
! 3434: }
! 3435:
! 3436: H = H0 = H1 = NULL;
! 3437: lb = lgef(B); b = u_allocpol(degpol(B), 0);
! 3438: bound = ZY_ZXY_ResBound(A,B);
! 3439: if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound);
! 3440:
! 3441: for(;;)
! 3442: {
! 3443: p += *d++; if (!*d) err(primer1);
! 3444:
! 3445: a = u_Fp_FpX(A, 0, p);
! 3446: for (i=2; i<lb; i++)
! 3447: b[i] = (long)u_Fp_FpX((GEN)B[i], 0, p);
! 3448: if (LERS)
! 3449: {
! 3450: if (!b[lb-1] || degpol(a) < degpol(A)) continue; /* p | lc(A)lc(B) */
! 3451: if (checksqfree)
! 3452: { /* find degree list for generic Euclidean Remainder Sequence */
! 3453: int goal = min(degpol(a), degpol(b)); /* longest possible */
! 3454: for (n=1; n <= goal; n++) dglist[n] = 0;
! 3455: setlg(dglist, 1);
! 3456: for (n=0; n <= dres; n++)
! 3457: {
! 3458: ev = vec_u_FpX_eval(b, n, p);
! 3459: (void)u_FpX_resultant_all(a, ev, NULL, NULL, dglist, p);
! 3460: if (lg(dglist)-1 == goal) break;
! 3461: }
! 3462: /* last pol in ERS has degree > 1 ? */
! 3463: goal = lg(dglist)-1;
! 3464: if (degpol(B) == 1) { if (!goal) goto INIT; }
! 3465: else
! 3466: {
! 3467: if (goal <= 1) goto INIT;
! 3468: if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
! 3469: }
! 3470: if (DEBUGLEVEL>4)
! 3471: fprintferr("Degree list for ERS (trials: %ld) = %Z\n",n+1,dglist);
! 3472: }
! 3473:
! 3474: for (i=0,n = 0; i <= dres; n++)
! 3475: {
! 3476: i++; ev = vec_u_FpX_eval(b, n, p);
! 3477: x[i] = n;
! 3478: y[i] = u_FpX_resultant_all(a, ev, C0+i, C1+i, dglist, p);
! 3479: if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
! 3480: }
! 3481: u_FpV_polint_all(x,y,C0,C1, p);
! 3482: Hp = (GEN)y[1];
! 3483: H0p= (GEN)C0[1];
! 3484: H1p= (GEN)C1[1];
! 3485: }
! 3486: else
! 3487: {
! 3488: ulong la = (ulong)leading_term(a);
! 3489: int drop;
! 3490: /* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X),
! 3491: * where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */
! 3492: for (i=0,n = 1; n <= nmax; n++)
! 3493: {
! 3494: ev = vec_u_FpX_eval_gen(b, n, p, &drop);
! 3495: i++; x[i] = n; y[i] = u_FpX_resultant(a, ev, p);
! 3496: if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p);
! 3497: ev = vec_u_FpX_eval_gen(b, p-n, p, &drop);
! 3498: i++; x[i] = p-n; y[i] = u_FpX_resultant(a, ev, p);
! 3499: if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p);
! 3500: }
! 3501: if (i == dres)
! 3502: {
! 3503: ev = vec_u_FpX_eval_gen(b, 0, p, &drop);
! 3504: i++; x[i] = 0; y[i] = u_FpX_resultant(a, ev, p);
! 3505: if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p);
! 3506: }
! 3507: Hp = u_FpV_polint(x,y, p);
! 3508: }
! 3509: if (!H && degpol(Hp) != dres) continue;
! 3510: if (checksqfree) {
! 3511: if (!u_FpX_is_squarefree(Hp, p)) goto INIT;
! 3512: if (DEBUGLEVEL>4) fprintferr("Final lambda = %ld\n",*lambda);
! 3513: checksqfree = 0;
! 3514: }
! 3515:
! 3516: if (!H)
! 3517: { /* initialize */
! 3518: q = utoi(p); stable = 0;
! 3519: H = ZX_init_CRT(Hp, p,vX);
! 3520: if (LERS) {
! 3521: H0= ZX_init_CRT(H0p, p,vX);
! 3522: H1= ZX_init_CRT(H1p, p,vX);
! 3523: }
! 3524: }
! 3525: else
! 3526: {
! 3527: GEN qp = muliu(q,p);
! 3528: stable = ZX_incremental_CRT(H, Hp, q,qp, p);
! 3529: if (LERS) {
! 3530: stable &= ZX_incremental_CRT(H0,H0p, q,qp, p);
! 3531: stable &= ZX_incremental_CRT(H1,H1p, q,qp, p);
! 3532: }
! 3533: q = qp;
! 3534: }
! 3535: /* could make it probabilistic for H ? [e.g if stable twice, etc]
! 3536: * Probabilistic anyway for H0, H1 */
! 3537: if (DEBUGLEVEL>5)
! 3538: msgtimer("resultant mod %ld (bound 2^%ld, stable=%ld)", p,expi(q),stable);
! 3539: if (stable && (ulong)expi(q) >= bound) break; /* DONE */
! 3540: if (low_stack(lim, stack_lim(av,2)))
! 3541: {
! 3542: GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1;
! 3543: if (DEBUGMEM>1) err(warnmem,"ZY_ZXY_resultant");
! 3544: gerepilemany(av2,gptr,LERS? 4: 2); b = u_allocpol(degpol(B), 0);
! 3545: }
! 3546: }
! 3547: END:
! 3548: setvarn(H, vX); if (delete) delete_var();
! 3549: if (cB) H = gmul(H, gpowgs(cB, degpol(A)));
! 3550: if (LERS)
! 3551: {
! 3552: GEN *gptr[2];
! 3553: GEN z = cgetg(3, t_VEC);
! 3554: z[1] = (long)H0;
! 3555: z[2] = (long)H1; *LERS = z;
! 3556: gptr[0]=&H; gptr[1]=LERS;
! 3557: gerepilemany(av, gptr, 2);
! 3558: return H;
! 3559: }
! 3560: if (!cB) H = gcopy(H);
! 3561: return gerepileupto(av, H);
! 3562: }
! 3563:
! 3564: GEN
! 3565: ZY_ZXY_resultant(GEN A, GEN B0, long *lambda)
! 3566: {
! 3567: return ZY_ZXY_resultant_all(A, B0, lambda, NULL);
! 3568: }
! 3569:
! 3570: /* If lambda = NULL, return caract(Mod(B, A)), A,B in Z[X].
! 3571: * Otherwise find a small lambda such that caract (Mod(B + lambda X, A)) is
! 3572: * squarefree */
! 3573: GEN
! 3574: ZX_caract_sqf(GEN A, GEN B, long *lambda, long v)
! 3575: {
! 3576: ulong av = avma;
! 3577: GEN B0, R, a;
! 3578: long dB;
! 3579: int delete;
! 3580:
! 3581: if (v < 0) v = 0;
! 3582: switch (typ(B))
! 3583: {
! 3584: case t_POL: dB = degpol(B); if (dB > 0) break;
! 3585: B = dB? (GEN)B[2]: gzero; /* fall through */
! 3586: default:
! 3587: if (lambda) { B = scalarpol(B,varn(A)); dB = 0; break;}
! 3588: return gerepileupto(av, gpowgs(gsub(polx[v], B), degpol(A)));
! 3589: }
! 3590: delete = 0;
! 3591: if (varn(A) == 0)
! 3592: {
! 3593: long v0 = fetch_var(); delete = 1;
! 3594: A = dummycopy(A); setvarn(A,v0);
! 3595: B = dummycopy(B); setvarn(B,v0);
! 3596: }
! 3597: B0 = cgetg(4, t_POL); B0[1] = evalsigne(1)|evalvarn(0)|evallgef(4);
! 3598: B0[2] = (long)gneg_i(B);
! 3599: B0[3] = un;
! 3600: R = ZY_ZXY_resultant(A, B0, lambda);
! 3601: if (delete) delete_var();
! 3602: setvarn(R, v); a = leading_term(A);
! 3603: if (!gcmp1(a)) R = gdiv(R, gpowgs(a, dB));
! 3604: return gerepileupto(av, R);
! 3605: }
! 3606:
! 3607: GEN
! 3608: ZX_caract(GEN A, GEN B, long v)
! 3609: {
! 3610: return ZX_caract_sqf(A, B, NULL, v);
! 3611: }
! 3612:
! 3613: static GEN
! 3614: trivial_case(GEN A, GEN B)
! 3615: {
! 3616: if (typ(A) == t_INT) return gpowgs(A, degpol(B));
! 3617: if (degpol(A) == 0) return trivial_case((GEN)A[2],B);
! 3618: return NULL;
! 3619: }
! 3620:
! 3621: GEN
! 3622: ZX_resultant_all(GEN A, GEN B, ulong bound)
! 3623: {
! 3624: ulong av = avma, av2, lim, Hp;
! 3625: int stable;
! 3626: GEN q, a, b, H;
! 3627: byteptr d = diffptr + 3000;
! 3628: ulong p = 27449; /* p = prime(3000) */
! 3629:
! 3630: if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
! 3631: q = H = NULL;
! 3632: av2 = avma; lim = stack_lim(av,2);
! 3633: if (!bound) bound = ZY_ZXY_ResBound(A,B);
! 3634: if (DEBUGLEVEL>4) fprintferr("bound for resultant: 2^%ld\n",bound);
! 3635:
! 3636: for(;;)
! 3637: {
! 3638: p += *d++; if (!*d) err(primer1);
! 3639: a = u_Fp_FpX(A, 0, p);
! 3640: b = u_Fp_FpX(B, 0, p);
! 3641: Hp= u_FpX_resultant(a, b, p);
! 3642: if (!H)
! 3643: {
! 3644: stable = 0; q = utoi(p);
! 3645: H = u_center(Hp, p, p>>1);
! 3646: }
! 3647: else /* could make it probabilistic ??? [e.g if stable twice, etc] */
! 3648: {
! 3649: GEN qp = muliu(q,p);
! 3650: stable = Z_incremental_CRT(&H, Hp, q,qp, p);
! 3651: q = qp;
! 3652: }
! 3653: if (DEBUGLEVEL>5)
! 3654: msgtimer("resultant mod %ld (bound 2^%ld, stable = %d)",p,expi(q),stable);
! 3655: if (stable && (ulong)expi(q) >= bound) break; /* DONE */
! 3656: if (low_stack(lim, stack_lim(av,2)))
! 3657: {
! 3658: GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q;
! 3659: if (DEBUGMEM>1) err(warnmem,"ZX_resultant");
! 3660: gerepilemany(av2,gptr, 2);
! 3661: }
! 3662: }
! 3663: return gerepileuptoint(av, icopy(H));
! 3664: }
! 3665:
! 3666: GEN
! 3667: ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,0); }
! 3668:
! 3669: /* assume x has integral coefficients */
! 3670: GEN
! 3671: ZX_disc_all(GEN x, ulong bound)
! 3672: {
! 3673: ulong av = avma;
! 3674: GEN l, d = ZX_resultant_all(x, derivpol(x), bound);
! 3675: l = leading_term(x); if (!gcmp1(l)) d = divii(d,l);
! 3676: if (degpol(x) & 2) d = negi(d);
! 3677: return gerepileuptoint(av,d);
! 3678: }
! 3679: GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
! 3680:
! 3681: int
! 3682: ZX_is_squarefree(GEN x)
! 3683: {
! 3684: ulong av = avma;
! 3685: int d = (lgef(modulargcd(x,derivpol(x))) == 3);
! 3686: avma = av; return d;
! 3687: }
! 3688:
! 3689: /* h integer, P in Z[X]. Return h^degpol(P) P(x / h) */
! 3690: GEN
! 3691: ZX_rescale_pol(GEN P, GEN h)
! 3692: {
! 3693: long i, l = lgef(P);
! 3694: GEN Q = cgetg(l,t_POL), hi = gun;
! 3695: Q[l-1] = P[l-1];
! 3696: for (i=l-2; i>=2; i--)
! 3697: {
! 3698: hi = gmul(hi,h); Q[i] = lmul((GEN)P[i], hi);
! 3699: }
! 3700: Q[1] = P[1]; return Q;
! 3701: }
! 3702:
! 3703: /* A0 and B0 in Q[X] */
! 3704: GEN
! 3705: modulargcd(GEN A0, GEN B0)
! 3706: {
! 3707: GEN a,b,Hp,D,A,B,q,qp,H,g,p1;
! 3708: long m,n;
! 3709: ulong p, av2, av = avma, avlim = stack_lim(av,1);
! 3710: byteptr d = diffptr;
! 3711:
! 3712: if ((typ(A0) | typ(B0)) !=t_POL) err(notpoler,"modulargcd");
! 3713: if (!signe(A0)) return gcopy(B0);
! 3714: if (!signe(B0)) return gcopy(A0);
! 3715: A = content(A0);
! 3716: B = content(B0); D = ggcd(A,B);
! 3717: A = gcmp1(A)? A0: gdiv(A0,A);
! 3718: B = gcmp1(B)? B0: gdiv(B0,B);
! 3719: /* A, B in Z[X] */
! 3720: g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */
! 3721: if (degpol(A) < degpol(B)) swap(A, B);
! 3722: n = 1 + degpol(B); /* > degree(gcd) */
! 3723:
! 3724: av2 = avma; H = NULL;
! 3725: d += 3000; /* 27449 = prime(3000) */
! 3726: for(p = 27449; ; p+= *d++)
! 3727: {
! 3728: if (!*d) err(primer1);
! 3729: if (!umodiu(g,p)) continue;
! 3730:
! 3731: a = u_Fp_FpX(A, 0, p);
! 3732: b = u_Fp_FpX(B, 0, p); Hp = u_FpX_gcd_i(a,b, p);
! 3733: m = degpol(Hp);
! 3734: if (m == 0) { H = polun[varn(A0)]; break; } /* coprime. DONE */
! 3735: if (m > n) continue; /* p | Res(A/G, B/G). Discard */
! 3736:
! 3737: if (is_pm1(g)) /* make sure lead(H) = g mod p */
! 3738: Hp = u_FpX_normalize(Hp, p);
! 3739: else
! 3740: {
! 3741: ulong t = mulssmod(umodiu(g, p), u_invmod(Hp[m+2],p), p);
! 3742: Hp = u_FpX_Fp_mul(Hp, t, p, 0);
! 3743: }
! 3744: if (m < n)
! 3745: { /* First time or degree drop [all previous p were as above; restart]. */
! 3746: H = ZX_init_CRT(Hp,p,varn(A0));
! 3747: q = utoi(p); n = m; continue;
! 3748: }
! 3749:
! 3750: qp = muliu(q,p);
! 3751: if (ZX_incremental_CRT(H, Hp, q, qp, p))
! 3752: { /* H stable: check divisibility */
! 3753: if (!is_pm1(g)) { p1 = content(H); if (!is_pm1(p1)) H = gdiv(H,p1); }
! 3754: if (!signe(gres(A,H)) && !signe(gres(B,H))) break; /* DONE */
! 3755:
! 3756: if (DEBUGLEVEL) fprintferr("modulargcd: trial division failed");
! 3757: }
! 3758: q = qp;
! 3759: if (low_stack(avlim, stack_lim(av,1)))
! 3760: {
! 3761: GEN *gptr[2]; gptr[0]=&H; gptr[1]=&q;
! 3762: if (DEBUGMEM>1) err(warnmem,"modulargcd");
! 3763: gerepilemany(av2,gptr,2);
! 3764: }
! 3765: }
! 3766: return gerepileupto(av, gmul(D,H));
! 3767: }
! 3768:
! 3769: /* lift(1 / Mod(A,B)) */
! 3770: GEN
! 3771: ZX_invmod(GEN A0, GEN B0)
! 3772: {
! 3773: GEN a,b,D,A,B,q,qp,Up,Vp,U,V,res;
! 3774: long stable;
! 3775: ulong p, av2, av = avma, avlim = stack_lim(av,1);
! 3776: byteptr d = diffptr;
! 3777:
! 3778: if (typ(B0) != t_POL) err(notpoler,"ZX_invmod");
! 3779: if (typ(A0) != t_POL)
! 3780: {
! 3781: if (is_scalar_t(typ(A0))) return ginv(A0);
! 3782: err(notpoler,"ZX_invmod");
! 3783: }
! 3784: A = content(A0); D = A;
! 3785: B = content(B0);
! 3786: A = gcmp1(A)? A0: gdiv(A0,A);
! 3787: B = gcmp1(B)? B0: gdiv(B0,B);
! 3788: /* A, B in Z[X] */
! 3789: av2 = avma; U = NULL;
! 3790: d += 3000; /* 27449 = prime(3000) */
! 3791: for(p = 27449; ; p+= *d++)
! 3792: {
! 3793: if (!*d) err(primer1);
! 3794: a = u_Fp_FpX(A, 0, p);
! 3795: b = u_Fp_FpX(B, 0, p);
! 3796: /* if p | Res(A/G, B/G), discard */
! 3797: if (!u_FpX_extresultant(b,a,p, &Vp,&Up)) continue;
! 3798:
! 3799: if (!U)
! 3800: { /* First time */
! 3801: U = ZX_init_CRT(Up,p,varn(A0));
! 3802: V = ZX_init_CRT(Vp,p,varn(A0));
! 3803: q = utoi(p); continue;
! 3804: }
! 3805: if (DEBUGLEVEL>5) msgtimer("ZX_invmod: mod %ld (bound 2^%ld)", p,expi(q));
! 3806: qp = muliu(q,p);
! 3807: stable = ZX_incremental_CRT(U, Up, q,qp, p);
! 3808: stable &= ZX_incremental_CRT(V, Vp, q,qp, p);
! 3809: if (stable)
! 3810: { /* all stable: check divisibility */
! 3811: res = gadd(gmul(A,U), gmul(B,V));
! 3812: if (degpol(res) == 0) break; /* DONE */
! 3813: if (DEBUGLEVEL) fprintferr("ZX_invmod: char 0 check failed");
! 3814: }
! 3815: q = qp;
! 3816: if (low_stack(avlim, stack_lim(av,1)))
! 3817: {
! 3818: GEN *gptr[3]; gptr[0]=&q; gptr[1]=&U; gptr[2]=&V;
! 3819: if (DEBUGMEM>1) err(warnmem,"ZX_invmod");
! 3820: gerepilemany(av2,gptr,3);
! 3821: }
! 3822: }
! 3823: D = gmul(D,res);
! 3824: return gerepileupto(av, gdiv(U,D));
! 3825: }
! 3826:
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-2.2 / src / basemath