Annotation of OpenXM_contrib/pari-2.2/src/basemath/polarit3.c, Revision 1.1.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