Annotation of OpenXM_contrib/pari-2.2/src/basemath/polarit3.c, Revision 1.2
1.2 ! noro 1: /* $Id: polarit3.c,v 1.147 2002/09/07 15:46:28 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /***********************************************************************/
17: /** **/
18: /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
19: /** (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:
1.2 ! noro 36: #define both_odd(x,y) ((x)&(y)&1)
! 37: extern GEN caractducos(GEN p, GEN x, int v);
! 38: extern double mylog2(GEN z);
! 39: extern GEN revpol(GEN x);
! 40:
1.1 noro 41: /*******************************************************************/
42: /* */
43: /* KARATSUBA (for polynomials) */
44: /* */
45: /*******************************************************************/
46:
47: #if 1 /* for tunings */
48: long SQR_LIMIT = 6;
49: long MUL_LIMIT = 10;
50: long u_SQR_LIMIT = 6;
51: long u_MUL_LIMIT = 10;
52:
53: void
54: setsqpol(long a) { SQR_LIMIT=a; u_SQR_LIMIT=a; }
55: void
56: setmulpol(long a) { MUL_LIMIT=a; u_MUL_LIMIT=a; }
57:
58: GEN
59: u_specpol(GEN x, long nx)
60: {
61: GEN z = cgetg(nx+2,t_POL);
62: long i;
63: for (i=0; i<nx; i++) z[i+2] = lstoi(x[i]);
64: z[1]=evalsigne(1)|evallgef(nx+2);
65: return z;
66: }
67:
68: GEN
69: specpol(GEN x, long nx)
70: {
71: GEN z = cgetg(nx+2,t_POL);
72: long i;
73: for (i=0; i<nx; i++) z[i+2] = x[i];
74: z[1]=evalsigne(1)|evallgef(nx+2);
75: return z;
76: }
77: #else
78: # define SQR_LIMIT 6
79: # define MUL_LIMIT 10
80: #endif
81:
82: #ifndef HIGHWORD
83: # define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG)
84: #endif
85:
86: /* multiplication in Fp[X], p small */
87:
1.2 ! noro 88: GEN
1.1 noro 89: u_normalizepol(GEN x, long lx)
90: {
91: long i;
92: for (i=lx-1; i>1; i--)
93: if (x[i]) break;
94: setlgef(x,i+1);
95: setsigne(x, (i > 1)); return x;
96: }
97:
98: GEN
99: u_FpX_deriv(GEN z, ulong p)
100: {
101: long i,l = lgef(z)-1;
102: GEN x;
103: if (l < 2) l = 2;
104: x = cgetg(l,t_VECSMALL); x[1] = z[1]; z++;
105: if (HIGHWORD(l | p))
106: for (i=2; i<l; i++) x[i] = mulssmod(i-1, z[i], p);
107: else
108: for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
109: return u_normalizepol(x,l);
110: }
111:
112: static GEN
113: u_FpX_gcd_i(GEN a, GEN b, ulong p)
114: {
115: GEN c;
116: if (lgef(b) > lgef(a)) swap(a, b);
117: while (signe(b))
118: {
119: c = u_FpX_rem(a,b,p);
120: a = b; b = c;
121: }
122: return a;
123: }
124:
125: GEN
126: u_FpX_gcd(GEN a, GEN b, ulong p)
127: {
1.2 ! noro 128: gpmem_t av = avma;
1.1 noro 129: return gerepileupto(av, u_FpX_gcd_i(a,b,p));
130: }
131:
132: int
133: u_FpX_is_squarefree(GEN z, ulong p)
134: {
1.2 ! noro 135: gpmem_t av = avma;
1.1 noro 136: GEN d = u_FpX_gcd_i(z, u_FpX_deriv(z,p) , p);
137: avma = av; return (degpol(d) == 0);
138: }
139:
140: static GEN
141: u_FpX_addspec(GEN x, GEN y, long p, long lx, long ly)
142: {
143: long i,lz;
144: GEN z;
145:
146: if (ly>lx) swapspec(x,y, lx,ly);
147: lz = lx+2; z = cgetg(lz,t_VECSMALL) + 2;
148: for (i=0; i<ly; i++) z[i] = addssmod(x[i],y[i],p);
149: for ( ; i<lx; i++) z[i] = x[i];
150: z -= 2; z[1]=0; return u_normalizepol(z, lz);
151: }
152:
153: #ifdef INLINE
154: INLINE
155: #endif
156: ulong
157: u_FpX_mullimb(GEN x, GEN y, ulong p, long a, long b)
158: {
159: ulong p1 = 0;
160: long i;
161: for (i=a; i<b; i++)
162: if (y[i])
163: {
164: p1 += y[i] * x[-i];
165: if (p1 & HIGHBIT) p1 %= p;
166: }
167: return p1 % p;
168: }
169:
170: ulong
171: u_FpX_mullimb_safe(GEN x, GEN y, ulong p, long a, long b)
172: {
173: ulong p1 = 0;
174: long i;
175: for (i=a; i<b; i++)
176: if (y[i])
177: p1 = addssmod(p1, mulssmod(y[i],x[-i],p), p);
178: return p1;
179: }
180:
181: /* assume nx >= ny > 0 */
182: static GEN
183: s_FpX_mulspec(GEN x, GEN y, ulong p, long nx, long ny)
184: {
185: long i,lz,nz;
186: GEN z;
187:
188: lz = nx+ny+1; nz = lz-2;
189: z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
190: if (u_OK_ULONG(p))
191: {
192: for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb(x+i,y,p,0,i+1);
193: for ( ; i<nx; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,0,ny);
194: for ( ; i<nz; i++) z[i] = (long)u_FpX_mullimb(x+i,y,p,i-nx+1,ny);
195: }
196: else
197: {
198: for (i=0; i<ny; i++)z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,i+1);
199: for ( ; i<nx; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,0,ny);
200: for ( ; i<nz; i++) z[i] = (long)u_FpX_mullimb_safe(x+i,y,p,i-nx+1,ny);
201: }
202: z -= 2; z[1]=0; return u_normalizepol(z, lz);
203: }
204:
205: /* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
206: GEN
207: u_FpX_addshift(GEN x, GEN y, ulong p, long d)
208: {
209: GEN xd,yd,zd = (GEN)avma;
210: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
211:
212: x += 2; y += 2; a = ny-d;
213: if (a <= 0)
214: {
215: lz = (a>nx)? ny+2: nx+d+2;
216: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
217: while (xd > x) *--zd = *--xd;
218: x = zd + a;
219: while (zd > x) *--zd = 0;
220: }
221: else
222: {
223: xd = new_chunk(d); yd = y+d;
224: x = u_FpX_addspec(x,yd,p, nx,a);
225: lz = (a>nx)? ny+2: lgef(x)+d;
226: x += 2; while (xd > x) *--zd = *--xd;
227: }
228: while (yd > y) *--zd = *--yd;
229: *--zd = evalsigne(1) | evallgef(lz);
230: *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
231: }
232:
233: #define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2)
234: #define u_FpX_sqr(x, p) u_FpX_Ksqr(x+2,p, lgef(x)-2)
235: #define u_FpX_add(x,y, p) u_FpX_addspec(x+2,y+2,p, lgef(x)-2,lgef(y)-2)
236:
1.2 ! noro 237: GEN
! 238: u_zeropol()
! 239: {
! 240: GEN z = cgetg(2, t_VECSMALL);
! 241: z[1] = evallgef(2) | evalvarn(0); return z;
! 242: }
! 243:
1.1 noro 244: static GEN
1.2 ! noro 245: u_mallocpol(long d)
1.1 noro 246: {
1.2 ! noro 247: GEN z = (GEN)gpmalloc((d+3) * sizeof(long));
1.1 noro 248: z[0] = evaltyp(t_VECSMALL) | evallg(d+3);
249: z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0);
250: return z;
251: }
1.2 ! noro 252: static GEN
! 253: u_getpol(long d)
! 254: {
! 255: GEN z = cgetg(d + 3, t_VECSMALL);
! 256: z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0);
! 257: return z;
! 258: }
1.1 noro 259:
260: static GEN
1.2 ! noro 261: u_scalarpol(ulong x)
1.1 noro 262: {
1.2 ! noro 263: GEN z;
! 264: if (!x) return u_zeropol();
! 265: z = u_getpol(0);
1.1 noro 266: z[2] = x; return z;
267: }
268:
269: static GEN
1.2 ! noro 270: u_FpX_neg(GEN x, ulong p)
! 271: {
! 272: long i, l = lgef(x);
! 273: GEN z = cgetg(l, t_VECSMALL);
! 274: z[1] = x[1];
! 275: for (i=2; i<l; i++) z[i] = x[i]? p - x[i]: 0;
! 276: return z;
! 277: }
! 278:
! 279: static GEN
1.1 noro 280: u_FpX_neg_i(GEN x, ulong p)
281: {
282: long i, l = lgef(x);
283: for (i=2; i<l; i++)
284: if (x[i]) x[i] = p - x[i];
285: return x;
286: }
287:
288: /* shift polynomial + gerepile */
289: static GEN
1.2 ! noro 290: u_FpX_shiftip(gpmem_t av, GEN x, long v)
1.1 noro 291: {
292: long i, lx = lgef(x), ly;
293: GEN y;
294: if (v <= 0 || !signe(x)) return gerepileupto(av, x);
295: avma = av; ly = lx + v;
296: x += lx; y = new_chunk(ly) + ly;
297: for (i = 2; i<lx; i++) *--y = *--x;
298: for (i = 0; i< v; i++) *--y = 0;
299: *--y = evalsigne(1) | evallgef(ly);
300: *--y = evaltyp(t_VECSMALL) | evallg(ly); return y;
301: }
302:
303: /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
304: * b+2 were sent instead. na, nb = number of terms of a, b.
305: * Only c, c0, c1, c2 are genuine GEN.
306: */
307: GEN
308: u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb)
309: {
310: GEN a0,c,c0;
1.2 ! noro 311: long n0, n0a, i, v = 0;
! 312: gpmem_t av;
1.1 noro 313:
314: while (na && !a[0]) { a++; na--; v++; }
315: while (nb && !b[0]) { b++; nb--; v++; }
316: if (na < nb) swapspec(a,b, na,nb);
1.2 ! noro 317: if (!nb) return u_zeropol();
1.1 noro 318:
319: av = avma;
320: if (nb < u_MUL_LIMIT)
321: return u_FpX_shiftip(av,s_FpX_mulspec(a,b,p,na,nb), v);
322: i=(na>>1); n0=na-i; na=i;
323: a0=a+n0; n0a=n0;
324: while (n0a && !a[n0a-1]) n0a--;
325:
326: if (nb > n0)
327: {
328: GEN b0,c1,c2;
329: long n0b;
330:
331: nb -= n0; b0 = b+n0; n0b = n0;
332: while (n0b && !b[n0b-1]) n0b--;
333: c = u_FpX_Kmul(a,b,p,n0a,n0b);
334: c0 = u_FpX_Kmul(a0,b0,p,na,nb);
335:
336: c2 = u_FpX_addspec(a0,a,p,na,n0a);
337: c1 = u_FpX_addspec(b0,b,p,nb,n0b);
338:
339: c1 = u_FpX_mul(c1,c2,p);
340: c2 = u_FpX_add(c0,c,p);
341:
342: c2 = u_FpX_neg_i(c2,p);
343: c2 = u_FpX_add(c1,c2,p);
344: c0 = u_FpX_addshift(c0,c2 ,p, n0);
345: }
346: else
347: {
348: c = u_FpX_Kmul(a,b,p,n0a,nb);
349: c0 = u_FpX_Kmul(a0,b,p,na,nb);
350: }
351: c0 = u_FpX_addshift(c0,c,p,n0);
352: return u_FpX_shiftip(av,c0, v);
353: }
354:
355: GEN
356: u_FpX_sqrpol(GEN x, ulong p, long nx)
357: {
358: long i,j,l,lz,nz;
359: ulong p1;
360: GEN z;
361:
1.2 ! noro 362: if (!nx) return u_zeropol();
1.1 noro 363: lz = (nx << 1) + 1, nz = lz-2;
364: z = cgetg(lz,t_VECSMALL) + 2;
365: for (i=0; i<nx; i++)
366: {
367: p1 = 0; l = (i+1)>>1;
368: for (j=0; j<l; j++)
369: {
370: p1 += x[j] * x[i-j];
371: if (p1 & HIGHBIT) p1 %= p;
372: }
373: p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
374: if ((i&1) == 0)
375: p1 += x[i>>1] * x[i>>1];
376: z[i] = p1 % p;
377: }
378: for ( ; i<nz; i++)
379: {
380: p1 = 0; l = (i+1)>>1;
381: for (j=i-nx+1; j<l; j++)
382: {
383: p1 += x[j] * x[i-j];
384: if (p1 & HIGHBIT) p1 %= p;
385: }
386: p1 <<= 1; if (p1 & HIGHBIT) p1 %= p;
387: if ((i&1) == 0)
388: p1 += x[i>>1] * x[i>>1];
389: z[i] = p1 % p;
390: }
391: z -= 2; z[1]=0; return u_normalizepol(z,lz);
392: }
393:
394: static GEN
395: u_Fp_gmul2_1(GEN x, ulong p)
396: {
397: long i, l = lgef(x);
398: GEN z = cgetg(l, t_VECSMALL);
399: for (i=2; i<l; i++)
400: {
401: ulong p1 = x[i]<<1;
402: if (p1 > p) p1 -= p;
403: z[i] = p1;
404: }
405: z[1] = x[1]; return z;
406: }
407:
408: GEN
409: u_FpX_Ksqr(GEN a, ulong p, long na)
410: {
411: GEN a0,c,c0,c1;
1.2 ! noro 412: long n0, n0a, i, v = 0;
! 413: gpmem_t av;
1.1 noro 414:
415: while (na && !a[0]) { a++; na--; v += 2; }
416: av = avma;
417: if (na < u_SQR_LIMIT) return u_FpX_shiftip(av, u_FpX_sqrpol(a,p,na), v);
418: i=(na>>1); n0=na-i; na=i;
419: a0=a+n0; n0a=n0;
420: while (n0a && !a[n0a-1]) n0a--;
421:
422: c = u_FpX_Ksqr(a,p,n0a);
423: c0= u_FpX_Ksqr(a0,p,na);
424: if (p == 2) n0 *= 2;
425: else
426: {
427: c1 = u_Fp_gmul2_1(u_FpX_Kmul(a0,a,p, na,n0a), p);
428: c0 = u_FpX_addshift(c0,c1,p,n0);
429: }
430: c0 = u_FpX_addshift(c0,c,p,n0);
431: return u_FpX_shiftip(av,c0,v);
432: }
433:
434: /* generic multiplication */
435:
436: static GEN
437: addpol(GEN x, GEN y, long lx, long ly)
438: {
439: long i,lz;
440: GEN z;
441:
442: if (ly>lx) swapspec(x,y, lx,ly);
443: lz = lx+2; z = cgetg(lz,t_POL) + 2;
444: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
445: for ( ; i<lx; i++) z[i]=x[i];
446: z -= 2; z[1]=0; return normalizepol_i(z, lz);
447: }
448:
449: static GEN
450: addpolcopy(GEN x, GEN y, long lx, long ly)
451: {
452: long i,lz;
453: GEN z;
454:
455: if (ly>lx) swapspec(x,y, lx,ly);
456: lz = lx+2; z = cgetg(lz,t_POL) + 2;
457: for (i=0; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
458: for ( ; i<lx; i++) z[i]=lcopy((GEN)x[i]);
459: z -= 2; z[1]=0; return normalizepol_i(z, lz);
460: }
461:
462: #ifdef INLINE
463: INLINE
464: #endif
465: GEN
466: mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b)
467: {
468: GEN p1 = NULL;
1.2 ! noro 469: long i;
! 470: gpmem_t av = avma;
1.1 noro 471: for (i=a; i<b; i++)
472: if (ynonzero[i])
473: {
474: GEN p2 = gmul((GEN)y[i],(GEN)x[-i]);
475: p1 = p1 ? gadd(p1, p2): p2;
476: }
477: return p1 ? gerepileupto(av, p1): gzero;
478: }
479:
480: /* assume nx >= ny > 0 */
481: static GEN
482: mulpol(GEN x, GEN y, long nx, long ny)
483: {
484: long i,lz,nz;
485: GEN z;
486: char *p1;
487:
488: lz = nx+ny+1; nz = lz-2;
489: z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
490: p1 = gpmalloc(ny);
491: for (i=0; i<ny; i++)
492: {
493: p1[i] = !isexactzero((GEN)y[i]);
494: z[i] = (long)mulpol_limb(x+i,y,p1,0,i+1);
495: }
496: for ( ; i<nx; i++) z[i] = (long)mulpol_limb(x+i,y,p1,0,ny);
497: for ( ; i<nz; i++) z[i] = (long)mulpol_limb(x+i,y,p1,i-nx+1,ny);
498: free(p1); z -= 2; z[1]=0; return normalizepol_i(z, lz);
499: }
500:
501: /* return (x * X^d) + y. Assume d > 0, x > 0 and y >= 0 */
502: GEN
503: addshiftw(GEN x, GEN y, long d)
504: {
505: GEN xd,yd,zd = (GEN)avma;
506: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
507:
508: x += 2; y += 2; a = ny-d;
509: if (a <= 0)
510: {
511: lz = (a>nx)? ny+2: nx+d+2;
512: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
513: while (xd > x) *--zd = *--xd;
514: x = zd + a;
515: while (zd > x) *--zd = zero;
516: }
517: else
518: {
519: xd = new_chunk(d); yd = y+d;
520: x = addpol(x,yd, nx,a);
521: lz = (a>nx)? ny+2: lgef(x)+d;
522: x += 2; while (xd > x) *--zd = *--xd;
523: }
524: while (yd > y) *--zd = *--yd;
525: *--zd = evalsigne(1) | evallgef(lz);
526: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
527: }
528:
529: GEN
530: addshiftpol(GEN x, GEN y, long d)
531: {
532: long v = varn(x);
533: if (!signe(x)) return y;
534: x = addshiftw(x,y,d);
535: setvarn(x,v); return x;
536: }
537:
538: /* as above, producing a clean malloc */
539: static GEN
540: addshiftwcopy(GEN x, GEN y, long d)
541: {
542: GEN xd,yd,zd = (GEN)avma;
543: long a,lz,ny = lgef(y)-2, nx = lgef(x)-2;
544:
545: x += 2; y += 2; a = ny-d;
546: if (a <= 0)
547: {
548: lz = nx+d+2;
549: (void)new_chunk(lz); xd = x+nx; yd = y+ny;
550: while (xd > x) *--zd = lcopy((GEN)*--xd);
551: x = zd + a;
552: while (zd > x) *--zd = zero;
553: }
554: else
555: {
556: xd = new_chunk(d); yd = y+d;
557: x = addpolcopy(x,yd, nx,a);
558: lz = (a>nx)? ny+2: lgef(x)+d;
559: x += 2; while (xd > x) *--zd = *--xd;
560: }
561: while (yd > y) *--zd = lcopy((GEN)*--yd);
562: *--zd = evalsigne(1) | evallgef(lz);
563: *--zd = evaltyp(t_POL) | evallg(lz); return zd;
564: }
565:
566: /* shift polynomial in place. assume v free cells have been left before x */
567: static GEN
568: shiftpol_ip(GEN x, long v)
569: {
570: long i, lx;
571: GEN y;
572: if (v <= 0 || !signe(x)) return x;
573: lx = lgef(x);
574: x += 2; y = x + v;
575: for (i = lx-3; i>=0; i--) y[i] = x[i];
576: for (i = 0 ; i< v; i++) x[i] = zero;
577: lx += v;
578: *--x = evalsigne(1) | evallgef(lx);
579: *--x = evaltyp(t_POL) | evallg(lx); return x;
580: }
581:
582: /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
583: * b+2 were sent instead. na, nb = number of terms of a, b.
584: * Only c, c0, c1, c2 are genuine GEN.
585: */
586: GEN
587: quickmul(GEN a, GEN b, long na, long nb)
588: {
589: GEN a0,c,c0;
1.2 ! noro 590: long n0, n0a, i, v = 0;
! 591: gpmem_t av;
1.1 noro 592:
593: while (na && isexactzero((GEN)a[0])) { a++; na--; v++; }
594: while (nb && isexactzero((GEN)b[0])) { b++; nb--; v++; }
595: if (na < nb) swapspec(a,b, na,nb);
596: if (!nb) return zeropol(0);
597:
598: if (v) (void)cgetg(v,t_STR); /* v gerepile-safe cells for shiftpol_ip */
599: if (nb < MUL_LIMIT)
600: return shiftpol_ip(mulpol(a,b,na,nb), v);
601: i=(na>>1); n0=na-i; na=i;
602: av=avma; a0=a+n0; n0a=n0;
603: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
604:
605: if (nb > n0)
606: {
607: GEN b0,c1,c2;
608: long n0b;
609:
610: nb -= n0; b0 = b+n0; n0b = n0;
611: while (n0b && isexactzero((GEN)b[n0b-1])) n0b--;
612: c = quickmul(a,b,n0a,n0b);
613: c0 = quickmul(a0,b0, na,nb);
614:
615: c2 = addpol(a0,a, na,n0a);
616: c1 = addpol(b0,b, nb,n0b);
617:
618: c1 = quickmul(c1+2,c2+2, lgef(c1)-2,lgef(c2)-2);
619: c2 = gneg_i(gadd(c0,c));
620: c0 = addshiftw(c0, gadd(c1,c2), n0);
621: }
622: else
623: {
624: c = quickmul(a,b,n0a,nb);
625: c0 = quickmul(a0,b,na,nb);
626: }
627: c0 = addshiftwcopy(c0,c,n0);
628: return shiftpol_ip(gerepileupto(av,c0), v);
629: }
630:
631: GEN
632: sqrpol(GEN x, long nx)
633: {
1.2 ! noro 634: long i, j, l, lz, nz;
! 635: gpmem_t av;
1.1 noro 636: GEN p1,z;
637: char *p2;
638:
639: if (!nx) return zeropol(0);
640: lz = (nx << 1) + 1, nz = lz-2;
641: z = cgetg(lz,t_POL) + 2;
642: p2 = gpmalloc(nx);
643: for (i=0; i<nx; i++)
644: {
645: p2[i] = !isexactzero((GEN)x[i]);
646: p1=gzero; av=avma; l=(i+1)>>1;
647: for (j=0; j<l; j++)
648: if (p2[j] && p2[i-j])
649: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
650: p1 = gshift(p1,1);
651: if ((i&1) == 0 && p2[i>>1])
652: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
653: z[i] = lpileupto(av,p1);
654: }
655: for ( ; i<nz; i++)
656: {
657: p1=gzero; av=avma; l=(i+1)>>1;
658: for (j=i-nx+1; j<l; j++)
659: if (p2[j] && p2[i-j])
660: p1 = gadd(p1, gmul((GEN)x[j],(GEN)x[i-j]));
661: p1 = gshift(p1,1);
662: if ((i&1) == 0 && p2[i>>1])
663: p1 = gadd(p1, gsqr((GEN)x[i>>1]));
664: z[i] = lpileupto(av,p1);
665: }
666: free(p2); z -= 2; z[1]=0; return normalizepol_i(z,lz);
667: }
668:
669: GEN
670: quicksqr(GEN a, long na)
671: {
672: GEN a0,c,c0,c1;
1.2 ! noro 673: long n0, n0a, i, v = 0;
! 674: gpmem_t av;
1.1 noro 675:
676: while (na && isexactzero((GEN)a[0])) { a++; na--; v += 2; }
677: if (v) (void)new_chunk(v);
678: if (na<SQR_LIMIT) return shiftpol_ip(sqrpol(a,na), v);
679: i=(na>>1); n0=na-i; na=i;
680: av=avma; a0=a+n0; n0a=n0;
681: while (n0a && isexactzero((GEN)a[n0a-1])) n0a--;
682:
683: c = quicksqr(a,n0a);
684: c0 = quicksqr(a0,na);
685: c1 = gmul2n(quickmul(a0,a, na,n0a), 1);
686: c0 = addshiftw(c0,c1, n0);
687: c0 = addshiftwcopy(c0,c,n0);
688: return shiftpol_ip(gerepileupto(av,c0), v);
689: }
690: /*****************************************
691: * Arithmetic in Z/pZ[X] *
692: *****************************************/
693:
694: /*********************************************************************
695: This functions supposes polynomials already reduced.
696: There are clean and memory efficient.
697: **********************************************************************/
698:
699: GEN
700: FpX_center(GEN T,GEN mod)
701: {/*OK centermod exists, but is not so clean*/
1.2 ! noro 702: gpmem_t av;
1.1 noro 703: long i, l=lg(T);
704: GEN P,mod2;
705: P=cgetg(l,t_POL);
706: P[1]=T[1];
707: av=avma;
708: mod2=gclone(shifti(mod,-1));/*clone*/
709: avma=av;
710: for(i=2;i<l;i++)
711: P[i]=cmpii((GEN)T[i],mod2)<0?licopy((GEN)T[i]):lsubii((GEN)T[i],mod);
712: gunclone(mod2);/*unclone*/
713: return P;
714: }
715:
716: GEN
717: FpX_neg(GEN x,GEN p)
718: {
719: long i,d=lgef(x);
720: GEN y;
721: y=cgetg(d,t_POL); y[1]=x[1];
722: for(i=2;i<d;i++)
723: if (signe(x[i])) y[i]=lsubii(p,(GEN)x[i]);
724: else y[i]=zero;
725: return y;
726: }
727: /**********************************************************************
728: Unclean functions, do not garbage collect.
729: This is a feature: The malloc is corrupted only by the call to FpX_red
730: so garbage collecting so often is not desirable.
731: FpX_red can sometime be avoided by passing NULL for p.
732: In this case the function is usually clean (see below for detail)
733: Added to help not using POLMOD of INTMOD which are deadly slow.
734: gerepileupto of the result is legible. Bill.
735: I don't like C++. I am wrong.
736: **********************************************************************/
737: /*
738: *If p is NULL no reduction is performed and the function is clean.
739: * for FpX_add,FpX_mul,FpX_sqr,FpX_Fp_mul
740: */
741: GEN
742: FpX_add(GEN x,GEN y,GEN p)
743: {
744: long lx,ly,i;
745: GEN z;
746: lx = lgef(x); ly = lgef(y); if (lx < ly) swapspec(x,y, lx,ly);
747: z = cgetg(lx,t_POL); z[1] = x[1];
748: for (i=2; i<ly; i++) z[i]=laddii((GEN)x[i],(GEN)y[i]);
749: for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]);
750: (void)normalizepol_i(z, lx);
1.2 ! noro 751: if (lgef(z) == 2) { avma = (gpmem_t)(z + lx); z = zeropol(varn(x)); }
1.1 noro 752: if (p) z= FpX_red(z, p);
753: return z;
754: }
755: GEN
756: FpX_sub(GEN x,GEN y,GEN p)
757: {
758: long lx,ly,i,lz;
759: GEN z;
760: lx = lgef(x); ly = lgef(y);
761: lz=max(lx,ly);
762: z = cgetg(lz,t_POL);
763: if (lx >= ly)
764: {
765: z[1] = x[1];
766: for (i=2; i<ly; i++) z[i]=lsubii((GEN)x[i],(GEN)y[i]);
767: for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]);
768: (void)normalizepol_i(z, lz);
769: }
770: else
771: {
772: z[1] = y[1];
773: for (i=2; i<lx; i++) z[i]=lsubii((GEN)x[i],(GEN)y[i]);
774: for ( ; i<ly; i++) z[i]=lnegi((GEN)y[i]);
775: /*polynomial is always normalized*/
776: }
1.2 ! noro 777: if (lgef(z) == 2) { avma = (gpmem_t)(z + lz); z = zeropol(varn(x)); }
1.1 noro 778: if (p) z= FpX_red(z, p);
779: return z;
780: }
781: GEN
782: FpX_mul(GEN x,GEN y,GEN p)
783: {
784: GEN z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2);
785: setvarn(z,varn(x));
786: if (!p) return z;
787: return FpX_red(z, p);
788: }
789: GEN
790: FpX_sqr(GEN x,GEN p)
791: {
792: GEN z = quicksqr(x+2, lgef(x)-2);
793: setvarn(z,varn(x));
794: if (!p) return z;
795: return FpX_red(z, p);
796: }
1.2 ! noro 797: #define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL)
1.1 noro 798:
799: /* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
800: static GEN
801: u_FpXQ_mul(GEN y,GEN x,GEN pol,ulong p)
802: {
803: GEN z = u_FpX_mul(y,x,p);
804: return u_FpX_rem(z,pol,p);
805: }
806: /* Square of y in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */
807: static GEN
808: u_FpXQ_sqr(GEN y,GEN pol,ulong p)
809: {
810: GEN z = u_FpX_sqr(y,p);
811: return u_FpX_rem(z,pol,p);
812: }
813:
814: /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
815: * return lift(1 / (x mod (p,pol))) */
816: GEN
817: FpXQ_invsafe(GEN x, GEN T, GEN p)
818: {
819: GEN z, U, V;
820:
821: if (typ(x) != t_POL) return mpinvmod(x, p);
822: z = FpX_extgcd(x, T, p, &U, &V);
1.2 ! noro 823: if (degpol(z)) return NULL;
1.1 noro 824: z = mpinvmod((GEN)z[2], p);
825: return FpX_Fp_mul(U, z, p);
826: }
827:
828: /* Product of y and x in Z/pZ[X]/(T)
829: * return lift(lift(Mod(x*y*Mod(1,p),T*Mod(1,p)))) */
1.2 ! noro 830: /* x and y must be poynomials in the same var as T.
! 831: * t_INT are not allowed. Use Fq_mul instead.
! 832: */
1.1 noro 833: GEN
834: FpXQ_mul(GEN y,GEN x,GEN T,GEN p)
835: {
836: GEN z;
1.2 ! noro 837: if (typ(x) == t_INT || typ(y) == t_INT)
! 838: err(bugparier,"FpXQ_mul, t_INT are absolutely forbidden");
1.1 noro 839: z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2); setvarn(z,varn(y));
840: z = FpX_red(z, p); return FpX_res(z,T, p);
841: }
842:
843: /* Square of y in Z/pZ[X]/(pol)
844: * return lift(lift(Mod(y^2*Mod(1,p),pol*Mod(1,p)))) */
845: GEN
846: FpXQ_sqr(GEN y,GEN pol,GEN p)
847: {
848: GEN z = quicksqr(y+2,lgef(y)-2); setvarn(z,varn(y));
849: z = FpX_red(z, p); return FpX_res(z,pol, p);
850: }
851: /*Modify y[2].
852: *No reduction if p is NULL
853: */
854: GEN
855: FpX_Fp_add(GEN y,GEN x,GEN p)
856: {
857: if (!signe(x)) return y;
858: if (!signe(y))
859: return scalarpol(x,varn(y));
860: y[2]=laddii((GEN)y[2],x);
861: if (p) y[2]=lmodii((GEN)y[2],p);
862: if (!signe(y[2]) && degpol(y) == 0) return zeropol(varn(y));
863: return y;
864: }
865: GEN
866: ZX_s_add(GEN y,long x)
867: {
868: if (!x) return y;
869: if (!signe(y))
870: return scalarpol(stoi(x),varn(y));
871: y[2] = laddis((GEN)y[2],x);
872: if (!signe(y[2]) && degpol(y) == 0) return zeropol(varn(y));
873: return y;
874: }
875: /* y is a polynomial in ZZ[X] and x an integer.
1.2 ! noro 876: * If p is NULL, no reduction is performed and return x*y
1.1 noro 877: *
878: * else the result is lift(y*x*Mod(1,p))
879: */
880: GEN
881: FpX_Fp_mul(GEN y,GEN x,GEN p)
882: {
883: GEN z;
884: int i;
885: if (!signe(x))
886: return zeropol(varn(y));
887: z=cgetg(lg(y),t_POL);
888: z[1]=y[1];
889: for(i=2;i<lgef(y);i++)
890: z[i]=lmulii((GEN)y[i],x);
891: if(!p) return z;
892: return FpX_red(z,p);
893: }
894: /*****************************************************************
895: * End of unclean functions. *
896: * *
897: *****************************************************************/
1.2 ! noro 898:
! 899: /* modular power */
! 900: ulong
! 901: powuumod(ulong x, ulong n0, ulong p)
! 902: {
! 903: ulong y, z, n;
! 904: if (n0 <= 2)
! 905: { /* frequent special cases */
! 906: if (n0 == 2) return mulssmod(x,x,p);
! 907: if (n0 == 1) return x;
! 908: if (n0 == 0) return 1;
! 909: }
! 910: y = 1; z = x; n = n0;
! 911: for(;;)
! 912: {
! 913: if (n&1) y = mulssmod(y,z,p);
! 914: n>>=1; if (!n) return y;
! 915: z = mulssmod(z,z,p);
! 916: }
! 917: }
! 918:
! 919: GEN
! 920: powgumod(GEN x, ulong n0, GEN p)
! 921: {
! 922: GEN y, z;
! 923: ulong n;
! 924: if (n0 <= 2)
! 925: { /* frequent special cases */
! 926: if (n0 == 2) return resii(sqri(x),p);
! 927: if (n0 == 1) return x;
! 928: if (n0 == 0) return gun;
! 929: }
! 930: y = gun; z = x; n = n0;
! 931: for(;;)
! 932: {
! 933: if (n&1) y = resii(mulii(y,z),p);
! 934: n>>=1; if (!n) return y;
! 935: z = resii(sqri(z),p);
! 936: }
! 937: }
! 938:
1.1 noro 939: /*****************************************************************
940: Clean and with no reduced hypothesis. Beware that some operations
941: will be much slower with big unreduced coefficient
942: *****************************************************************/
943: /* Inverse of x in Z[X] / (p,T)
944: * return lift(lift(Mod(x*Mod(1,p), T*Mod(1,p))^-1)); */
945: GEN
946: FpXQ_inv(GEN x,GEN T,GEN p)
947: {
1.2 ! noro 948: gpmem_t av;
1.1 noro 949: GEN U;
950:
951: if (!T) return mpinvmod(x,p);
952: av = avma;
953: U = FpXQ_invsafe(x, T, p);
954: if (!U) err(talker,"non invertible polynomial in FpXQ_inv");
955: return gerepileupto(av, U);
956: }
957:
958: GEN
1.2 ! noro 959: FpXQ_div(GEN x,GEN y,GEN T,GEN p)
! 960: {
! 961: gpmem_t av = avma;
! 962: return gerepileupto(av, FpXQ_mul(x,FpXQ_inv(y,T,p),T,p));
! 963: }
! 964:
! 965: GEN
! 966: FpXV_FpV_innerprod(GEN V, GEN W, GEN p)
1.1 noro 967: {
1.2 ! noro 968: gpmem_t ltop=avma;
1.1 noro 969: long i;
970: GEN z = FpX_Fp_mul((GEN)V[1],(GEN)W[1],NULL);
971: for(i=2;i<lg(V);i++)
972: z=FpX_add(z,FpX_Fp_mul((GEN)V[i],(GEN)W[i],NULL),NULL);
973: return gerepileupto(ltop,FpX_red(z,p));
974: }
975:
976: /* generates the list of powers of x of degree 0,1,2,...,l*/
977: GEN
978: FpXQ_powers(GEN x, long l, GEN T, GEN p)
979: {
980: GEN V=cgetg(l+2,t_VEC);
981: long i;
982: V[1] = (long) scalarpol(gun,varn(T));
983: if (l==0) return V;
984: V[2] = lcopy(x);
985: if (l==1) return V;
986: V[3] = (long) FpXQ_sqr(x,T,p);
987: for(i=4;i<l+2;i++)
988: V[i] = (long) FpXQ_mul((GEN) V[i-1],x,T,p);
989: #if 0
990: TODO: Karim proposes to use squaring:
991: V[i] = (long) ((i&1)?FpXQ_sqr((GEN) V[(i+1)>>1],T,p)
992: :FpXQ_mul((GEN) V[i-1],x,T,p));
993: Please profile it.
994: #endif
995: return V;
996: }
997: #if 0
998: static long brent_kung_nbmul(long d, long n, long p)
999: {
1000: return p+n*((d+p-1)/p);
1001: }
1002: TODO: This the the optimal parameter for the purpose of reducing
1003: multiplications, but profiling need to be done to ensure it is not slower
1004: than the other option in practice
1005: /*Return optimal parameter l for the evaluation of n polynomials of degree d*/
1006: long brent_kung_optpow(long d, long n)
1007: {
1008: double r;
1009: long f,c,pr;
1010: if (n>=d ) return d;
1011: pr=n*d;
1012: if (pr<=1) return 1;
1013: r=d/sqrt(pr);
1014: c=(long)ceil(d/ceil(r));
1015: f=(long)floor(d/floor(r));
1016: return (brent_kung_nbmul(d, n, c) <= brent_kung_nbmul(d, n, f))?c:f;
1017: }
1018: #endif
1019: /*Return optimal parameter l for the evaluation of n polynomials of degree d*/
1020: long brent_kung_optpow(long d, long n)
1021: {
1022: double r;
1023: long l,pr;
1024: if (n>=d ) return d;
1025: pr=n*d;
1026: if (pr<=1) return 1;
1027: r=d/sqrt(pr);
1028: l=(long)r;
1029: return (d+l-1)/l;
1030: }
1031:
1.2 ! noro 1032: /*Close to FpXV_FpV_innerprod*/
1.1 noro 1033:
1034: static GEN
1035: spec_compo_powers(GEN P, GEN V, long a, long n)
1036: {
1037: long i;
1038: GEN z;
1039: z = scalarpol((GEN)P[2+a],varn(P));
1040: for(i=1;i<=n;i++)
1041: z=FpX_add(z,FpX_Fp_mul((GEN)V[i+1],(GEN)P[2+a+i],NULL),NULL);
1042: return z;
1043: }
1044: /*Try to implement algorithm in Brent & Kung (Fast algorithms for
1045: *manipulating formal power series, JACM 25:581-595, 1978)
1046:
1047: V must be as output by FpXQ_powers.
1048: For optimal performance, l (of FpXQ_powers) must be as output by
1049: brent_kung_optpow
1050: */
1051:
1052: GEN
1053: FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p)
1054: {
1.2 ! noro 1055: gpmem_t ltop=avma;
1.1 noro 1056: long l=lg(V)-1;
1057: GEN z,u;
1058: long d=degpol(P),cnt=0;
1059: if (d==-1) return zeropol(varn(T));
1060: if (d<l)
1061: {
1062: z=spec_compo_powers(P,V,0,d);
1063: return gerepileupto(ltop,FpX_red(z,p));
1064: }
1065: if (l<=1)
1066: err(talker,"powers is only [] or [1] in FpX_FpXQV_compo");
1067: z=spec_compo_powers(P,V,d-l+1,l-1);
1068: d-=l;
1069: while(d>=l-1)
1070: {
1071: u=spec_compo_powers(P,V,d-l+2,l-2);
1072: z=FpX_add(u,FpXQ_mul(z,(GEN)V[l],T,p),NULL);
1073: d-=l-1;
1074: cnt++;
1075: }
1076: u=spec_compo_powers(P,V,0,d);
1077: z=FpX_add(u,FpXQ_mul(z,(GEN)V[d+2],T,p),NULL);
1078: cnt++;
1079: if (DEBUGLEVEL>=8) fprintferr("FpX_FpXQV_compo: %d FpXQ_mul [%d]\n",cnt,l-1);
1080: return gerepileupto(ltop,FpX_red(z,p));
1081: }
1082:
1083: /* T in Z[X] and x in Z/pZ[X]/(pol)
1084: * return lift(lift(subst(T,variable(T),Mod(x*Mod(1,p),pol*Mod(1,p)))));
1085: */
1086: GEN
1087: FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p)
1088: {
1.2 ! noro 1089: gpmem_t ltop=avma;
1.1 noro 1090: GEN z;
1091: long d=degpol(T),rtd;
1092: if (!signe(T)) return zeropol(varn(T));
1093: rtd = (long) sqrt((double)d);
1094: z = FpX_FpXQV_compo(T,FpXQ_powers(x,rtd,pol,p),pol,p);
1095: return gerepileupto(ltop,z);
1096: }
1097:
1098: /* Evaluation in Fp
1099: * x in Z[X] and y in Z return x(y) mod p
1100: *
1101: * If p is very large (several longs) and x has small coefficients(<<p),
1102: * then Brent & Kung algorithm is faster. */
1103: GEN
1104: FpX_eval(GEN x,GEN y,GEN p)
1105: {
1.2 ! noro 1106: gpmem_t av;
1.1 noro 1107: GEN p1,r,res;
1108: long i,j;
1109: i=lgef(x)-1;
1110: if (i<=2)
1111: return (i==2)? modii((GEN)x[2],p): gzero;
1112: res=cgetg(lgefint(p),t_INT);
1113: av=avma; p1=(GEN)x[i];
1114: /* specific attention to sparse polynomials (see poleval)*/
1115: /*You've guess it! It's a copy-paste(tm)*/
1116: for (i--; i>=2; i=j-1)
1117: {
1118: for (j=i; !signe((GEN)x[j]); j--)
1119: if (j==2)
1120: {
1.2 ! noro 1121: if (i!=j) y = powgumod(y,i-j+1,p);
1.1 noro 1122: p1=mulii(p1,y);
1123: goto fppoleval;/*sorry break(2) no implemented*/
1124: }
1.2 ! noro 1125: r = (i==j)? y: powgumod(y,i-j+1,p);
1.1 noro 1126: p1 = modii(addii(mulii(p1,r), (GEN)x[j]),p);
1127: }
1128: fppoleval:
1129: modiiz(p1,p,res);
1130: avma=av;
1131: return res;
1132: }
1133: /* Tz=Tx*Ty where Tx and Ty coprime
1134: * return lift(chinese(Mod(x*Mod(1,p),Tx*Mod(1,p)),Mod(y*Mod(1,p),Ty*Mod(1,p))))
1135: * if Tz is NULL it is computed
1136: * As we do not return it, and the caller will frequently need it,
1137: * it must compute it and pass it.
1138: */
1139: GEN
1140: FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p)
1141: {
1.2 ! noro 1142: gpmem_t av = avma;
1.1 noro 1143: GEN ax,p1;
1144: ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p);
1145: p1=FpX_mul(ax, FpX_sub(y,x,p),p);
1146: p1 = FpX_add(x,p1,p);
1147: if (!Tz) Tz=FpX_mul(Tx,Ty,p);
1148: p1 = FpX_res(p1,Tz,p);
1149: return gerepileupto(av,p1);
1150: }
1151:
1.2 ! noro 1152: typedef struct {
! 1153: GEN pol, p;
! 1154: } FpX_muldata;
! 1155: typedef struct {
! 1156: GEN pol;
! 1157: ulong p;
! 1158: } u_FpX_muldata;
! 1159:
! 1160: static GEN
! 1161: _u_sqr(void *data, GEN x)
! 1162: {
! 1163: u_FpX_muldata *D = (u_FpX_muldata*)data;
! 1164: return u_FpXQ_sqr(x, D->pol, D->p);
! 1165: }
! 1166: static GEN
! 1167: _u_mul(void *data, GEN x, GEN y)
! 1168: {
! 1169: u_FpX_muldata *D = (u_FpX_muldata*)data;
! 1170: return u_FpXQ_mul(x,y, D->pol, D->p);
! 1171: }
! 1172: static GEN
! 1173: _sqr(void *data, GEN x)
! 1174: {
! 1175: FpX_muldata *D = (FpX_muldata*)data;
! 1176: return FpXQ_sqr(x, D->pol, D->p);
! 1177: }
! 1178: static GEN
! 1179: _mul(void *data, GEN x, GEN y)
! 1180: {
! 1181: FpX_muldata *D = (FpX_muldata*)data;
! 1182: return FpXQ_mul(x,y, D->pol, D->p);
! 1183: }
! 1184:
1.1 noro 1185: /* assume n > 0 */
1186: GEN
1187: u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p)
1188: {
1.2 ! noro 1189: gpmem_t av = avma;
! 1190: u_FpX_muldata D;
! 1191: GEN y;
! 1192: D.pol = pol;
! 1193: D.p = p;
! 1194: y = leftright_pow(x, n, (void*)&D, &_u_sqr, &_u_mul);
1.1 noro 1195: return gerepileupto(av, y);
1196: }
1197:
1198: /* x,pol in Z[X], p in Z, n in Z, compute lift(x^n mod (p, pol)) */
1199: GEN
1200: FpXQ_pow(GEN x, GEN n, GEN pol, GEN p)
1201: {
1.2 ! noro 1202: gpmem_t av = avma;
! 1203: FpX_muldata D;
! 1204: long vx = varn(x);
! 1205: GEN y;
1.1 noro 1206: if (!signe(n)) return polun[vx];
1207: if (signe(n) < 0)
1208: {
1209: x=FpXQ_inv(x,pol,p);
1210: if (is_pm1(n)) return x; /* n = -1*/
1211: }
1212: else
1213: if (is_pm1(n)) return gcopy(x); /* n = 1 */
1214: if (OK_ULONG(p))
1215: {
1216: ulong pp = p[2];
1.2 ! noro 1217: pol = u_Fp_FpX(pol, pp);
! 1218: x = u_Fp_FpX(x, pp);
1.1 noro 1219: y = u_FpXQ_pow(x, n, pol, pp);
1220: y = small_to_pol(y,vx);
1221: }
1222: else
1223: {
1224: av = avma;
1.2 ! noro 1225: D.pol = pol;
! 1226: D.p = p;
! 1227: y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul);
1.1 noro 1228: }
1.2 ! noro 1229: return gerepileupto(av, y);
! 1230: }
! 1231:
! 1232: GEN
! 1233: Fq_pow(GEN x, GEN n, GEN pol, GEN p)
! 1234: {
! 1235: if (typ(x) == t_INT) return powmodulo(x,n,p);
! 1236: return FpXQ_pow(x,n,pol,p);
1.1 noro 1237: }
1238:
1239: /*******************************************************************/
1240: /* */
1241: /* Fp[X][Y] */
1242: /* */
1243: /*******************************************************************/
1244: /*Polynomials whose coefficients are either polynomials or integers*/
1245:
1246: GEN
1247: FpXX_red(GEN z, GEN p)
1248: {
1.2 ! noro 1249: GEN res;
! 1250: int i;
! 1251: res=cgetg(lgef(z),t_POL);
! 1252: res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z));
! 1253: for(i=2;i<lgef(res);i++)
! 1254: if (typ(z[i])!=t_INT)
! 1255: {
! 1256: gpmem_t av=avma;
! 1257: res[i]=(long)FpX_red((GEN)z[i],p);
! 1258: if (lgef(res[i])<=3)
1.1 noro 1259: {
1.2 ! noro 1260: if (lgef(res[i])==2) {avma=av;res[i]=zero;}
! 1261: else res[i]=lpilecopy(av,gmael(res,i,2));
1.1 noro 1262: }
1.2 ! noro 1263: }
! 1264: else
! 1265: res[i]=lmodii((GEN)z[i],p);
! 1266: res=normalizepol_i(res,lgef(res));
! 1267: return res;
1.1 noro 1268: }
1269:
1270: /*******************************************************************/
1271: /* */
1272: /* (Fp[X]/(Q))[Y] */
1273: /* */
1274: /*******************************************************************/
1.2 ! noro 1275: extern GEN to_Kronecker(GEN P, GEN Q);
1.1 noro 1276:
1277: /*Not malloc nor warn-clean.*/
1.2 ! noro 1278: GEN
! 1279: FpXQX_from_Kronecker(GEN Z, GEN T, GEN p)
1.1 noro 1280: {
1.2 ! noro 1281: long i,j,lx,l = lgef(Z), N = (degpol(T)<<1) + 1;
! 1282: GEN x, t = cgetg(N,t_POL), z = FpX_red(Z, p);
! 1283: t[1] = T[1] & VARNBITS;
! 1284: lx = (l-2) / (N-2);
! 1285: x = cgetg(lx+3,t_POL);
1.1 noro 1286: for (i=2; i<lx+2; i++)
1287: {
1288: for (j=2; j<N; j++) t[j] = z[j];
1289: z += (N-2);
1.2 ! noro 1290: x[i] = (long)FpX_res(normalizepol_i(t,N), T,p);
1.1 noro 1291: }
1292: N = (l-2) % (N-2) + 2;
1293: for (j=2; j<N; j++) t[j] = z[j];
1.2 ! noro 1294: x[i] = (long)FpX_res(normalizepol_i(t,N), T,p);
1.1 noro 1295: return normalizepol_i(x, i+1);
1296: }
1.2 ! noro 1297:
! 1298: GEN
! 1299: FpVQX_red(GEN z, GEN T, GEN p)
! 1300: {
! 1301: GEN res;
! 1302: int i, l = lg(z);
! 1303: res=cgetg(l,typ(z));
! 1304: for(i=1;i<l;i++)
! 1305: if (typ(z[i]) != t_INT)
! 1306: {
! 1307: if (T)
! 1308: res[i]=(long)FpX_res((GEN)z[i],T,p);
! 1309: else
! 1310: res[i]=(long)FpX_red((GEN)z[i],p);
! 1311: }
! 1312: else
! 1313: res[i] = lmodii((GEN)z[i],p);
! 1314: return res;
! 1315: }
1.1 noro 1316: GEN
1317: FpXQX_red(GEN z, GEN T, GEN p)
1318: {
1319: GEN res;
1.2 ! noro 1320: int i, l = lgef(z);
! 1321: res=cgetg(l,t_POL);
1.1 noro 1322: res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z));
1.2 ! noro 1323: for(i=2;i<l;i++)
! 1324: if (typ(z[i]) != t_INT)
1.1 noro 1325: {
1326: if (T)
1327: res[i]=(long)FpX_res((GEN)z[i],T,p);
1328: else
1329: res[i]=(long)FpX_red((GEN)z[i],p);
1330: }
1331: else
1332: res[i]=lmodii((GEN)z[i],p);
1333: res=normalizepol_i(res,lgef(res));
1334: return res;
1335: }
1336: /* if T = NULL, call FpX_mul */
1337: GEN
1338: FpXQX_mul(GEN x, GEN y, GEN T, GEN p)
1339: {
1.2 ! noro 1340: gpmem_t ltop;
1.1 noro 1341: GEN z,kx,ky;
1342: long vx;
1343: if (!T) return FpX_mul(x,y,p);
1344: ltop = avma;
1345: vx = min(varn(x),varn(y));
1346: kx= to_Kronecker(x,T);
1347: ky= to_Kronecker(y,T);
1348: z = quickmul(ky+2, kx+2, lgef(ky)-2, lgef(kx)-2);
1349: z = FpXQX_from_Kronecker(z,T,p);
1350: setvarn(z,vx);/*quickmul and Fq_from_Kronecker are not varn-clean*/
1351: return gerepileupto(ltop,z);
1352: }
1353: GEN/*Unused/untested*/
1354: FpXQX_sqr(GEN x, GEN T, GEN p)
1355: {
1.2 ! noro 1356: gpmem_t ltop=avma;
1.1 noro 1357: GEN z,kx;
1358: long vx=varn(x);
1359: kx= to_Kronecker(x,T);
1360: z = quicksqr(kx+2, lgef(kx)-2);
1361: z = FpXQX_from_Kronecker(z,T,p);
1362: setvarn(z,vx);/*quickmul and Fq_from_Kronecker are nor varn-clean*/
1363: return gerepileupto(ltop,z);
1364: }
1.2 ! noro 1365:
1.1 noro 1366: GEN
1367: FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p)
1368: {
1369: int i, lP = lgef(P);
1370: GEN res = cgetg(lP,t_POL);
1371: res[1] = evalsigne(1) | evalvarn(varn(P)) | evallgef(lP);
1372: for(i=2; i<lP; i++)
1.2 ! noro 1373: res[i] = (long)Fq_mul(U,(GEN)P[i], T,p);
1.1 noro 1374: return normalizepol_i(res,lgef(res));
1375: }
1376:
1377: /* a X^degpol, assume degpol >= 0 */
1378: static GEN
1379: monomial(GEN a, int degpol, int v)
1380: {
1381: long i, lP = degpol+3;
1382: GEN P = cgetg(lP, t_POL);
1383: P[1] = evalsigne(1) | evalvarn(v) | evallgef(lP);
1384: lP--; P[lP] = (long)a;
1385: for (i=2; i<lP; i++) P[i] = zero;
1386: return P;
1387: }
1388:
1389: /* return NULL if Euclidean remainder sequence fails (==> T reducible mod p)
1390: * last non-zero remainder otherwise */
1391: GEN
1392: FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p)
1393: {
1.2 ! noro 1394: gpmem_t btop, ltop = avma, st_lim;
1.1 noro 1395: long dg, vx = varn(P);
1396: GEN U, q;
1.2 ! noro 1397: P = FpXX_red(P, p); btop = avma;
1.1 noro 1398: Q = FpXX_red(Q, p);
1399: if (!signe(P)) return gerepileupto(ltop, Q);
1.2 ! noro 1400: if (!signe(Q)) { avma = btop; return P; }
1.1 noro 1401: T = FpX_red(T, p);
1402:
1403: btop = avma; st_lim = stack_lim(btop, 1);
1404: dg = lgef(P)-lgef(Q);
1405: if (dg < 0) { swap(P, Q); dg = -dg; }
1406: for(;;)
1407: {
1408: U = FpXQ_invsafe(leading_term(Q), T, p);
1409: if (!U) { avma = ltop; return NULL; }
1410: do /* set P := P % Q */
1411: {
1.2 ! noro 1412: q = Fq_mul(U, gneg(leading_term(P)), T, p);
1.1 noro 1413: P = gadd(P, FpXQX_mul(monomial(q, dg, vx), Q, T, p));
1414: P = FpXQX_red(P, T, p); /* wasteful, but negligible */
1415: dg = lgef(P)-lgef(Q);
1416: } while (dg >= 0);
1417: if (!signe(P)) break;
1418:
1419: if (low_stack(st_lim, stack_lim(btop, 1)))
1420: {
1421: GEN *bptr[2]; bptr[0]=&P; bptr[1]=&Q;
1422: if (DEBUGLEVEL>1) err(warnmem,"FpXQX_safegcd");
1423: gerepilemany(btop, bptr, 2);
1424: }
1425: swap(P, Q); dg = -dg;
1426: }
1427: Q = FpXQX_FpXQ_mul(Q, U, T, p); /* normalize GCD */
1428: return gerepileupto(ltop, Q);
1429: }
1430:
1431: /*******************************************************************/
1432: /* */
1.2 ! noro 1433: /* (Fp[X]/T(X))[Y] / S(Y) */
! 1434: /* */
! 1435: /*******************************************************************/
! 1436: /* TODO: merge these 2 (I don't understand the first one) -- KB */
! 1437:
! 1438: /*Preliminary implementation to speed up Fp_isom*/
! 1439: typedef struct {
! 1440: GEN S, T, p;
! 1441: } FpXQYQ_muldata;
! 1442:
! 1443: static GEN
! 1444: FpXQYQ_redswap(GEN x, GEN S, GEN T, GEN p)
! 1445: {
! 1446: gpmem_t ltop=avma;
! 1447: long n=degpol(S);
! 1448: long m=degpol(T);
! 1449: long v=varn(T),w=varn(S);
! 1450: GEN V = swap_polpol(x,n,w);
! 1451: setvarn(T,w);
! 1452: V = FpXQX_red(V,T,p);
! 1453: setvarn(T,v);
! 1454: V = swap_polpol(V,m,w);
! 1455: return gerepilecopy(ltop,V);
! 1456: }
! 1457: static GEN
! 1458: FpXQYQ_sqr(void *data, GEN x)
! 1459: {
! 1460: FpXQYQ_muldata *D = (FpXQYQ_muldata*)data;
! 1461: return FpXQYQ_redswap(FpXQX_sqr(x, D->S, D->p),D->S,D->T,D->p);
! 1462:
! 1463: }
! 1464: static GEN
! 1465: FpXQYQ_mul(void *data, GEN x, GEN y)
! 1466: {
! 1467: FpXQYQ_muldata *D = (FpXQYQ_muldata*)data;
! 1468: return FpXQYQ_redswap(FpXQX_mul(x,y, D->S, D->p),D->S,D->T,D->p);
! 1469: }
! 1470:
! 1471: /* x in Z[X,Y], S in Z[X] over Fq = Z[Y]/(p,T); compute lift(x^n mod (S,T,p)) */
! 1472: GEN
! 1473: FpXQYQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
! 1474: {
! 1475: gpmem_t av = avma;
! 1476: FpXQYQ_muldata D;
! 1477: GEN y;
! 1478: D.S = S;
! 1479: D.T = T;
! 1480: D.p = p;
! 1481: y = leftright_pow(x, n, (void*)&D, &FpXQYQ_sqr, &FpXQYQ_mul);
! 1482: return gerepileupto(av, y);
! 1483: }
! 1484:
! 1485: typedef struct {
! 1486: GEN T, p, S;
! 1487: long v;
! 1488: } kronecker_muldata;
! 1489:
! 1490: static GEN
! 1491: _FqXQ_red(void *data, GEN x)
! 1492: {
! 1493: kronecker_muldata *D = (kronecker_muldata*)data;
! 1494: GEN t = FpXQX_from_Kronecker(x, D->T,D->p);
! 1495: setvarn(t, D->v);
! 1496: t = FpXQX_divres(t, D->S,D->T,D->p, ONLY_REM);
! 1497: return to_Kronecker(t,D->T);
! 1498: }
! 1499: static GEN
! 1500: _FqXQ_mul(void *data, GEN x, GEN y) {
! 1501: return _FqXQ_red(data, FpX_mul(x,y,NULL));
! 1502: }
! 1503: static GEN
! 1504: _FqXQ_sqr(void *data, GEN x) {
! 1505: return _FqXQ_red(data, FpX_sqr(x,NULL));
! 1506: }
! 1507:
! 1508: /* x over Fq, return lift(x^n) mod S */
! 1509: GEN
! 1510: FqXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
! 1511: {
! 1512: gpmem_t av0 = avma;
! 1513: long v = varn(x);
! 1514: GEN y;
! 1515: kronecker_muldata D;
! 1516:
! 1517: x = to_Kronecker(x,T);
! 1518: D.S = S;
! 1519: D.T = T;
! 1520: D.p = p;
! 1521: D.v = v;
! 1522: y = leftright_pow(x, n, (void*)&D, &_FqXQ_sqr, &_FqXQ_mul);
! 1523: y = FpXQX_from_Kronecker(y, T,p);
! 1524: setvarn(y, v); return gerepileupto(av0, y);
! 1525: }
! 1526: /*******************************************************************/
! 1527: /* */
1.1 noro 1528: /* Fq[X] */
1529: /* */
1530: /*******************************************************************/
1531:
1532: /*Well nothing, for now. So we reuse FpXQX*/
1533: #define FqX_mul FpXQX_mul
1534: #define FqX_sqr FpXQX_sqr
1535: #define FqX_red FpXQX_red
1536: static GEN
1.2 ! noro 1537: Fq_neg(GEN x, GEN T/*unused*/, GEN p)
1.1 noro 1538: {
1539: switch(typ(x)==t_POL)
1540: {
1541: case 0: return signe(x)?subii(p,x):gzero;
1542: case 1: return FpX_neg(x,p);
1543: }
1544: return NULL;
1545: }
1546: static GEN modulo,Tmodulo;
1547: static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodulo,modulo);}
1548: GEN
1549: FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v)
1550: {
1.2 ! noro 1551: gpmem_t ltop=avma;
1.1 noro 1552: long k;
1553: GEN W=cgetg(lg(V),t_VEC);
1554: for(k=1;k<lg(V);k++)
1555: W[k]=(long)deg1pol(gun,Fq_neg((GEN)V[k],Tp,p),v);
1556: modulo=p;Tmodulo=Tp;
1557: W=divide_conquer_prod(W,&fgmul);
1558: return gerepileupto(ltop,W);
1559: }
1560:
1561: /*******************************************************************/
1562: /* */
1563: /* n-th ROOT in Fq */
1564: /* */
1565: /*******************************************************************/
1566: /*NO clean malloc*/
1567: static GEN fflgen(GEN l, long e, GEN r, GEN T ,GEN p, GEN *zeta)
1568: {
1.2 ! noro 1569: gpmem_t av1 = avma;
! 1570: GEN z, m, m1;
! 1571: const long pp = is_bigint(p)? VERYBIGINT: itos(p);
! 1572: long x=varn(T),k,u,v,i;
! 1573:
! 1574: for (k=0; ; k++)
! 1575: {
! 1576: z = (degpol(T)==1)? polun[x]: polx[x];
! 1577: u = k/pp; v=2; /* FpX_Fp_add modify y */
! 1578: z = gaddgs(z, k%pp);
! 1579: while(u)
1.1 noro 1580: {
1.2 ! noro 1581: z = FpX_add(z, monomial(stoi(u%pp),v,x), NULL);
! 1582: u /= pp; v++;
1.1 noro 1583: }
1.2 ! noro 1584: if ( DEBUGLEVEL>=6 ) fprintferr("FF l-Gen:next %Z\n",z);
1.1 noro 1585: m1 = m = FpXQ_pow(z,r,T,p);
1.2 ! noro 1586: if (gcmp1(m)) { avma = av1; continue; }
! 1587:
1.1 noro 1588: for (i=1; i<e; i++)
1.2 ! noro 1589: if (gcmp1(m = FpXQ_pow(m,l,T,p))) break;
1.1 noro 1590: if (i==e) break;
1591: avma = av1;
1592: }
1.2 ! noro 1593: *zeta = m; return m1;
1.1 noro 1594: }
1.2 ! noro 1595: /* Solve x^l = a mod (p,T)
! 1596: * l must be prime
! 1597: * q = p^degpol(T)-1 = (l^e)*r, with e>=1 and pgcd(r,l)=1
! 1598: * m = y^(q/l)
! 1599: * y not an l-th power [ m != 1 ] */
1.1 noro 1600: GEN
1601: ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, GEN r, GEN y, GEN m)
1602: {
1.2 ! noro 1603: gpmem_t av = avma, lim;
1.1 noro 1604: long i,k;
1605: GEN p1,p2,u1,u2,v,w,z;
1606:
1.2 ! noro 1607: if (gcmp1(a)) return gcopy(a);
! 1608:
! 1609: (void)bezout(r,l,&u1,&u2); /* result is 1 */
! 1610: v = FpXQ_pow(a,u2,T,p);
! 1611: w = FpXQ_pow(a, modii(mulii(negi(u1),r),q), T,p);
1.1 noro 1612: lim = stack_lim(av,1);
1613: while (!gcmp1(w))
1614: {
1.2 ! noro 1615: k = 0;
! 1616: p1 = w;
! 1617: do { /* if p is not prime, loop will not end */
! 1618: z = p1;
! 1619: p1 = FpXQ_pow(p1,l,T,p);
1.1 noro 1620: k++;
1.2 ! noro 1621: } while (!gcmp1(p1));
1.1 noro 1622: if (k==e) { avma=av; return NULL; }
1623: p2 = FpXQ_mul(z,m,T,p);
1.2 ! noro 1624: for (i=1; !gcmp1(p2); i++) p2 = FpXQ_mul(p2,m,T,p);/*TODO: BS/GS instead*/
! 1625: p1= FpXQ_pow(y, modii(mulsi(i,gpowgs(l,e-k-1)), q), T,p);
1.1 noro 1626: m = FpXQ_pow(m,stoi(i),T,p);
1627: e = k;
1628: v = FpXQ_mul(p1,v,T,p);
1629: y = FpXQ_pow(p1,l,T,p);
1630: w = FpXQ_mul(y,w,T,p);
1631: if (low_stack(lim, stack_lim(av,1)))
1632: {
1633: if(DEBUGMEM>1) err(warnmem,"ffsqrtlmod");
1.2 ! noro 1634: gerepileall(av,4, &y,&v,&w,&m);
1.1 noro 1635: }
1636: }
1.2 ! noro 1637: return gerepilecopy(av, v);
1.1 noro 1638: }
1.2 ! noro 1639: /* Solve x^n = a mod p: n integer, a in Fp[X]/(T) [ p prime, T irred. mod p ]
! 1640: *
! 1641: * 1) if no solution, return NULL and (if zetan != NULL) set zetan to gzero.
! 1642: *
! 1643: * 2) If there is a solution, there are exactly m=gcd(p-1,n) of them.
! 1644: * If zetan != NULL, it is set to a primitive mth root of unity so that the set
! 1645: * of solutions is {x*zetan^k;k=0 to m-1}
! 1646: *
! 1647: * If a = 0, return 0 and (if zetan != NULL) set zetan = gun */
1.1 noro 1648: GEN ffsqrtnmod(GEN a, GEN n, GEN T, GEN p, GEN *zetan)
1649: {
1.2 ! noro 1650: gpmem_t ltop=avma, av1, lim;
1.1 noro 1651: long i,j,e;
1652: GEN m,u1,u2;
1653: GEN q,r,zeta,y,l,z;
1.2 ! noro 1654:
1.1 noro 1655: if (typ(a) != t_POL || typ(n) != t_INT || typ(T) != t_POL || typ(p)!=t_INT)
1656: err(typeer,"ffsqrtnmod");
1.2 ! noro 1657: if (!degpol(T)) err(constpoler,"ffsqrtnmod");
! 1658: if (!signe(n)) err(talker,"1/0 exponent in ffsqrtnmod");
! 1659: if (gcmp1(n)) {if (zetan) *zetan=gun;return gcopy(a);}
! 1660: if (gcmp0(a)) {if (zetan) *zetan=gun;return gzero;}
! 1661:
! 1662: q = addsi(-1, gpowgs(p,degpol(T)));
! 1663: m = bezout(n,q,&u1,&u2);
! 1664: if (!egalii(m,n)) a = FpXQ_pow(a, modii(u1,q), T,p);
! 1665: if (zetan) z = polun[varn(T)];
! 1666: lim = stack_lim(ltop,1);
1.1 noro 1667: if (!gcmp1(m))
1668: {
1.2 ! noro 1669: m = decomp(m); av1 = avma;
1.1 noro 1670: for (i = lg(m[1])-1; i; i--)
1671: {
1.2 ! noro 1672: l = gcoeff(m,i,1);
! 1673: j = itos(gcoeff(m,i,2));
! 1674: e = pvaluation(q,l,&r);
! 1675: if(DEBUGLEVEL>=6) (void)timer2();
! 1676: y = fflgen(l,e,r,T,p,&zeta);
! 1677: if(DEBUGLEVEL>=6) msgtimer("fflgen");
! 1678: if (zetan) z = FpXQ_mul(z, FpXQ_pow(y,gpowgs(l,e-j),T,p), T,p);
! 1679: for (; j; j--)
1.1 noro 1680: {
1.2 ! noro 1681: a = ffsqrtlmod(a,l,T,p,q,e,r,y,zeta);
! 1682: if (!a) {avma=ltop; return NULL;}
! 1683: }
1.1 noro 1684: if (low_stack(lim, stack_lim(ltop,1)))
1.2 ! noro 1685: { /* n can have lots of prime factors */
1.1 noro 1686: if(DEBUGMEM>1) err(warnmem,"ffsqrtnmod");
1.2 ! noro 1687: gerepileall(av1,zetan? 2: 1, &a,&z);
1.1 noro 1688: }
1689: }
1690: }
1691: if (zetan)
1692: {
1.2 ! noro 1693: *zetan = z;
! 1694: gerepileall(ltop,2,&a,zetan);
1.1 noro 1695: }
1696: else
1.2 ! noro 1697: a = gerepileupto(ltop, a);
1.1 noro 1698: return a;
1699: }
1700: /*******************************************************************/
1701: /* Isomorphisms between finite fields */
1702: /* */
1703: /*******************************************************************/
1704: GEN
1705: matrixpow(long n, long m, GEN y, GEN P,GEN l)
1706: {
1707: return vecpol_to_mat(FpXQ_powers(y,m-1,P,l),n);
1708: }
1709: /* compute the reciprocical isomorphism of S mod T,p, i.e. V such that
1710: V(S)=X mod T,p*/
1711: GEN
1712: Fp_inv_isom(GEN S,GEN T, GEN p)
1713: {
1.2 ! noro 1714: gpmem_t lbot, ltop = avma;
1.1 noro 1715: GEN M, V;
1716: int n, i;
1717: long x;
1718: x = varn(T);
1719: n = degpol(T);
1720: M = matrixpow(n,n,S,T,p);
1721: V = cgetg(n + 1, t_COL);
1722: for (i = 1; i <= n; i++)
1723: V[i] = zero;
1724: V[2] = un;
1725: V = FpM_invimage(M,V,p);
1726: lbot = avma;
1727: V = gtopolyrev(V, x);
1728: return gerepile(ltop, lbot, V);
1729: }
1.2 ! noro 1730:
! 1731: /* Let M the matrix of the x^p Frobenius automorphism.
! 1732: * Compute x^(p^i) for i=0..r
1.1 noro 1733: */
1734: static GEN
1735: polfrobenius(GEN M, GEN p, long r, long v)
1736: {
1737: GEN V,W;
1738: long i;
1739: V = cgetg(r+2,t_VEC);
1740: V[1] = (long) polx[v];
1741: if (r == 0) return V;
1.2 ! noro 1742: V[2] = (long) vec_to_pol((GEN)M[2],v);
1.1 noro 1743: W = (GEN) M[2];
1744: for (i = 3; i <= r+1; ++i)
1745: {
1746: W = FpM_FpV_mul(M,W,p);
1.2 ! noro 1747: V[i] = (long) vec_to_pol(W,v);
1.1 noro 1748: }
1749: return V;
1750: }
1751:
1.2 ! noro 1752: /* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in
! 1753: * FFp[X]/T. Compute P(M)
! 1754: * V=polfrobenius(M, p, degpol(P), v)
! 1755: * not stack clean
! 1756: */
! 1757:
1.1 noro 1758: static GEN
1759: matpolfrobenius(GEN V, GEN P, GEN T, GEN p)
1760: {
1.2 ! noro 1761: gpmem_t btop;
1.1 noro 1762: long i;
1763: long l=degpol(T);
1764: long v=varn(T);
1.2 ! noro 1765: GEN M,W,Mi;
1.1 noro 1766: GEN PV=gtovec(P);
1.2 ! noro 1767: GEN *gptr[2];
! 1768: long lV=lg(V);
1.1 noro 1769: PV=cgetg(degpol(P)+2,t_VEC);
1770: for(i=1;i<lg(PV);i++)
1771: PV[i]=P[1+i];
1772: M=cgetg(l+1,t_VEC);
1773: M[1]=(long)scalarpol(poleval(P,gun),v);
1.2 ! noro 1774: M[2]=(long)FpXV_FpV_innerprod(V,PV,p);
! 1775: btop=avma;
! 1776: gptr[0]=&Mi;
! 1777: gptr[1]=&W;
! 1778: W=cgetg(lV,t_VEC);
! 1779: for(i=1;i<lV;i++)
1.1 noro 1780: W[i]=V[i];
1781: for(i=3;i<=l;i++)
1782: {
1783: long j;
1.2 ! noro 1784: gpmem_t bbot;
! 1785: GEN W2=cgetg(lV,t_VEC);
! 1786: for(j=1;j<lV;j++)
! 1787: W2[j]=(long)FpXQ_mul((GEN)W[j],(GEN)V[j],T,p);
! 1788: bbot=avma;
! 1789: Mi=FpXV_FpV_innerprod(W2,PV,p);
! 1790: W=gcopy(W2);
! 1791: gerepilemanysp(btop,bbot,gptr,2);
! 1792: btop=(gpmem_t)W;
! 1793: M[i]=(long)Mi;
1.1 noro 1794: }
1795: return vecpol_to_mat(M,l);
1796: }
1797:
1.2 ! noro 1798: /* Let M the matrix of the Frobenius automorphism of Fp[X]/(T).
! 1799: * Compute M^d
! 1800: * TODO: use left-right binary (tricky!)
! 1801: */
! 1802: GEN
! 1803: FpM_frobenius_pow(GEN M, long d, GEN T, GEN p)
! 1804: {
! 1805: gpmem_t ltop=avma;
! 1806: long i,l=degpol(T);
! 1807: GEN R, W = (GEN) M[2];
! 1808: for (i = 2; i <= d; ++i)
! 1809: W = FpM_FpV_mul(M,W,p);
! 1810: R=matrixpow(l,l,vec_to_pol(W,varn(T)),T,p);
! 1811: return gerepilecopy(ltop,R);
! 1812: }
! 1813:
1.1 noro 1814: /* Essentially we want to compute
1.2 ! noro 1815: * FqM_ker(MA-polx[MAXVARN],U,l)
1.1 noro 1816: * To avoid use of matrix in Fq we procede as follows:
1.2 ! noro 1817: * We compute FpM_ker(U(MA),l) and then we recover
1.1 noro 1818: * the eigen value by Galois action, see formula.
1819: */
1820: static GEN
1.2 ! noro 1821: intersect_ker(GEN P, GEN MA, GEN U, GEN l)
1.1 noro 1822: {
1.2 ! noro 1823: gpmem_t ltop=avma;
1.1 noro 1824: long vp=varn(P);
1.2 ! noro 1825: long vu=varn(U), r=degpol(U);
1.1 noro 1826: long i;
1827: GEN A, R, M, ib0, V;
1.2 ! noro 1828: if (DEBUGLEVEL>=4) (void)timer2();
! 1829: V=polfrobenius(MA,l,r,vu);
1.1 noro 1830: if (DEBUGLEVEL>=4) msgtimer("pol[frobenius]");
1.2 ! noro 1831: M=matpolfrobenius(V,U,P,l);
1.1 noro 1832: if (DEBUGLEVEL>=4) msgtimer("matrix cyclo");
1833: A=FpM_ker(M,l);
1834: if (DEBUGLEVEL>=4) msgtimer("kernel");
1835: A=gerepileupto(ltop,A);
1836: if (lg(A)!=r+1)
1837: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
1838: ,l,polx[vp],P);
1839: /*The formula is
1840: * a_{r-1}=-\phi(a_0)/b_0
1841: * a_{i-1}=\phi(a_i)+b_ia_{r-1} i=r-1 to 1
1842: * Where a_0=A[1] and b_i=U[i+2]
1843: */
1.2 ! noro 1844: ib0=negi(mpinvmod((GEN)U[2],l));
1.1 noro 1845: R=cgetg(r+1,t_MAT);
1846: R[1]=A[1];
1847: R[r]=(long)FpM_FpV_mul(MA,gmul((GEN)A[1],ib0),l);
1848: for(i=r-1;i>1;i--)
1849: R[i]=(long)FpV_red(gadd(FpM_FpV_mul(MA,(GEN)R[i+1],l),
1.2 ! noro 1850: gmul((GEN)U[i+2],(GEN)R[r])),l);
1.1 noro 1851: R=gtrans_i(R);
1852: for(i=1;i<lg(R);i++)
1.2 ! noro 1853: R[i]=(long)vec_to_pol((GEN)R[i],vu);
1.1 noro 1854: A=gtopolyrev(R,vp);
1855: return gerepileupto(ltop,A);
1856: }
1857:
1858: /* n must divide both the degree of P and Q. Compute SP and SQ such
1859: that the subfield of FF_l[X]/(P) generated by SP and the subfield of
1860: FF_l[X]/(Q) generated by SQ are isomorphic of degree n. P and Q do
1861: not need to be of the same variable. if MA (resp. MB) is not NULL,
1862: must be the matrix of the frobenius map in FF_l[X]/(P) (resp.
1863: FF_l[X]/(Q) ). */
1864: /* Note on the implementation choice:
1865: * We assume the prime p is very large
1866: * so we handle Frobenius as matrices.
1867: */
1868: void
1869: Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN *SQ, GEN MA, GEN MB)
1870: {
1.2 ! noro 1871: gpmem_t lbot, ltop=avma;
1.1 noro 1872: long vp,vq,np,nq,e,pg;
1873: GEN q;
1874: GEN A,B,Ap,Bp;
1875: GEN *gptr[2];
1876: vp=varn(P);vq=varn(Q);
1877: np=degpol(P);nq=degpol(Q);
1878: if (np<=0 || nq<=0 || n<=0 || np%n!=0 || nq%n!=0)
1879: err(talker,"bad degrees in Fp_intersect: %d,%d,%d",n,degpol(P),degpol(Q));
1880: e=pvaluation(stoi(n),l,&q);
1881: pg=itos(q);
1882: avma=ltop;
1883: if(!MA) MA=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l);
1884: if(!MB) MB=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l);
1885: A=Ap=zeropol(vp);
1886: B=Bp=zeropol(vq);
1.2 ! noro 1887: if (pg > 1)
1.1 noro 1888: {
1.2 ! noro 1889: if (smodis(l,pg) == 1)
1.1 noro 1890: /*We do not need to use relative extension in this setting, so
1891: we don't. (Well,now that we don't in the other case also, it is more
1892: dubious to treat cases apart...)*/
1893: {
1894: GEN L,An,Bn,ipg,z;
1895: z=rootmod(cyclo(pg,-1),l);
1896: if (lg(z)<2) err(talker,"%Z is not a prime in Fp_intersect",l);
1897: z=negi(lift((GEN)z[1]));
1898: ipg=stoi(pg);
1.2 ! noro 1899: if (DEBUGLEVEL>=4) (void)timer2();
1.1 noro 1900: A=FpM_ker(gaddmat(z, MA),l);
1901: if (lg(A)!=2)
1902: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
1903: ,l,polx[vp],P);
1.2 ! noro 1904: A=vec_to_pol((GEN)A[1],vp);
1.1 noro 1905: B=FpM_ker(gaddmat(z, MB),l);
1906: if (lg(B)!=2)
1907: err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect"
1908: ,l,polx[vq],Q);
1.2 ! noro 1909: B=vec_to_pol((GEN)B[1],vq);
1.1 noro 1910: if (DEBUGLEVEL>=4) msgtimer("FpM_ker");
1911: An=(GEN) FpXQ_pow(A,ipg,P,l)[2];
1912: Bn=(GEN) FpXQ_pow(B,ipg,Q,l)[2];
1.2 ! noro 1913: if (!invmod(Bn,l,&z))
! 1914: err(talker,"Polynomials not irreducible in Fp_intersect");
! 1915: z=modii(mulii(An,z),l);
1.1 noro 1916: L=mpsqrtnmod(z,ipg,l,NULL);
1917: if ( !L )
1918: err(talker,"Polynomials not irreducible in Fp_intersect");
1919: if (DEBUGLEVEL>=4) msgtimer("mpsqrtnmod");
1920: B=FpX_Fp_mul(B,L,l);
1921: }
1922: else
1923: {
1.2 ! noro 1924: GEN L,An,Bn,ipg,U,z;
! 1925: U=lift(gmael(factmod(cyclo(pg,MAXVARN),l),1,1));
1.1 noro 1926: ipg=stoi(pg);
1.2 ! noro 1927: A=intersect_ker(P, MA, U, l);
! 1928: B=intersect_ker(Q, MB, U, l);
! 1929: if (DEBUGLEVEL>=4) (void)timer2();
! 1930: An=(GEN) FpXQYQ_pow(A,stoi(pg),U,P,l)[2];
! 1931: Bn=(GEN) FpXQYQ_pow(B,stoi(pg),U,Q,l)[2];
1.1 noro 1932: if (DEBUGLEVEL>=4) msgtimer("pows [P,Q]");
1.2 ! noro 1933: z=FpXQ_inv(Bn,U,l);
! 1934: z=FpXQ_mul(An,z,U,l);
! 1935: L=ffsqrtnmod(z,ipg,U,l,NULL);
1.1 noro 1936: if (DEBUGLEVEL>=4) msgtimer("ffsqrtn");
1937: if ( !L )
1938: err(talker,"Polynomials not irreducible in Fp_intersect");
1.2 ! noro 1939: B=FpXQX_FpXQ_mul(B,L,U,l);
! 1940: B=gsubst(B,MAXVARN,gzero);
! 1941: A=gsubst(A,MAXVARN,gzero);
1.1 noro 1942: }
1943: }
1944: if (e!=0)
1945: {
1946: GEN VP,VQ,moinsun,Ay,By,lmun;
1947: int i,j;
1948: moinsun=stoi(-1);
1949: lmun=addis(l,-1);
1950: MA=gaddmat(moinsun,MA);
1951: MB=gaddmat(moinsun,MB);
1952: Ay=polun[vp];
1953: By=polun[vq];
1954: VP=cgetg(np+1,t_COL);
1955: VP[1]=un;
1956: for(i=2;i<=np;i++) VP[i]=zero;
1957: if (np==nq) VQ=VP;/*save memory*/
1958: else
1959: {
1960: VQ=cgetg(nq+1,t_COL);
1961: VQ[1]=un;
1962: for(i=2;i<=nq;i++) VQ[i]=zero;
1963: }
1964: for(j=0;j<e;j++)
1965: {
1966: if (j)
1967: {
1968: Ay=FpXQ_mul(Ay,FpXQ_pow(Ap,lmun,P,l),P,l);
1969: for(i=1;i<lgef(Ay)-1;i++) VP[i]=Ay[i+1];
1970: for(;i<=np;i++) VP[i]=zero;
1971: }
1972: Ap=FpM_invimage(MA,VP,l);
1.2 ! noro 1973: Ap=vec_to_pol(Ap,vp);
1.1 noro 1974: if (j)
1975: {
1976: By=FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l);
1977: for(i=1;i<lgef(By)-1;i++) VQ[i]=By[i+1];
1978: for(;i<=nq;i++) VQ[i]=zero;
1979: }
1980: Bp=FpM_invimage(MB,VQ,l);
1.2 ! noro 1981: Bp=vec_to_pol(Bp,vq);
1.1 noro 1982: if (DEBUGLEVEL>=4) msgtimer("FpM_invimage");
1983: }
1984: }/*FpX_add is not clean, so we must do it *before* lbot=avma*/
1985: A=FpX_add(A,Ap,NULL);
1986: B=FpX_add(B,Bp,NULL);
1987: lbot=avma;
1988: *SP=FpX_red(A,l);
1989: *SQ=FpX_red(B,l);
1990: gptr[0]=SP;gptr[1]=SQ;
1991: gerepilemanysp(ltop,lbot,gptr,2);
1992: }
1993: /* Let l be a prime number, P, Q in ZZ[X]. P and Q are both
1994: * irreducible modulo l and degree(P) divides degree(Q). Output a
1995: * monomorphism between FF_l[X]/(P) and FF_l[X]/(Q) as a polynomial R such
1996: * that Q | P(R) mod l. If P and Q have the same degree, it is of course an
1997: * isomorphism. */
1998: GEN Fp_isom(GEN P,GEN Q,GEN l)
1999: {
1.2 ! noro 2000: gpmem_t ltop=avma;
1.1 noro 2001: GEN SP,SQ,R;
2002: Fp_intersect(degpol(P),P,Q,l,&SP,&SQ,NULL,NULL);
2003: R=Fp_inv_isom(SP,P,l);
2004: R=FpX_FpXQ_compo(R,SQ,Q,l);
2005: return gerepileupto(ltop,R);
2006: }
2007: GEN
1.2 ! noro 2008: Fp_factorgalois(GEN P,GEN l, long d, long w, GEN MP)
1.1 noro 2009: {
1.2 ! noro 2010: gpmem_t ltop=avma;
! 2011: GEN R,V,Tl,z,M;
1.1 noro 2012: long n,k,m;
2013: long v;
2014: v=varn(P);
2015: n=degpol(P);
2016: m=n/d;
1.2 ! noro 2017: if (DEBUGLEVEL>=4) (void)timer2();
! 2018: M=FpM_frobenius_pow(MP,d,P,l);
! 2019: if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: frobenius power");
1.1 noro 2020: Tl=gcopy(P); setvarn(Tl,w);
2021: V=cgetg(m+1,t_VEC);
2022: V[1]=lpolx[w];
1.2 ! noro 2023: z=pol_to_vec((GEN)V[1],n);
1.1 noro 2024: for(k=2;k<=m;k++)
1.2 ! noro 2025: {
! 2026: z=FpM_FpV_mul(M,z,l);
! 2027: V[k]=(long)vec_to_pol(z,w);
! 2028: }
! 2029: if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: roots");
1.1 noro 2030: R=FqV_roots_to_pol(V,Tl,l,v);
1.2 ! noro 2031: if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: pol");
1.1 noro 2032: return gerepileupto(ltop,R);
2033: }
2034: /* P,Q irreducible over F_l. Factor P over FF_l[X] / Q [factors are ordered as
2035: * a Frobenius cycle] */
2036: GEN
2037: Fp_factor_irred(GEN P,GEN l, GEN Q)
2038: {
1.2 ! noro 2039: gpmem_t ltop=avma, av;
! 2040: GEN SP,SQ,MP,MQ,M,FP,FQ,E,V,IR,res;
1.1 noro 2041: long np=degpol(P),nq=degpol(Q);
2042: long i,d=cgcd(np,nq);
2043: long vp=varn(P),vq=varn(Q);
2044: if (d==1)
2045: {
2046: res=cgetg(2,t_COL);
2047: res[1]=lcopy(P);
2048: return res;
2049: }
1.2 ! noro 2050: if (DEBUGLEVEL>=4) (void)timer2();
! 2051: FP=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l);
! 2052: FQ=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l);
! 2053: if (DEBUGLEVEL>=4) msgtimer("matrixpows");
! 2054: Fp_intersect(d,P,Q,l,&SP,&SQ,FP,FQ);
1.1 noro 2055: av=avma;
1.2 ! noro 2056: E=Fp_factorgalois(P,l,d,vq,FP);
1.1 noro 2057: E=polpol_to_mat(E,np);
2058: MP = matrixpow(np,d,SP,P,l);
2059: IR = (GEN)FpM_sindexrank(MP,l)[1];
2060: E = rowextract_p(E, IR);
2061: M = rowextract_p(MP,IR);
2062: M = FpM_inv(M,l);
2063: MQ = matrixpow(nq,d,SQ,Q,l);
2064: M = FpM_mul(MQ,M,l);
2065: M = FpM_mul(M,E,l);
2066: M = gerepileupto(av,M);
2067: V = cgetg(d+1,t_VEC);
2068: V[1]=(long)M;
2069: for(i=2;i<=d;i++)
1.2 ! noro 2070: V[i]=(long)FpM_mul(FQ,(GEN)V[i-1],l);
1.1 noro 2071: res=cgetg(d+1,t_COL);
2072: for(i=1;i<=d;i++)
2073: res[i]=(long)mat_to_polpol((GEN)V[i],vp,vq);
2074: return gerepilecopy(ltop,res);
2075: }
2076: GEN Fp_factor_rel0(GEN P,GEN l, GEN Q)
2077: {
1.2 ! noro 2078: gpmem_t ltop=avma, tetpil;
1.1 noro 2079: GEN V,ex,F,y,R;
2080: long n,nbfact=0,nmax=lgef(P)-2;
2081: long i;
2082: F=factmod0(P,l);
2083: n=lg((GEN)F[1]);
2084: V=cgetg(nmax,t_VEC);
2085: ex=cgetg(nmax,t_VECSMALL);
2086: for(i=1;i<n;i++)
2087: {
2088: int j,r;
2089: R=Fp_factor_irred(gmael(F,1,i),l,Q);
2090: r=lg(R);
2091: for (j=1;j<r;j++)
2092: {
2093: V[++nbfact]=R[j];
2094: ex[nbfact]=mael(F,2,i);
2095: }
2096: }
2097: setlg(V,nbfact+1);
2098: setlg(ex,nbfact+1);
2099: tetpil=avma; y=cgetg(3,t_VEC);
2100: y[1]=lcopy(V);
2101: y[2]=lcopy(ex);
2102: (void)sort_factor(y,cmp_pol);
2103: return gerepile(ltop,tetpil,y);
2104: }
2105: GEN Fp_factor_rel(GEN P, GEN l, GEN Q)
2106: {
1.2 ! noro 2107: gpmem_t tetpil, av=avma;
1.1 noro 2108: long nbfact;
2109: long j;
2110: GEN y,u,v;
2111: GEN z=Fp_factor_rel0(P,l,Q),t = (GEN) z[1],ex = (GEN) z[2];
2112: nbfact=lg(t);
2113: tetpil=avma; y=cgetg(3,t_MAT);
2114: u=cgetg(nbfact,t_COL); y[1]=(long)u;
2115: v=cgetg(nbfact,t_COL); y[2]=(long)v;
2116: for (j=1; j<nbfact; j++)
2117: {
2118: u[j] = lcopy((GEN)t[j]);
2119: v[j] = lstoi(ex[j]);
2120: }
2121: return gerepile(av,tetpil,y);
2122: }
2123:
2124: /*******************************************************************/
2125: extern int ff_poltype(GEN *x, GEN *p, GEN *pol);
2126:
2127: static GEN
2128: to_intmod(GEN x, GEN p)
2129: {
2130: GEN a = cgetg(3,t_INTMOD);
2131: a[1] = (long)p;
2132: a[2] = lmodii(x,p); return a;
2133: }
2134:
2135: /* z in Z[X], return z * Mod(1,p), normalized*/
2136: GEN
2137: FpX(GEN z, GEN p)
2138: {
2139: long i,l = lgef(z);
2140: GEN x = cgetg(l,t_POL);
2141: if (isonstack(p)) p = icopy(p);
2142: for (i=2; i<l; i++) x[i] = (long)to_intmod((GEN)z[i], p);
2143: x[1] = z[1]; return normalizepol_i(x,l);
2144: }
2145:
2146: /* z in Z^n, return z * Mod(1,p), normalized*/
2147: GEN
2148: FpV(GEN z, GEN p)
2149: {
2150: long i,l = lg(z);
2151: GEN x = cgetg(l,typ(z));
2152: if (isonstack(p)) p = icopy(p);
2153: for (i=1; i<l; i++) x[i] = (long)to_intmod((GEN)z[i], p);
2154: return x;
2155: }
2156: /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
2157: GEN
2158: FpM(GEN z, GEN p)
2159: {
2160: long i,j,l = lg(z), m = lg((GEN)z[1]);
2161: GEN x,y,zi;
2162: if (isonstack(p)) p = icopy(p);
2163: x = cgetg(l,t_MAT);
2164: for (i=1; i<l; i++)
2165: {
2166: x[i] = lgetg(m,t_COL);
2167: y = (GEN)x[i];
2168: zi= (GEN)z[i];
2169: for (j=1; j<m; j++) y[j] = (long)to_intmod((GEN)zi[j], p);
2170: }
2171: return x;
2172: }
2173: /* z in Z[X] or Z, return lift(z * Mod(1,p)), normalized*/
2174: GEN
2175: FpX_red(GEN z, GEN p)
2176: {
2177: long i,l;
2178: GEN x;
2179: if (typ(z) == t_INT) return modii(z,p);
2180: l = lgef(z);
2181: x = cgetg(l,t_POL);
2182: for (i=2; i<l; i++) x[i] = lmodii((GEN)z[i],p);
2183: x[1] = z[1]; return normalizepol_i(x,l);
2184: }
2185:
2186: GEN
2187: FpXV_red(GEN z, GEN p)
2188: {
2189: long i,l = lg(z);
2190: GEN x = cgetg(l, t_VEC);
2191: for (i=1; i<l; i++) x[i] = (long)FpX_red((GEN)z[i], p);
2192: return x;
2193: }
2194:
2195: /* z in Z^n, return lift(z * Mod(1,p)) */
2196: GEN
2197: FpV_red(GEN z, GEN p)
2198: {
2199: long i,l = lg(z);
2200: GEN x = cgetg(l,typ(z));
2201: for (i=1; i<l; i++) x[i] = lmodii((GEN)z[i],p);
2202: return x;
2203: }
2204:
2205: /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
2206: GEN
2207: FpM_red(GEN z, GEN p)
2208: {
1.2 ! noro 2209: long i, l = lg(z);
! 2210: GEN x = cgetg(l,t_MAT);
! 2211: for (i=1; i<l; i++) x[i] = (long)FpV_red((GEN)z[i], p);
1.1 noro 2212: return x;
2213: }
2214:
2215: /* no garbage collection, divide by leading coeff, mod p */
2216: GEN
2217: FpX_normalize(GEN z, GEN p)
2218: {
2219: GEN p1 = leading_term(z);
2220: if (gcmp1(p1)) return z;
2221: return FpX_Fp_mul(z, mpinvmod(p1,p), p);
2222: }
2223:
2224: GEN
2225: FpXQX_normalize(GEN z, GEN T, GEN p)
2226: {
2227: GEN p1 = leading_term(z);
2228: if (gcmp1(p1)) return z;
2229: if (!T) return FpX_normalize(z,p);
2230: return FpXQX_FpXQ_mul(z, FpXQ_inv(p1,T,p), T,p);
2231: }
2232:
2233: GEN
1.2 ! noro 2234: small_to_vec(GEN z)
! 2235: {
! 2236: long i, l = lg(z);
! 2237: GEN x = cgetg(l,t_VEC);
! 2238: for (i=1; i<l; i++) x[i] = lstoi(z[i]);
! 2239: return x;
! 2240: }
! 2241:
! 2242: GEN
! 2243: small_to_col(GEN z)
! 2244: {
! 2245: long i, l = lg(z);
! 2246: GEN x = cgetg(l,t_COL);
! 2247: for (i=1; i<l; i++) x[i] = lstoi(z[i]);
! 2248: return x;
! 2249: }
! 2250:
! 2251: GEN
1.1 noro 2252: small_to_mat(GEN z)
2253: {
2254: long i, l = lg(z);
2255: GEN x = cgetg(l,t_MAT);
1.2 ! noro 2256: for (i=1; i<l; i++)
! 2257: x[i] = (long)small_to_col((GEN)z[i]);
1.1 noro 2258: return x;
2259: }
2260:
2261: GEN
2262: small_to_pol_i(GEN z, long l)
2263: {
2264: long i;
2265: GEN x = cgetg(l,t_POL);
2266: for (i=2; i<l; i++) x[i] = lstoi(z[i]);
2267: x[1] = z[1]; return x;
2268: }
2269:
2270: /* assume z[i] % p has been done */
2271: GEN
2272: small_to_pol(GEN z, long v)
2273: {
2274: GEN x = small_to_pol_i(z, lgef(z));
2275: setvarn(x,v); return x;
2276: }
1.2 ! noro 2277: /* same. In place */
! 2278: GEN
! 2279: small_to_pol_ip(GEN z, long v)
! 2280: {
! 2281: long i, l = lgef(z);
! 2282: for (i=2; i<l; i++) z[i] = lstoi(z[i]);
! 2283: settyp(z, t_POL); setvarn(z, v); return z;
! 2284: }
1.1 noro 2285:
2286: GEN
2287: pol_to_small(GEN x)
2288: {
2289: long i, lx = lgef(x);
1.2 ! noro 2290: GEN a = u_getpol(lx-3);
1.1 noro 2291: for (i=2; i<lx; i++) a[i] = itos((GEN)x[i]);
2292: return a;
2293: }
2294:
1.2 ! noro 2295: /* z in Z[X,Y] representing an elt in F_p[X,Y] mod T(Y) i.e F_q[X])
1.1 noro 2296: * in Kronecker form. Recover the "real" z, normalized
2297: * If p = NULL, use generic functions and the coeff. ring implied by the
2298: * coefficients of z */
2299: GEN
1.2 ! noro 2300: FqX_from_Kronecker(GEN z, GEN T, GEN p)
1.1 noro 2301: {
1.2 ! noro 2302: long i,j,lx,l = lgef(z), N = (degpol(T)<<1) + 1;
1.1 noro 2303: GEN a,x, t = cgetg(N,t_POL);
1.2 ! noro 2304: t[1] = T[1] & VARNBITS;
1.1 noro 2305: lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
1.2 ! noro 2306: if (isonstack(T)) T = gcopy(T);
1.1 noro 2307: for (i=2; i<lx+2; i++)
2308: {
2309: a = cgetg(3,t_POLMOD); x[i] = (long)a;
1.2 ! noro 2310: a[1] = (long)T;
1.1 noro 2311: for (j=2; j<N; j++) t[j] = z[j];
2312: z += (N-2);
1.2 ! noro 2313: a[2] = (long)FpX_res(normalizepol_i(t,N), T,p);
1.1 noro 2314: }
2315: a = cgetg(3,t_POLMOD); x[i] = (long)a;
1.2 ! noro 2316: a[1] = (long)T;
1.1 noro 2317: N = (l-2) % (N-2) + 2;
2318: for (j=2; j<N; j++) t[j] = z[j];
1.2 ! noro 2319: a[2] = (long)FpX_res(normalizepol_i(t,N), T,p);
1.1 noro 2320: return normalizepol_i(x, i+1);
2321: }
2322:
1.2 ! noro 2323: #if 0
! 2324: /* z in Kronecker form representing a polynomial in FqX. Reduce mod (p,T) */
! 2325: GEN
! 2326: FqX_Kronecker_red(GEN z, GEN T, GEN p)
! 2327: {
! 2328: long i,j,lx,l = lgef(z), lT = lgef(T), N = ((dT-3)<<1) + 1;
! 2329: GEN a,x,y, t = cgetg(N,t_POL);
! 2330: t[1] = T[1] & VARNBITS;
! 2331: lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
! 2332: x = y = z = FpX_red(z, p);
! 2333: for (i=2; i<lx+2; i++)
! 2334: {
! 2335: for (j=2; j<N; j++) t[j] = z[j];
! 2336: a = FpX_res(normalizepol_i(t,N), T,p);
! 2337: for (j=2; j<lT; j++) y[j] = a[j];
! 2338: for ( ; j<N; j++) y[j] = zero;
! 2339: z += (N-2);
! 2340: y += (N-2);
! 2341: }
! 2342: N = (l-2) % (N-2) + 2;
! 2343: for (j=2; j<N; j++) t[j] = z[j];
! 2344: a = FpX_res(normalizepol_i(t,N), T,p);
! 2345: for (j=2; j<lT; j++) y[j] = a[j];
! 2346: for ( ; j<N; j++) y[j] = zero;
! 2347: return normalizepol_i(x, l);
! 2348: }
! 2349: #endif
! 2350:
1.1 noro 2351: /* z in ?[X,Y] mod Q(Y) in Kronecker form ((subst(lift(z), x, y^(2deg(z)-1)))
2352: * Recover the "real" z, normalized */
2353: GEN
2354: from_Kronecker(GEN z, GEN pol)
2355: {
2356: return FqX_from_Kronecker(z,pol,NULL);
2357: }
2358:
2359: /*******************************************************************/
2360: /* */
2361: /* MODULAR GCD */
2362: /* */
2363: /*******************************************************************/
1.2 ! noro 2364: extern ulong xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s);
1.1 noro 2365: /* 1 / Mod(x,p) , or 0 if inverse doesn't exist */
2366: ulong
2367: u_invmod(ulong x, ulong p)
2368: {
2369: long s;
2370: ulong xv, xv1, g = xgcduu(p, x, 1, &xv, &xv1, &s);
2371: if (g != 1UL) return 0UL;
2372: xv = xv1 % p; if (s < 0) xv = p - xv;
2373: return xv;
2374: }
2375:
1.2 ! noro 2376: int
! 2377: u_val(ulong n, ulong p)
! 2378: {
! 2379: ulong dummy;
! 2380: return svaluation(n,p,&dummy);
! 2381: }
! 2382:
! 2383: /* assume p^k is SMALL */
! 2384: int
! 2385: u_pow(int p, int k)
! 2386: {
! 2387: int i, pk;
! 2388:
! 2389: if (!k) return 1;
! 2390: if (p == 2) return 1<<k;
! 2391: pk = p; for (i=2; i<=k; i++) pk *= p;
! 2392: return pk;
! 2393: }
! 2394:
! 2395: #if 0
! 2396: static ulong
! 2397: umodratu(GEN a, ulong p)
1.1 noro 2398: {
2399: if (typ(a) == t_INT)
2400: return umodiu(a,p);
2401: else { /* assume a a t_FRAC */
2402: ulong num = umodiu((GEN)a[1],p);
2403: ulong den = umodiu((GEN)a[2],p);
2404: return (ulong)mulssmod(num, u_invmod(den,p), p);
2405: }
2406: }
2407: #endif
2408:
2409: /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL/VECSMALL/INT */
2410: GEN
1.2 ! noro 2411: u_Fp_FpX(GEN x, ulong p)
1.1 noro 2412: {
2413: long i, lx;
2414: GEN a;
2415:
2416: switch (typ(x))
2417: {
2418: case t_VECSMALL: return x;
1.2 ! noro 2419: case t_INT: a = u_getpol(0);
1.1 noro 2420: a[2] = umodiu(x, p); return a;
2421: }
1.2 ! noro 2422: lx = lgef(x); a = u_getpol(lx-3);
1.1 noro 2423: for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p);
2424: return u_normalizepol(a,lx);
2425: }
2426:
2427: GEN
1.2 ! noro 2428: u_Fp_FpV(GEN x, ulong p)
! 2429: {
! 2430: long i, n = lg(x);
! 2431: GEN y = cgetg(n,t_VECSMALL);
! 2432: for (i=1; i<n; i++) y[i] = (long)umodiu((GEN)x[i], p);
! 2433: return y;
! 2434: }
! 2435:
! 2436: GEN
1.1 noro 2437: u_Fp_FpM(GEN x, ulong p)
2438: {
1.2 ! noro 2439: long j,n = lg(x);
1.1 noro 2440: GEN y = cgetg(n,t_MAT);
2441: if (n == 1) return y;
1.2 ! noro 2442: for (j=1; j<n; j++) y[j] = (long)u_Fp_FpV((GEN)x[j], p);
1.1 noro 2443: return y;
2444: }
2445:
2446: static GEN
1.2 ! noro 2447: u_FpX_Fp_mul(GEN y, ulong x,ulong p)
1.1 noro 2448: {
2449: GEN z;
2450: int i, l;
1.2 ! noro 2451: if (!x) return u_zeropol();
! 2452: l = lgef(y); z = u_getpol(l-3);
1.1 noro 2453: if (HIGHWORD(x | p))
2454: for(i=2; i<l; i++) z[i] = mulssmod(y[i], x, p);
2455: else
2456: for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
2457: return z;
2458: }
2459:
2460: GEN
2461: u_FpX_normalize(GEN z, ulong p)
2462: {
2463: long l = lgef(z)-1;
2464: ulong p1 = z[l]; /* leading term */
2465: if (p1 == 1) return z;
1.2 ! noro 2466: return u_FpX_Fp_mul(z, u_invmod(p1,p), p);
1.1 noro 2467: }
2468:
2469: static GEN
1.2 ! noro 2470: u_copy(GEN x)
1.1 noro 2471: {
2472: long i, l = lgef(x);
1.2 ! noro 2473: GEN z = u_getpol(l-3);
1.1 noro 2474: for (i=2; i<l; i++) z[i] = x[i];
2475: return z;
2476: }
2477:
1.2 ! noro 2478: /* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES
! 2479: * if relevant, *pr is the last object on stack */
1.1 noro 2480: GEN
1.2 ! noro 2481: u_FpX_divrem(GEN x, GEN y, ulong p, GEN *pr)
1.1 noro 2482: {
2483: GEN z,q,c;
2484: long dx,dy,dz,i,j;
2485: ulong p1,inv;
2486:
1.2 ! noro 2487: if (pr == ONLY_REM) return u_FpX_rem(x, y, p);
1.1 noro 2488: dy = degpol(y);
2489: if (!dy)
2490: {
1.2 ! noro 2491: if (y[2] == 1UL)
! 2492: q = u_copy(x);
! 2493: else
! 2494: q = u_FpX_Fp_mul(x, u_invmod(y[2], p), p);
! 2495: if (pr) *pr = u_zeropol();
! 2496: return q;
1.1 noro 2497: }
2498: dx = degpol(x);
2499: dz = dx-dy;
2500: if (dz < 0)
2501: {
1.2 ! noro 2502: q = u_zeropol();
! 2503: if (pr) *pr = u_copy(x);
! 2504: return q;
1.1 noro 2505: }
2506: x += 2;
2507: y += 2;
1.2 ! noro 2508: z = u_getpol(dz) + 2;
1.1 noro 2509: inv = y[dy];
2510: if (inv != 1UL) inv = u_invmod(inv,p);
2511:
2512: if (u_OK_ULONG(p))
2513: {
2514: z[dz] = (inv*x[dx]) % p;
2515: for (i=dx-1; i>=dy; --i)
2516: {
2517: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
2518: for (j=i-dy+1; j<=i && j<=dz; j++)
2519: {
2520: p1 += z[j]*y[i-j];
1.2 ! noro 2521: if (p1 & HIGHBIT) p1 %= p;
1.1 noro 2522: }
2523: p1 %= p;
2524: z[i-dy] = p1? ((p - p1)*inv) % p: 0;
2525: }
2526: }
2527: else
2528: {
2529: z[dz] = mulssmod(inv, x[dx], p);
2530: for (i=dx-1; i>=dy; --i)
2531: {
2532: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
2533: for (j=i-dy+1; j<=i && j<=dz; j++)
2534: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
2535: z[i-dy] = p1? mulssmod(p - p1, inv, p): 0;
2536: }
2537: }
2538: q = u_normalizepol(z-2, dz+3);
2539: if (!pr) return q;
2540:
1.2 ! noro 2541: c = u_getpol(dy) + 2;
1.1 noro 2542: if (u_OK_ULONG(p))
2543: {
2544: for (i=0; i<dy; i++)
2545: {
2546: p1 = z[0]*y[i];
2547: for (j=1; j<=i && j<=dz; j++)
2548: {
2549: p1 += z[j]*y[i-j];
1.2 ! noro 2550: if (p1 & HIGHBIT) p1 %= p;
1.1 noro 2551: }
2552: c[i] = subssmod(x[i], p1%p, p);
2553: }
2554: }
2555: else
2556: {
2557: for (i=0; i<dy; i++)
2558: {
2559: p1 = mulssmod(z[0],y[i],p);
2560: for (j=1; j<=i && j<=dz; j++)
2561: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
2562: c[i] = subssmod(x[i], p1, p);
2563: }
2564: }
2565: i=dy-1; while (i>=0 && !c[i]) i--;
2566: c = u_normalizepol(c-2, i+3);
2567: *pr = c; return q;
2568: }
2569:
1.2 ! noro 2570: /*FIXME: Unify the following 3 divrem routines. Treat the case x,y (lifted) in
! 2571: * R[X], y non constant. Given: (lifted) [inv(), mul()], (delayed) red() in R */
! 2572:
1.1 noro 2573: /* x and y in Z[X]. Possibly x in Z */
2574: GEN
2575: FpX_divres(GEN x, GEN y, GEN p, GEN *pr)
2576: {
1.2 ! noro 2577: long vx, dx, dy, dz, i, j, sx, lrem;
! 2578: gpmem_t av0, av, tetpil;
1.1 noro 2579: GEN z,p1,rem,lead;
2580:
2581: if (!p) return poldivres(x,y,pr);
2582: if (!signe(y)) err(talker,"division by zero in FpX_divres");
1.2 ! noro 2583: vx = varn(x);
! 2584: dy = degpol(y);
! 2585: dx = (typ(x)==t_INT)? 0: degpol(x);
1.1 noro 2586: if (dx < dy)
2587: {
2588: if (pr)
2589: {
2590: av0 = avma; x = FpX_red(x, p);
2591: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
2592: if (pr == ONLY_REM) return x;
2593: *pr = x;
2594: }
2595: return zeropol(vx);
2596: }
2597: lead = leading_term(y);
2598: if (!dy) /* y is constant */
2599: {
2600: if (pr && pr != ONLY_DIVIDES)
2601: {
2602: if (pr == ONLY_REM) return zeropol(vx);
2603: *pr = zeropol(vx);
2604: }
2605: if (gcmp1(lead)) return gcopy(x);
2606: av0 = avma; x = gmul(x, mpinvmod(lead,p)); tetpil = avma;
2607: return gerepile(av0,tetpil,FpX_red(x,p));
2608: }
2609: av0 = avma; dz = dx-dy;
2610: if (OK_ULONG(p))
2611: { /* assume ab != 0 mod p */
2612: ulong pp = (ulong)p[2];
1.2 ! noro 2613: GEN a = u_Fp_FpX(x, pp);
! 2614: GEN b = u_Fp_FpX(y, pp);
! 2615: z = u_FpX_divrem(a,b,pp, pr);
! 2616: avma = av0; /* HACK: assume pr last on stack, then z */
! 2617: setlg(z, lgef(z)); z = dummycopy(z);
1.1 noro 2618: if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
2619: {
1.2 ! noro 2620: setlg(*pr, lgef(*pr)); *pr = dummycopy(*pr);
! 2621: *pr = small_to_pol_ip(*pr,vx);
1.1 noro 2622: }
1.2 ! noro 2623: return small_to_pol_ip(z,vx);
1.1 noro 2624: }
2625: lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p));
2626: avma = av0;
2627: z=cgetg(dz+3,t_POL);
2628: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
2629: x += 2; y += 2; z += 2;
2630:
2631: p1 = (GEN)x[dx]; av = avma;
2632: z[dz] = lead? lpileupto(av, modii(mulii(p1,lead), p)): licopy(p1);
2633: for (i=dx-1; i>=dy; i--)
2634: {
2635: av=avma; p1=(GEN)x[i];
2636: for (j=i-dy+1; j<=i && j<=dz; j++)
2637: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
2638: if (lead) p1 = mulii(p1,lead);
2639: tetpil=avma; z[i-dy] = lpile(av,tetpil,modii(p1, p));
2640: }
2641: if (!pr) { if (lead) gunclone(lead); return z-2; }
2642:
1.2 ! noro 2643: rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3);
1.1 noro 2644: for (sx=0; ; i--)
2645: {
2646: p1 = (GEN)x[i];
2647: for (j=0; j<=i && j<=dz; j++)
2648: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
2649: tetpil=avma; p1 = modii(p1,p); if (signe(p1)) { sx = 1; break; }
2650: if (!i) break;
2651: avma=av;
2652: }
2653: if (pr == ONLY_DIVIDES)
2654: {
2655: if (lead) gunclone(lead);
2656: if (sx) { avma=av0; return NULL; }
1.2 ! noro 2657: avma = (gpmem_t)rem; return z-2;
1.1 noro 2658: }
2659: lrem=i+3; rem -= lrem;
2660: rem[0]=evaltyp(t_POL) | evallg(lrem);
2661: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
1.2 ! noro 2662: p1 = gerepile((gpmem_t)rem,tetpil,p1);
1.1 noro 2663: rem += 2; rem[i]=(long)p1;
2664: for (i--; i>=0; i--)
2665: {
2666: av=avma; p1 = (GEN)x[i];
2667: for (j=0; j<=i && j<=dz; j++)
2668: p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));
2669: tetpil=avma; rem[i]=lpile(av,tetpil, modii(p1,p));
2670: }
2671: rem -= 2;
2672: if (lead) gunclone(lead);
1.2 ! noro 2673: if (!sx) (void)normalizepol_i(rem, lrem);
1.1 noro 2674: if (pr == ONLY_REM) return gerepileupto(av0,rem);
2675: *pr = rem; return z-2;
2676: }
2677:
2678: /* x and y in Z[Y][X]. Assume T irreducible mod p */
2679: GEN
2680: FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr)
2681: {
1.2 ! noro 2682: long vx, dx, dy, dz, i, j, sx, lrem;
! 2683: gpmem_t av0, av, tetpil;
1.1 noro 2684: GEN z,p1,rem,lead;
2685:
2686: if (!p) return poldivres(x,y,pr);
2687: if (!T) return FpX_divres(x,y,p,pr);
2688: if (!signe(y)) err(talker,"division by zero in FpX_divres");
2689: vx=varn(x); dy=degpol(y); dx=degpol(x);
2690: if (dx < dy)
2691: {
2692: if (pr)
2693: {
2694: av0 = avma; x = FpXQX_red(x, T, p);
2695: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
2696: if (pr == ONLY_REM) return x;
2697: *pr = x;
2698: }
2699: return zeropol(vx);
2700: }
2701: lead = leading_term(y);
2702: if (!dy) /* y is constant */
2703: {
2704: if (pr && pr != ONLY_DIVIDES)
2705: {
2706: if (pr == ONLY_REM) return zeropol(vx);
2707: *pr = zeropol(vx);
2708: }
2709: if (gcmp1(lead)) return gcopy(x);
2710: av0 = avma; x = gmul(x, FpXQ_inv(lead,T,p)); tetpil = avma;
2711: return gerepile(av0,tetpil,FpXQX_red(x,T,p));
2712: }
2713: av0 = avma; dz = dx-dy;
2714: #if 0 /* FIXME: to be done */
2715: if (OK_ULONG(p))
2716: { /* assume ab != 0 mod p */
2717: }
2718: #endif
2719: lead = gcmp1(lead)? NULL: gclone(FpXQ_inv(lead,T,p));
2720: avma = av0;
2721: z=cgetg(dz+3,t_POL);
2722: z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
2723: x += 2; y += 2; z += 2;
2724:
2725: p1 = (GEN)x[dx]; av = avma;
2726: z[dz] = lead? lpileupto(av, FpX_res(gmul(p1,lead), T, p)): lcopy(p1);
2727: for (i=dx-1; i>=dy; i--)
2728: {
2729: av=avma; p1=(GEN)x[i];
2730: for (j=i-dy+1; j<=i && j<=dz; j++)
2731: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
2732: if (lead) p1 = gmul(FpX_res(p1, T, p), lead);
2733: tetpil=avma; z[i-dy] = lpile(av,tetpil,FpX_res(p1, T, p));
2734: }
2735: if (!pr) { if (lead) gunclone(lead); return z-2; }
2736:
1.2 ! noro 2737: rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3);
1.1 noro 2738: for (sx=0; ; i--)
2739: {
2740: p1 = (GEN)x[i];
2741: for (j=0; j<=i && j<=dz; j++)
2742: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
2743: tetpil=avma; p1 = FpX_res(p1, T, p); if (signe(p1)) { sx = 1; break; }
2744: if (!i) break;
2745: avma=av;
2746: }
2747: if (pr == ONLY_DIVIDES)
2748: {
2749: if (lead) gunclone(lead);
2750: if (sx) { avma=av0; return NULL; }
1.2 ! noro 2751: avma = (gpmem_t)rem; return z-2;
1.1 noro 2752: }
2753: lrem=i+3; rem -= lrem;
2754: rem[0]=evaltyp(t_POL) | evallg(lrem);
2755: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
1.2 ! noro 2756: p1 = gerepile((gpmem_t)rem,tetpil,p1);
1.1 noro 2757: rem += 2; rem[i]=(long)p1;
2758: for (i--; i>=0; i--)
2759: {
2760: av=avma; p1 = (GEN)x[i];
2761: for (j=0; j<=i && j<=dz; j++)
2762: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
2763: tetpil=avma; rem[i]=lpile(av,tetpil, FpX_res(p1, T, p));
2764: }
2765: rem -= 2;
2766: if (lead) gunclone(lead);
1.2 ! noro 2767: if (!sx) (void)normalizepol_i(rem, lrem);
! 2768: if (pr == ONLY_REM) return gerepileupto(av0,rem);
! 2769: *pr = rem; return z-2;
! 2770: }
! 2771:
! 2772: /* R = any base ring */
! 2773: GEN
! 2774: RXQX_red(GEN P, GEN T)
! 2775: {
! 2776: long i, l = lgef(P);
! 2777: GEN Q = cgetg(l, t_POL);
! 2778: Q[1] = P[1];
! 2779: for (i=2; i<l; i++) Q[i] = lres((GEN)P[i], T);
! 2780: return Q;
! 2781: }
! 2782:
! 2783: /* R = any base ring */
! 2784: GEN
! 2785: RXQV_red(GEN P, GEN T)
! 2786: {
! 2787: long i, l = lg(P);
! 2788: GEN Q = cgetg(l, typ(P));
! 2789: for (i=1; i<l; i++) Q[i] = lres((GEN)P[i], T);
! 2790: return Q;
! 2791: }
! 2792:
! 2793: /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
! 2794: GEN
! 2795: RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
! 2796: {
! 2797: long vx, dx, dy, dz, i, j, sx, lrem;
! 2798: gpmem_t av0, av, tetpil;
! 2799: GEN z,p1,rem,lead;
! 2800:
! 2801: if (!signe(y)) err(talker,"division by zero in RXQX_divrem");
! 2802: vx = varn(x);
! 2803: dx = degpol(x);
! 2804: dy = degpol(y);
! 2805: if (dx < dy)
! 2806: {
! 2807: if (pr)
! 2808: {
! 2809: av0 = avma; x = RXQX_red(x, T);
! 2810: if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; }
! 2811: if (pr == ONLY_REM) return x;
! 2812: *pr = x;
! 2813: }
! 2814: return zeropol(vx);
! 2815: }
! 2816: lead = leading_term(y);
! 2817: if (!dy) /* y is constant */
! 2818: {
! 2819: if (pr && pr != ONLY_DIVIDES)
! 2820: {
! 2821: if (pr == ONLY_REM) return zeropol(vx);
! 2822: *pr = zeropol(vx);
! 2823: }
! 2824: if (gcmp1(lead)) return gcopy(x);
! 2825: av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
! 2826: return gerepile(av0,tetpil,RXQX_red(x,T));
! 2827: }
! 2828: av0 = avma; dz = dx-dy;
! 2829: lead = gcmp1(lead)? NULL: gclone(ginvmod(lead,T));
! 2830: avma = av0;
! 2831: z = cgetg(dz+3,t_POL);
! 2832: z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx);
! 2833: x += 2; y += 2; z += 2;
! 2834:
! 2835: p1 = (GEN)x[dx]; av = avma;
! 2836: z[dz] = lead? lpileupto(av, gres(gmul(p1,lead), T)): lcopy(p1);
! 2837: for (i=dx-1; i>=dy; i--)
! 2838: {
! 2839: av=avma; p1=(GEN)x[i];
! 2840: for (j=i-dy+1; j<=i && j<=dz; j++)
! 2841: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2842: if (lead) p1 = gmul(gres(p1, T), lead);
! 2843: tetpil=avma; z[i-dy] = lpile(av,tetpil, gres(p1, T));
! 2844: }
! 2845: if (!pr) { if (lead) gunclone(lead); return z-2; }
! 2846:
! 2847: rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3);
! 2848: for (sx=0; ; i--)
! 2849: {
! 2850: p1 = (GEN)x[i];
! 2851: for (j=0; j<=i && j<=dz; j++)
! 2852: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2853: tetpil=avma; p1 = gres(p1, T); if (signe(p1)) { sx = 1; break; }
! 2854: if (!i) break;
! 2855: avma=av;
! 2856: }
! 2857: if (pr == ONLY_DIVIDES)
! 2858: {
! 2859: if (lead) gunclone(lead);
! 2860: if (sx) { avma=av0; return NULL; }
! 2861: avma = (gpmem_t)rem; return z-2;
! 2862: }
! 2863: lrem=i+3; rem -= lrem;
! 2864: rem[0]=evaltyp(t_POL) | evallg(lrem);
! 2865: rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);
! 2866: p1 = gerepile((gpmem_t)rem,tetpil,p1);
! 2867: rem += 2; rem[i]=(long)p1;
! 2868: for (i--; i>=0; i--)
! 2869: {
! 2870: av=avma; p1 = (GEN)x[i];
! 2871: for (j=0; j<=i && j<=dz; j++)
! 2872: p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j]));
! 2873: tetpil=avma; rem[i]=lpile(av,tetpil, gres(p1, T));
! 2874: }
! 2875: rem -= 2;
! 2876: if (lead) gunclone(lead);
! 2877: if (!sx) (void)normalizepol_i(rem, lrem);
1.1 noro 2878: if (pr == ONLY_REM) return gerepileupto(av0,rem);
2879: *pr = rem; return z-2;
2880: }
2881:
2882: static GEN
2883: FpX_gcd_long(GEN x, GEN y, GEN p)
2884: {
1.2 ! noro 2885: ulong pp = (ulong)p[2];
! 2886: gpmem_t av = avma;
1.1 noro 2887: GEN a,b;
2888:
2889: (void)new_chunk((lgef(x) + lgef(y)) << 2); /* scratch space */
1.2 ! noro 2890: a = u_Fp_FpX(x, pp);
1.1 noro 2891: if (!signe(a)) { avma = av; return FpX_red(y,p); }
1.2 ! noro 2892: b = u_Fp_FpX(y, pp);
1.1 noro 2893: a = u_FpX_gcd_i(a,b, pp);
2894: avma = av; return small_to_pol(a, varn(x));
2895: }
2896:
2897: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)) */
2898: GEN
2899: FpX_gcd(GEN x, GEN y, GEN p)
2900: {
2901: GEN a,b,c;
1.2 ! noro 2902: gpmem_t av0, av;
1.1 noro 2903:
2904: if (OK_ULONG(p)) return FpX_gcd_long(x,y,p);
2905: av0=avma;
2906: a = FpX_red(x, p); av = avma;
2907: b = FpX_red(y, p);
2908: while (signe(b))
2909: {
2910: av = avma; c = FpX_res(a,b,p); a=b; b=c;
2911: }
2912: avma = av; return gerepileupto(av0, a);
2913: }
2914:
1.2 ! noro 2915: /*Return 1 if gcd can be computed
! 2916: * else return a factor of p*/
! 2917:
! 2918: GEN
! 2919: FpX_gcd_check(GEN x, GEN y, GEN p)
! 2920: {
! 2921: GEN a,b,c;
! 2922: gpmem_t av=avma;
! 2923:
! 2924: a = FpX_red(x, p);
! 2925: b = FpX_red(y, p);
! 2926: while (signe(b))
! 2927: {
! 2928: GEN lead = leading_term(b);
! 2929: GEN g = mppgcd(lead,p);
! 2930: if (!is_pm1(g)) return gerepileupto(av,g);
! 2931: c = FpX_res(a,b,p); a=b; b=c;
! 2932: }
! 2933: avma = av; return gun;
! 2934: }
! 2935:
1.1 noro 2936: GEN
2937: u_FpX_sub(GEN x, GEN y, ulong p)
2938: {
2939: long i,lz,lx = lgef(x), ly = lgef(y);
2940: GEN z;
2941:
2942: if (ly <= lx)
2943: {
2944: lz = lx; z = cgetg(lz,t_VECSMALL);
2945: for (i=2; i<ly; i++) z[i] = subssmod(x[i],y[i],p);
2946: for ( ; i<lx; i++) z[i] = x[i];
2947: }
2948: else
2949: {
2950: lz = ly; z = cgetg(lz,t_VECSMALL);
2951: for (i=2; i<lx; i++) z[i] = subssmod(x[i],y[i],p);
2952: for ( ; i<ly; i++) z[i] = y[i]? p - y[i]: y[i];
2953: }
2954: z[1]=0; return u_normalizepol(z, lz);
2955: }
2956:
2957: /* list of u_FpX in gptr, return then as GEN */
2958: static void
1.2 ! noro 2959: u_gerepilemany(gpmem_t av, GEN* gptr[], long n, long v)
1.1 noro 2960: {
2961: GEN *l = (GEN*)gpmalloc(n*sizeof(GEN));
2962: long i;
2963:
2964: /* copy objects */
2965: for (i=0; i<n; i++) l[i] = gclone(*(gptr[i]));
2966: avma = av;
2967: /* copy them back, kill clones */
2968: for (--i; i>=0; i--)
2969: {
2970: *(gptr[i]) = small_to_pol(l[i],v);
2971: gunclone(l[i]);
2972: }
2973: free(l);
2974: }
2975:
2976: static GEN
2977: u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
2978: {
2979: GEN q,z,u,v, x = a, y = b;
2980:
1.2 ! noro 2981: u = u_zeropol();
! 2982: v= u_scalarpol(1); /* v = 1 */
1.1 noro 2983: while (signe(y))
2984: {
1.2 ! noro 2985: q = u_FpX_divrem(x,y,p,&z);
1.1 noro 2986: x = y; y = z; /* (x,y) = (y, x - q y) */
2987: z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
2988: u = v; v = z; /* (u,v) = (v, u - q v) */
2989: }
2990: z = u_FpX_sub(x, u_FpX_mul(b,u,p), p);
2991: z = u_FpX_div(z,a,p);
2992: *ptu = z;
2993: *ptv = u; return x;
2994: }
2995:
2996: static GEN
2997: FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
2998: {
2999: ulong pp = (ulong)p[2];
1.2 ! noro 3000: gpmem_t av = avma;
1.1 noro 3001: GEN a, b, d;
3002:
1.2 ! noro 3003: a = u_Fp_FpX(x, pp);
! 3004: b = u_Fp_FpX(y, pp);
1.1 noro 3005: d = u_FpX_extgcd(a,b, pp, ptu,ptv);
3006: {
3007: GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv;
3008: u_gerepilemany(av, gptr, 3, varn(x));
3009: }
3010: return d;
3011: }
3012:
3013: /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
3014: * ux + vy = gcd (mod p) */
1.2 ! noro 3015: /*TODO: Document the typ() of *ptu and *ptv*/
1.1 noro 3016: GEN
3017: FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
3018: {
3019: GEN a,b,q,r,u,v,d,d1,v1;
1.2 ! noro 3020: gpmem_t ltop, lbot;
1.1 noro 3021:
3022: if (OK_ULONG(p)) return FpX_extgcd_long(x,y,p,ptu,ptv);
3023: ltop=avma;
3024: a = FpX_red(x, p);
3025: b = FpX_red(y, p);
3026: d = a; d1 = b; v = gzero; v1 = gun;
3027: while (signe(d1))
3028: {
3029: q = FpX_divres(d,d1,p, &r);
3030: v = gadd(v, gneg_i(gmul(q,v1)));
3031: v = FpX_red(v,p);
3032: u=v; v=v1; v1=u;
3033: u=r; d=d1; d1=u;
3034: }
3035: u = gadd(d, gneg_i(gmul(b,v)));
3036: u = FpX_red(u, p);
3037: lbot = avma;
3038: u = FpX_div(u,a,p);
3039: d = gcopy(d);
3040: v = gcopy(v);
3041: {
3042: GEN *gptr[3]; gptr[0] = &d; gptr[1] = &u; gptr[2] = &v;
3043: gerepilemanysp(ltop,lbot,gptr,3);
3044: }
3045: *ptu = u; *ptv = v; return d;
3046: }
3047:
3048: /* x and y in Z[Y][X], return lift(gcd(x mod T,p, y mod T,p)). Set u and v st
3049: * ux + vy = gcd (mod T,p) */
3050: GEN
3051: FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
3052: {
3053: GEN a,b,q,r,u,v,d,d1,v1;
1.2 ! noro 3054: gpmem_t ltop, lbot;
1.1 noro 3055:
3056: #if 0 /* FIXME To be done...*/
3057: if (OK_ULONG(p)) return FpXQX_extgcd_long(x,y,T,p,ptu,ptv);
3058: #endif
3059: if (!T) return FpX_extgcd(x,y,p,ptu,ptv);
3060: ltop=avma;
3061: a = FpXQX_red(x, T, p);
3062: b = FpXQX_red(y, T, p);
3063: d = a; d1 = b; v = gzero; v1 = gun;
3064: while (signe(d1))
3065: {
3066: q = FpXQX_divres(d,d1,T,p, &r);
3067: v = gadd(v, gneg_i(gmul(q,v1)));
3068: v = FpXQX_red(v,T,p);
3069: u=v; v=v1; v1=u;
3070: u=r; d=d1; d1=u;
3071: }
3072: u = gadd(d, gneg_i(gmul(b,v)));
3073: u = FpXQX_red(u,T, p);
3074: lbot = avma;
3075: u = FpXQX_divres(u,a,T,p,NULL);
3076: d = gcopy(d);
3077: v = gcopy(v);
3078: {
3079: GEN *gptr[3]; gptr[0] = &d; gptr[1] = &u; gptr[2] = &v;
3080: gerepilemanysp(ltop,lbot,gptr,3);
3081: }
3082: *ptu = u; *ptv = v; return d;
3083: }
3084:
1.2 ! noro 3085: /*x must be reduced*/
1.1 noro 3086: GEN
3087: FpXQ_charpoly(GEN x, GEN T, GEN p)
3088: {
1.2 ! noro 3089: gpmem_t ltop=avma;
! 3090: long v=varn(T);
! 3091: GEN R;
! 3092: T=gcopy(T); setvarn(T,MAXVARN);
! 3093: x=gcopy(x); setvarn(x,MAXVARN);
! 3094: R=FpY_FpXY_resultant(T,deg1pol(gun,FpX_neg(x,p),v),p);
1.1 noro 3095: return gerepileupto(ltop,R);
3096: }
3097:
3098: GEN
3099: FpXQ_minpoly(GEN x, GEN T, GEN p)
3100: {
1.2 ! noro 3101: gpmem_t ltop=avma;
1.1 noro 3102: GEN R=FpXQ_charpoly(x, T, p);
3103: GEN G=FpX_gcd(R,derivpol(R),p);
3104: G=FpX_normalize(G,p);
3105: G=FpX_div(R,G,p);
3106: return gerepileupto(ltop,G);
3107: }
3108:
3109: /* return z = a mod q, b mod p (p,q) = 1. qinv = 1/q mod p */
3110: static GEN
3111: u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq)
3112: {
1.2 ! noro 3113: ulong d, amod = umodiu(a, p);
! 3114: gpmem_t av = avma;
1.1 noro 3115: GEN ax;
3116:
3117: if (b == amod) return NULL;
3118: d = (b > amod)? b - amod: p - (amod - b); /* (b - a) mod p */
3119: (void)new_chunk(lgefint(pq)<<1); /* HACK */
3120: ax = mului(mulssmod(d,qinv,p), q); /* d mod p, 0 mod q */
3121: avma = av; return addii(a, ax); /* in ]-q, pq[ assuming a in -]-q,q[ */
3122: }
3123:
3124: /* centerlift(u mod p) */
1.2 ! noro 3125: long
1.1 noro 3126: u_center(ulong u, ulong p, ulong ps2)
3127: {
1.2 ! noro 3128: return (long) (u > ps2)? u - p: u;
1.1 noro 3129: }
3130:
3131: GEN
3132: ZX_init_CRT(GEN Hp, ulong p, long v)
3133: {
3134: long i, l = lgef(Hp), lim = (long)(p>>1);
3135: GEN H = cgetg(l, t_POL);
3136: H[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
3137: for (i=2; i<l; i++)
1.2 ! noro 3138: H[i] = lstoi(u_center(Hp[i], p, lim));
1.1 noro 3139: return H;
3140: }
3141:
3142: /* assume lg(Hp) > 1 */
3143: GEN
3144: ZM_init_CRT(GEN Hp, ulong p)
3145: {
3146: long i,j, m = lg(Hp[1]), l = lg(Hp), lim = (long)(p>>1);
3147: GEN c,cp,H = cgetg(l, t_MAT);
3148: for (j=1; j<l; j++)
3149: {
3150: cp = (GEN)Hp[j];
3151: c = cgetg(m, t_COL);
3152: H[j] = (long)c;
1.2 ! noro 3153: for (i=1; i<l; i++) c[i] = lstoi(u_center(cp[i],p, lim));
1.1 noro 3154: }
3155: return H;
3156: }
3157:
3158: int
3159: Z_incremental_CRT(GEN *H, ulong Hp, GEN q, GEN qp, ulong p)
3160: {
3161: GEN h, lim = shifti(qp,-1);
3162: ulong qinv = u_invmod(umodiu(q,p), p);
3163: int stable = 1;
3164: h = u_chrem_coprime(*H,Hp,q,p,qinv,qp);
3165: if (h)
3166: {
3167: if (cmpii(h,lim) > 0) h = subii(h,qp);
3168: *H = h; stable = 0;
3169: }
3170: return stable;
3171: }
3172:
3173: int
1.2 ! noro 3174: ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
1.1 noro 3175: {
1.2 ! noro 3176: GEN H = *ptH, h, lim = shifti(qp,-1);
1.1 noro 3177: ulong qinv = u_invmod(umodiu(q,p), p);
3178: long i, l = lgef(H), lp = lgef(Hp);
3179: int stable = 1;
1.2 ! noro 3180:
! 3181: if (l < lp)
! 3182: { /* degree increases */
! 3183: GEN x = cgetg(lp, t_POL);
! 3184: for (i=1; i<l; i++) x[i] = H[i];
! 3185: for ( ; i<lp; i++)
! 3186: {
! 3187: h = stoi(Hp[i]);
! 3188: if (cmpii(h,lim) > 0) h = subii(h,qp);
! 3189: x[i] = (long)h;
! 3190: }
! 3191: setlgef(x,lp); *ptH = H = x;
! 3192: stable = 0; lp = l;
! 3193: }
1.1 noro 3194: for (i=2; i<lp; i++)
3195: {
3196: h = u_chrem_coprime((GEN)H[i],Hp[i],q,p,qinv,qp);
3197: if (h)
3198: {
3199: if (cmpii(h,lim) > 0) h = subii(h,qp);
3200: H[i] = (long)h; stable = 0;
3201: }
3202: }
3203: for ( ; i<l; i++)
3204: {
3205: h = u_chrem_coprime((GEN)H[i], 0,q,p,qinv,qp);
3206: if (h)
3207: {
3208: if (cmpii(h,lim) > 0) h = subii(h,qp);
3209: H[i] = (long)h; stable = 0;
3210: }
3211: }
3212: return stable;
3213: }
3214:
3215: int
3216: ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p)
3217: {
3218: GEN h, lim = shifti(qp,-1);
3219: ulong qinv = u_invmod(umodiu(q,p), p);
3220: long i,j, l = lg(H), m = lg(H[1]);
3221: int stable = 1;
3222: for (j=1; j<l; j++)
3223: for (i=1; i<m; i++)
3224: {
3225: h = u_chrem_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp);
3226: if (h)
3227: {
3228: if (cmpii(h,lim) > 0) h = subii(h,qp);
3229: coeff(H,i,j) = (long)h; stable = 0;
3230: }
3231: }
3232: return stable;
3233: }
3234:
3235: /* returns a polynomial in variable v, whose coeffs correspond to the digits
3236: * of m (in base p) */
3237: GEN
3238: stopoly(long m, long p, long v)
3239: {
3240: GEN y = cgetg(BITS_IN_LONG + 2, t_POL);
3241: long l=2;
3242:
3243: do { y[l++] = lstoi(m%p); m=m/p; } while (m);
3244: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
3245: return y;
3246: }
3247:
3248: GEN
3249: stopoly_gen(GEN m, GEN p, long v)
3250: {
3251: GEN y = cgetg(bit_accuracy(lgefint(m))+2, t_POL);
3252: long l=2;
3253:
3254: do { y[l++] = lmodii(m,p); m=divii(m,p); } while (signe(m));
3255: y[1] = evalsigne(1)|evallgef(l)|evalvarn(v);
3256: return y;
3257: }
3258:
3259: /* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0 */
3260: GEN
3261: u_FpX_rem(GEN x, GEN y, ulong p)
3262: {
3263: GEN z, c;
3264: long dx,dy,dz,i,j;
3265: ulong p1,inv;
3266:
1.2 ! noro 3267: dy = degpol(y); if (!dy) return u_zeropol();
1.1 noro 3268: dx = degpol(x);
1.2 ! noro 3269: dz = dx-dy; if (dz < 0) return u_copy(x);
1.1 noro 3270: x += 2;
3271: y += 2;
1.2 ! noro 3272: z = u_mallocpol(dz) + 2;
1.1 noro 3273: inv = y[dy];
3274: if (inv != 1UL) inv = u_invmod(inv,p);
3275:
1.2 ! noro 3276: c = u_getpol(dy) + 2;
1.1 noro 3277: if (u_OK_ULONG(p))
3278: {
3279: z[dz] = (inv*x[dx]) % p;
3280: for (i=dx-1; i>=dy; --i)
3281: {
3282: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
3283: for (j=i-dy+1; j<=i && j<=dz; j++)
3284: {
3285: p1 += z[j]*y[i-j];
1.2 ! noro 3286: if (p1 & HIGHBIT) p1 %= p;
1.1 noro 3287: }
3288: p1 %= p;
3289: z[i-dy] = p1? ((p - p1)*inv) % p: 0;
3290: }
3291: for (i=0; i<dy; i++)
3292: {
3293: p1 = z[0]*y[i];
3294: for (j=1; j<=i && j<=dz; j++)
3295: {
3296: p1 += z[j]*y[i-j];
1.2 ! noro 3297: if (p1 & HIGHBIT) p1 %= p;
1.1 noro 3298: }
3299: c[i] = subssmod(x[i], p1%p, p);
3300: }
3301: }
3302: else
3303: {
3304: z[dz] = mulssmod(inv, x[dx], p);
3305: for (i=dx-1; i>=dy; --i)
3306: {
3307: p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
3308: for (j=i-dy+1; j<=i && j<=dz; j++)
3309: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
3310: z[i-dy] = p1? mulssmod(p - p1, inv, p): 0;
3311: }
3312: for (i=0; i<dy; i++)
3313: {
3314: p1 = mulssmod(z[0],y[i],p);
3315: for (j=1; j<=i && j<=dz; j++)
3316: p1 = addssmod(p1, mulssmod(z[j],y[i-j],p), p);
3317: c[i] = subssmod(x[i], p1, p);
3318: }
3319: }
3320: i = dy-1; while (i>=0 && !c[i]) i--;
3321: free(z-2); return u_normalizepol(c-2, i+3);
3322: }
3323:
3324: ulong
3325: u_FpX_resultant(GEN a, GEN b, ulong p)
3326: {
3327: long da,db,dc,cnt;
1.2 ! noro 3328: ulong lb, res = 1UL;
! 3329: gpmem_t av;
1.1 noro 3330: GEN c;
3331:
3332: if (!signe(a) || !signe(b)) return 0;
3333: da = degpol(a);
3334: db = degpol(b);
3335: if (db > da)
3336: {
3337: swapspec(a,b, da,db);
1.2 ! noro 3338: if (both_odd(da,db)) res = p-res;
1.1 noro 3339: }
3340: if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
3341: cnt = 0; av = avma;
3342: while (db)
3343: {
3344: lb = b[db+2];
3345: c = u_FpX_rem(a,b, p);
3346: a = b; b = c; dc = degpol(c);
3347: if (dc < 0) { avma = av; return 0; }
3348:
1.2 ! noro 3349: if (both_odd(da,db)) res = p - res;
1.1 noro 3350: if (lb != 1) res = mulssmod(res, powuumod(lb, da - dc, p), p);
3351: if (++cnt == 4) { cnt = 0; avma = av; }
3352: da = db; /* = degpol(a) */
3353: db = dc; /* = degpol(b) */
3354: }
3355: avma = av; return mulssmod(res, powuumod(b[2], da, p), p);
3356: }
3357:
1.2 ! noro 3358: static GEN
! 3359: muliimod(GEN x, GEN y, GEN p)
! 3360: {
! 3361: return modii(mulii(x,y), p);
! 3362: }
! 3363:
! 3364: #define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM)
! 3365:
! 3366: /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab+a+r), with R=A%B, a=deg(A) ...*/
! 3367: GEN
! 3368: FpX_resultant(GEN a, GEN b, GEN p)
! 3369: {
! 3370: long da,db,dc,cnt;
! 3371: gpmem_t av, lim;
! 3372: GEN c,lb, res = gun;
! 3373:
! 3374: if (!signe(a) || !signe(b)) return gzero;
! 3375: da = degpol(a);
! 3376: db = degpol(b);
! 3377: if (db > da)
! 3378: {
! 3379: swapspec(a,b, da,db);
! 3380: if (both_odd(da,db)) res = subii(p, res);
! 3381: }
! 3382: if (!da) return gun; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
! 3383: cnt = 0; av = avma; lim = stack_lim(av,2);
! 3384: while (db)
! 3385: {
! 3386: lb = (GEN)b[db+2];
! 3387: c = FpX_rem(a,b, p);
! 3388: a = b; b = c; dc = degpol(c);
! 3389: if (dc < 0) { avma = av; return 0; }
! 3390:
! 3391: if (both_odd(da,db)) res = subii(p, res);
! 3392: if (!gcmp1(lb)) res = muliimod(res, powgumod(lb, da - dc, p), p);
! 3393: if (low_stack(lim,stack_lim(av,2)))
! 3394: {
! 3395: if (DEBUGMEM>1) err(warnmem,"FpX_resultant (da = %ld)",da);
! 3396: gerepileall(av,3, &a,&b,&res);
! 3397: }
! 3398: da = db; /* = degpol(a) */
! 3399: db = dc; /* = degpol(b) */
! 3400: }
! 3401: res = muliimod(res, powgumod((GEN)b[2], da, p), p);
! 3402: return gerepileuptoint(av, res);
! 3403: }
! 3404:
1.1 noro 3405: /* If resultant is 0, *ptU and *ptU are not set */
3406: ulong
3407: u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
3408: {
3409: GEN z,q,u,v, x = a, y = b;
1.2 ! noro 3410: ulong lb, res = 1UL;
! 3411: gpmem_t av = avma;
1.1 noro 3412: long dx,dy,dz;
3413:
3414: if (!signe(x) || !signe(y)) return 0;
3415: dx = degpol(x);
3416: dy = degpol(y);
3417: if (dy > dx)
3418: {
3419: swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
3420: a = x; b = y;
1.2 ! noro 3421: if (both_odd(dx,dy)) res = p-res;
1.1 noro 3422: }
1.2 ! noro 3423: u = u_zeropol();
! 3424: v = u_scalarpol(1); /* v = 1 */
1.1 noro 3425: while (dy)
3426: { /* b u = x (a), b v = y (a) */
3427: lb = y[dy+2];
1.2 ! noro 3428: q = u_FpX_divrem(x,y, p, &z);
1.1 noro 3429: x = y; y = z; /* (x,y) = (y, x - q y) */
3430: dz = degpol(z); if (dz < 0) { avma = av; return 0; }
3431: z = u_FpX_sub(u, u_FpX_mul(q,v, p), p);
3432: u = v; v = z; /* (u,v) = (v, u - q v) */
3433:
1.2 ! noro 3434: if (both_odd(dx,dy)) res = p - res;
1.1 noro 3435: if (lb != 1) res = mulssmod(res, powuumod(lb, dx-dz, p), p);
3436: dx = dy; /* = degpol(x) */
3437: dy = dz; /* = degpol(y) */
3438: }
3439: res = mulssmod(res, powuumod(y[2], dx, p), p);
3440: lb = mulssmod(res, u_invmod(y[2],p), p);
1.2 ! noro 3441: v = gerepileupto(av, u_FpX_Fp_mul(v, lb, p));
1.1 noro 3442: av = avma;
1.2 ! noro 3443: u = u_FpX_sub(u_scalarpol(res), u_FpX_mul(b,v,p), p);
1.1 noro 3444: u = gerepileupto(av, u_FpX_div(u,a,p)); /* = (res - b v) / a */
3445: *ptU = u;
3446: *ptV = v; return res;
3447: }
3448:
3449: /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with "generic"
3450: * degree sequence as given by dglist, set *Ci and return resultant(a,b) */
3451: ulong
3452: u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
3453: {
3454: long da,db,dc,cnt,ind;
1.2 ! noro 3455: ulong lb, cx = 1, res = 1UL;
! 3456: gpmem_t av;
1.1 noro 3457: GEN c;
3458:
3459: if (C0) { *C0 = 1; *C1 = 0; }
3460: if (!signe(a) || !signe(b)) return 0;
3461: da = degpol(a);
3462: db = degpol(b);
3463: if (db > da)
3464: {
3465: swapspec(a,b, da,db);
1.2 ! noro 3466: if (both_odd(da,db)) res = p-res;
1.1 noro 3467: }
3468: /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
3469: if (!da) return 1;
3470: cnt = ind = 0; av = avma;
3471: while (db)
3472: {
3473: lb = b[db+2];
3474: c = u_FpX_rem(a,b, p);
3475: a = b; b = c; dc = degpol(c);
3476: if (dc < 0) { avma = av; return 0; }
3477:
3478: ind++;
3479: if (C0)
3480: { /* check that Euclidean remainder sequence doesn't degenerate */
3481: if (dc != dglist[ind]) { avma = av; return 0; }
3482: /* update resultant */
1.2 ! noro 3483: if (both_odd(da,db)) res = p-res;
1.1 noro 3484: if (lb != 1)
3485: {
3486: ulong t = powuumod(lb, da - dc, p);
3487: res = mulssmod(res, t, p);
3488: if (dc) cx = mulssmod(cx, t, p);
3489: }
3490: }
3491: else
3492: {
3493: if (dc > dglist[ind]) dglist[ind] = dc;
3494: }
3495: if (++cnt == 4) { cnt = 0; avma = av; }
3496: da = db; /* = degpol(a) */
3497: db = dc; /* = degpol(b) */
3498: }
3499: if (!C0)
3500: {
3501: if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
3502: return 0;
3503: }
3504:
3505: if (da == 1) /* last non-constant polynomial has degree 1 */
3506: {
3507: *C0 = mulssmod(cx, a[2], p);
3508: *C1 = mulssmod(cx, a[3], p);
3509: lb = b[2];
3510: } else lb = powuumod(b[2], da, p);
3511: avma = av; return mulssmod(res, lb, p);
3512: }
3513:
3514: static ulong global_pp;
3515: static GEN
3516: _u_FpX_mul(GEN a, GEN b)
3517: {
3518: return u_FpX_mul(a,b, global_pp);
3519: }
3520:
3521: /* compute prod (x - a[i]) */
3522: GEN
3523: u_FpV_roots_to_pol(GEN a, ulong p)
3524: {
3525: long i,k,lx = lg(a);
3526: GEN p1,p2;
3527: if (lx == 1) return polun[0];
3528: p1 = cgetg(lx, t_VEC); global_pp = p;
3529: for (k=1,i=1; i<lx-1; i+=2)
3530: {
3531: p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2;
3532: p2[2] = mulssmod(a[i], a[i+1], p);
1.2 ! noro 3533: p2[3] = a[i] + a[i+1];
! 3534: if (p2[3] >= p) p2[3] -= p;
! 3535: if (p2[3]) p2[3] = p - p2[3]; /* - (a[i] + a[i+1]) mod p */
1.1 noro 3536: p2[4] = 1; p2[1] = evallgef(5);
3537: }
3538: if (i < lx)
3539: {
3540: p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
3541: p2[1] = evallgef(4);
3542: p2[2] = p - a[i];
3543: p2[3] = 1;
3544: }
3545: setlg(p1, k); return divide_conquer_prod(p1, _u_FpX_mul);
3546: }
3547:
3548:
3549: /* u P(X) + v P(-X) */
3550: static GEN
3551: pol_comp(GEN P, GEN u, GEN v)
3552: {
3553: long i, l = lgef(P);
3554: GEN y = cgetg(l, t_POL);
3555: for (i=2; i<l; i++)
3556: {
3557: GEN t = (GEN)P[i];
3558: y[i] = gcmp0(t)? zero:
3559: (i&1)? lmul(t, gsub(u,v)) /* odd degree */
3560: : lmul(t, gadd(u,v));/* even degree */
3561: }
3562: y[1] = P[1]; return normalizepol_i(y,l);
3563: }
3564:
3565: static GEN
3566: u_pol_comp(GEN P, ulong u, ulong v, ulong p)
3567: {
3568: long i, l = lgef(P);
1.2 ! noro 3569: GEN y = u_getpol(l-3);
1.1 noro 3570: for (i=2; i<l; i++)
3571: {
3572: ulong t = P[i];
3573: y[i] = (t == 0)? 0:
3574: (i&1)? mulssmod(t, u + (p - v), p)
3575: : mulssmod(t, u + v, p);
3576: }
3577: return u_normalizepol(y,l);
3578: }
3579:
1.2 ! noro 3580: extern GEN roots_to_pol(GEN a, long v);
1.1 noro 3581:
3582: GEN
3583: polint_triv(GEN xa, GEN ya)
3584: {
3585: GEN P = NULL, Q = roots_to_pol(xa,0);
1.2 ! noro 3586: long i, n = lg(xa);
! 3587: gpmem_t av = avma, lim = stack_lim(av, 2);
1.1 noro 3588: for (i=1; i<n; i++)
3589: {
3590: GEN T,dP;
3591: if (gcmp0((GEN)ya[i])) continue;
3592: T = gdeuc(Q, gsub(polx[0], (GEN)xa[i]));
3593: if (i < n-1 && absi_equal((GEN)xa[i], (GEN)xa[i+1]))
3594: { /* x_i = -x_{i+1} */
3595: T = gdiv(T, poleval(T, (GEN)xa[i]));
3596: dP = pol_comp(T, (GEN)ya[i], (GEN)ya[i+1]);
3597: i++;
3598: }
3599: else
3600: {
3601: dP = gmul((GEN)ya[i], T);
3602: dP = gdiv(dP, poleval(T,(GEN)xa[i]));
3603: }
3604: P = P? gadd(P, dP): dP;
3605: if (low_stack(lim,stack_lim(av,2)))
3606: {
3607: if (DEBUGMEM>1) err(warnmem,"polint_triv2 (i = %ld)",i);
3608: P = gerepileupto(av, P);
3609: }
3610: }
3611: return P? P: zeropol(0);
3612: }
3613:
3614: ulong
3615: u_FpX_eval(GEN x, ulong y, ulong p)
3616: {
3617: ulong p1,r;
3618: long i,j;
3619: i=lgef(x)-1;
3620: if (i<=2)
3621: return (i==2)? x[2]: 0;
3622: p1 = x[i];
3623: /* specific attention to sparse polynomials (see poleval)*/
3624: if (u_OK_ULONG(p))
3625: {
3626: for (i--; i>=2; i=j-1)
3627: {
3628: for (j=i; !x[j]; j--)
3629: if (j==2)
3630: {
3631: if (i != j) y = powuumod(y, i-j+1, p);
3632: return (p1 * y) % p;
3633: }
3634: r = (i==j)? y: powuumod(y, i-j+1, p);
3635: p1 = ((p1*r) + x[j]) % p;
3636: }
3637: }
3638: else
3639: {
3640: for (i--; i>=2; i=j-1)
3641: {
3642: for (j=i; !x[j]; j--)
3643: if (j==2)
3644: {
3645: if (i != j) y = powuumod(y, i-j+1, p);
3646: return mulssmod(p1, y, p);
3647: }
3648: r = (i==j)? y: powuumod(y, i-j+1, p);
3649: p1 = addssmod(x[j], mulssmod(p1,r,p), p);
3650: }
3651: }
3652: return p1;
3653: }
3654:
3655: static GEN
3656: u_FpX_div_by_X_x(GEN a, ulong x, ulong p)
3657: {
3658: long l = lgef(a), i;
1.2 ! noro 3659: GEN z = u_getpol(l-4), a0, z0;
1.1 noro 3660: a0 = a + l-1;
3661: z0 = z + l-2; *z0 = *a0--;
3662: if (u_OK_ULONG(p))
3663: {
3664: for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
3665: {
3666: ulong t = (*a0-- + x * *z0--) % p;
3667: *z0 = t;
3668: }
3669: }
3670: else
3671: {
3672: for (i=l-3; i>1; i--)
3673: {
3674: ulong t = addssmod(*a0--, mulssmod(x, *z0--, p), p);
3675: *z0 = t;
3676: }
3677: }
3678: return z;
3679: }
3680:
1.2 ! noro 3681: static GEN
! 3682: FpX_div_by_X_x(GEN a, GEN x, GEN p)
! 3683: {
! 3684: long l = lgef(a), i;
! 3685: GEN z = cgetg(l-1, t_POL), a0, z0;
! 3686: z[1] = evalsigne(1)|evalvarn(0)|evallgef(l-1);
! 3687: a0 = a + l-1;
! 3688: z0 = z + l-2; *z0 = *a0--;
! 3689: for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
! 3690: {
! 3691: GEN t = addii((GEN)*a0--, muliimod(x, (GEN)*z0--, p));
! 3692: *z0 = (long)t;
! 3693: }
! 3694: return z;
! 3695: }
! 3696:
1.1 noro 3697: /* xa, ya = t_VECSMALL */
3698: GEN
3699: u_FpV_polint(GEN xa, GEN ya, ulong p)
3700: {
3701: GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p);
3702: long i, n = lg(xa);
1.2 ! noro 3703: ulong inv;
! 3704: gpmem_t av = avma;
1.1 noro 3705: for (i=1; i<n; i++)
3706: {
3707: if (!ya[i]) continue;
3708: T = u_FpX_div_by_X_x(Q, xa[i], p);
3709: inv = u_invmod(u_FpX_eval(T,xa[i], p), p);
3710: if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p)
3711: {
3712: dP = u_pol_comp(T, mulssmod(ya[i],inv,p), mulssmod(ya[i+1],inv,p), p);
3713: i++; /* x_i = -x_{i+1} */
3714: }
3715: else
1.2 ! noro 3716: dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p);
1.1 noro 3717: P = P? u_FpX_add(P, dP, p): dP;
3718: }
1.2 ! noro 3719: return P? gerepileupto(av, P): u_zeropol();
! 3720: }
! 3721:
! 3722: GEN
! 3723: FpV_polint(GEN xa, GEN ya, GEN p)
! 3724: {
! 3725: GEN inv,T,dP, P = NULL, Q = FpV_roots_to_pol(xa, p, 0);
! 3726: long i, n = lg(xa);
! 3727: gpmem_t av, lim;
! 3728: av = avma; lim = stack_lim(av,2);
! 3729: for (i=1; i<n; i++)
! 3730: {
! 3731: if (!signe(ya[i])) continue;
! 3732: T = FpX_div_by_X_x(Q, (GEN)xa[i], p);
! 3733: inv = mpinvmod(FpX_eval(T,(GEN)xa[i], p), p);
! 3734: if (i < n-1 && egalii(addii((GEN)xa[i], (GEN)xa[i+1]), p))
! 3735: {
! 3736: dP = pol_comp(T, muliimod((GEN)ya[i], inv,p),
! 3737: muliimod((GEN)ya[i+1],inv,p));
! 3738: i++; /* x_i = -x_{i+1} */
! 3739: }
! 3740: else
! 3741: dP = FpX_Fp_mul(T, muliimod((GEN)ya[i],inv,p), p);
! 3742: P = P? FpX_add(P, dP, p): dP;
! 3743: if (low_stack(lim, stack_lim(av,2)))
! 3744: {
! 3745: if (DEBUGMEM>1) err(warnmem,"FpV_polint");
! 3746: if (!P) avma = av; else P = gerepileupto(av, P);
! 3747: }
! 3748: }
! 3749: return P? P: zeropol(0);
1.1 noro 3750: }
3751:
3752: static void
3753: u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p)
3754: {
3755: GEN T,Q = u_FpV_roots_to_pol(xa, p);
3756: GEN dP = NULL, P = NULL;
3757: GEN dP0 = NULL, P0= NULL;
3758: GEN dP1 = NULL, P1= NULL;
3759: long i, n = lg(xa);
3760: ulong inv;
3761: for (i=1; i<n; i++)
3762: {
3763: T = u_FpX_div_by_X_x(Q, xa[i], p);
3764: inv = u_invmod(u_FpX_eval(T,xa[i], p), p);
3765:
1.2 ! noro 3766: if (ya[i])
! 3767: {
! 3768: dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p);
! 3769: P = P ? u_FpX_add(P , dP , p): dP;
! 3770: }
! 3771: if (C0[i])
1.1 noro 3772: {
1.2 ! noro 3773: dP0= u_FpX_Fp_mul(T, mulssmod(C0[i],inv,p), p);
! 3774: P0= P0? u_FpX_add(P0, dP0, p): dP0;
1.1 noro 3775: }
1.2 ! noro 3776: if (C1[i])
1.1 noro 3777: {
1.2 ! noro 3778: dP1= u_FpX_Fp_mul(T, mulssmod(C1[i],inv,p), p);
! 3779: P1= P1? u_FpX_add(P1, dP1, p): dP1;
1.1 noro 3780: }
3781: }
1.2 ! noro 3782: ya[1] = (long) (P ? P : u_zeropol());
! 3783: C0[1] = (long) (P0? P0: u_zeropol());
! 3784: C1[1] = (long) (P1? P1: u_zeropol());
1.1 noro 3785: }
3786:
3787: /* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
3788: * Return 0 in case of degree drop. */
3789: static GEN
1.2 ! noro 3790: u_vec_FpX_eval(GEN b, ulong x, ulong p)
1.1 noro 3791: {
3792: GEN z;
3793: long i, lb = lgef(b);
3794: ulong leadz = u_FpX_eval((GEN)b[lb-1], x, p);
1.2 ! noro 3795: if (!leadz) return u_zeropol();
1.1 noro 3796:
1.2 ! noro 3797: z = u_getpol(lb-3);
1.1 noro 3798: for (i=2; i<lb-1; i++)
3799: z[i] = u_FpX_eval((GEN)b[i], x, p);
3800: z[i] = leadz; return z;
3801: }
3802:
3803: /* as above, but don't care about degree drop */
3804: static GEN
1.2 ! noro 3805: u_vec_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop)
1.1 noro 3806: {
3807: GEN z;
3808: long i, lb = lgef(b);
1.2 ! noro 3809: z = u_getpol(lb-3);
1.1 noro 3810: for (i=2; i<lb; i++)
3811: z[i] = u_FpX_eval((GEN)b[i], x, p);
3812: z = u_normalizepol(z, lb);
3813: *drop = lb - lgef(z);
3814: return z;
3815: }
3816:
1.2 ! noro 3817: static GEN
! 3818: vec_FpX_eval_gen(GEN b, GEN x, GEN p, int *drop)
! 3819: {
! 3820: GEN z;
! 3821: long i, lb = lgef(b);
! 3822: z = cgetg(lb, t_POL);
! 3823: z[1] = b[1];
! 3824: for (i=2; i<lb; i++)
! 3825: z[i] = (long)FpX_eval((GEN)b[i], x, p);
! 3826: z = normalizepol_i(z, lb);
! 3827: *drop = lb - lgef(z);
! 3828: return z;
! 3829: }
! 3830:
1.1 noro 3831: /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1.2 ! noro 3832: * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
! 3833: * N_2(A) = sqrt(sum (N_1(Ai))^2)
! 3834: * Return e such that Res(A, B) < 2^e */
1.1 noro 3835: ulong
3836: ZY_ZXY_ResBound(GEN A, GEN B)
3837: {
1.2 ! noro 3838: gpmem_t av = avma;
1.1 noro 3839: GEN a = gzero, b = gzero, run = realun(DEFAULTPREC);
3840: long i , lA = lgef(A), lB = lgef(B);
3841: for (i=2; i<lA; i++) a = addii(a, sqri((GEN)A[i]));
3842: for (i=2; i<lB; i++)
3843: {
3844: GEN t = (GEN)B[i];
3845: if (typ(t) == t_POL) t = gnorml1(t, 0);
3846: b = addii(b, sqri(t));
3847: }
1.2 ! noro 3848: a = mulir(a,run);
! 3849: b = mulir(b,run);
! 3850: b = gmul(gpowgs(a, degpol(B)), gpowgs(b, degpol(A)));
1.1 noro 3851: avma = av; return 1 + (gexpo(b)>>1);
3852: }
3853:
1.2 ! noro 3854: /* return Res(a(Y), b(n,Y)) over Fp. la = leading_term(a) [for efficiency] */
! 3855: static ulong
! 3856: u_FpX_resultant_after_eval(GEN a, GEN b, ulong n, ulong p, ulong la)
! 3857: {
! 3858: int drop;
! 3859: GEN ev = u_vec_FpX_eval_gen(b, n, p, &drop);
! 3860: ulong r = u_FpX_resultant(a, ev, p);
! 3861: if (drop && la != 1) r = mulssmod(r, powuumod(la, drop,p),p);
! 3862: return r;
! 3863: }
! 3864: static GEN
! 3865: FpX_resultant_after_eval(GEN a, GEN b, GEN n, GEN p, GEN la)
! 3866: {
! 3867: int drop;
! 3868: GEN ev = vec_FpX_eval_gen(b, n, p, &drop);
! 3869: GEN r = FpX_resultant(a, ev, p);
! 3870: if (drop && !gcmp1(la)) r = muliimod(r, powgumod(la, drop,p),p);
! 3871: return r;
! 3872: }
! 3873:
! 3874: /* assume dres := deg(Res_Y(a,b), X) <= deg(a,Y) * deg(b,X) < p */
! 3875: static GEN
! 3876: u_FpY_FpXY_resultant(GEN a, GEN b, ulong p, long dres)
! 3877: {
! 3878: ulong la;
! 3879: long i,n,nmax;
! 3880: GEN x,y;
! 3881:
! 3882: nmax = (dres+1)>>1;
! 3883: la = (ulong)leading_term(a);
! 3884: x = cgetg(dres+2, t_VECSMALL);
! 3885: y = cgetg(dres+2, t_VECSMALL);
! 3886: /* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X),
! 3887: * where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */
! 3888: for (i=0,n = 1; n <= nmax; n++)
! 3889: {
! 3890: i++; x[i] = n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
! 3891: i++; x[i] = p-n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
! 3892: }
! 3893: if (i == dres)
! 3894: {
! 3895: i++; x[i] = 0; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
! 3896: }
! 3897: return u_FpV_polint(x,y, p);
! 3898: }
! 3899:
! 3900: /* x^n mod p */
! 3901: ulong
! 3902: u_powmod(ulong x, long n, ulong p)
! 3903: {
! 3904: ulong y = 1, z;
! 3905: long m;
! 3906:
! 3907: if (n < 0) { n = -n; x = u_invmod(x, p); }
! 3908: m = n;
! 3909: z = x;
! 3910: for (;;)
! 3911: {
! 3912: if (m&1) y = mulssmod(y,z, p);
! 3913: m >>= 1; if (!m) return y;
! 3914: z = mulssmod(z,z, p);
! 3915: }
! 3916: }
! 3917:
! 3918: /* x^n mod p, assume n > 0 */
! 3919: static GEN
! 3920: u_FpX_pow(GEN x, long n, ulong p)
! 3921: {
! 3922: GEN y = u_scalarpol(1), z;
! 3923: long m;
! 3924: m = n;
! 3925: z = x;
! 3926: for (;;)
! 3927: {
! 3928: if (m&1) y = u_FpX_mul(y,z, p);
! 3929: m >>= 1; if (!m) return y;
! 3930: z = u_FpX_mul(z,z, p);
! 3931: }
! 3932: }
! 3933:
! 3934: static GEN
! 3935: u_FpXX_pseudorem(GEN x, GEN y, ulong p)
! 3936: {
! 3937: long vx = varn(x), dx, dy, dz, i, lx, dp;
! 3938: gpmem_t av = avma, av2, lim;
! 3939:
! 3940: if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");
! 3941: (void)new_chunk(2);
! 3942: dx=degpol(x); x = revpol(x);
! 3943: dy=degpol(y); y = revpol(y); dz=dx-dy; dp = dz+1;
! 3944: av2 = avma; lim = stack_lim(av2,1);
! 3945: for (;;)
! 3946: {
! 3947: x[0] = (long)u_FpX_neg((GEN)x[0], p); dp--;
! 3948: for (i=1; i<=dy; i++)
! 3949: x[i] = (long)u_FpX_add( u_FpX_mul((GEN)y[0], (GEN)x[i], p),
! 3950: u_FpX_mul((GEN)x[0], (GEN)y[i], p), p );
! 3951: for ( ; i<=dx; i++)
! 3952: x[i] = (long)u_FpX_mul((GEN)y[0], (GEN)x[i], p);
! 3953: do { x++; dx--; } while (dx >= 0 && !signe((GEN)x[0]));
! 3954: if (dx < dy) break;
! 3955: if (low_stack(lim,stack_lim(av2,1)))
! 3956: {
! 3957: if(DEBUGMEM>1) err(warnmem,"pseudorem dx = %ld >= %ld",dx,dy);
! 3958: gerepilemanycoeffs(av2,x,dx+1);
! 3959: }
! 3960: }
! 3961: if (dx < 0) return u_zeropol();
! 3962: lx = dx+3; x -= 2;
! 3963: x[0]=evaltyp(t_POL) | evallg(lx);
! 3964: x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);
! 3965: x = revpol(x) - 2;
! 3966: if (dp)
! 3967: { /* multiply by y[0]^dp [beware dummy vars from FpY_FpXY_resultant] */
! 3968: GEN t = u_FpX_pow((GEN)y[0], dp, p);
! 3969: for (i=2; i<lx; i++)
! 3970: x[i] = (long)u_FpX_mul((GEN)x[i], t, p);
! 3971: }
! 3972: return gerepilecopy(av, x);
! 3973: }
! 3974:
! 3975: static GEN
! 3976: u_FpX_divexact(GEN x, GEN y, ulong p)
! 3977: {
! 3978: long i, l;
! 3979: GEN z;
! 3980: if (degpol(y) == 0)
! 3981: {
! 3982: ulong t = (ulong)y[2];
! 3983: if (t == 1) return x;
! 3984: t = u_invmod(t, p);
! 3985: l = lgef(x); z = cgetg(l, t_POL); z[1] = x[1];
! 3986: for (i=2; i<l; i++) z[i] = (long)u_FpX_Fp_mul((GEN)x[i],t,p);
! 3987: }
! 3988: else
! 3989: {
! 3990: l = lgef(x); z = cgetg(l, t_POL); z[1] = x[1];
! 3991: for (i=2; i<l; i++) z[i] = (long)u_FpX_div((GEN)x[i],y,p);
! 3992: }
! 3993: return z;
! 3994: }
! 3995:
! 3996: static GEN
! 3997: u_FpYX_subres(GEN u, GEN v, ulong p)
! 3998: {
! 3999: gpmem_t av = avma, av2, lim;
! 4000: long degq,dx,dy,du,dv,dr,signh;
! 4001: GEN z,g,h,r,p1;
! 4002:
! 4003: dx=degpol(u); dy=degpol(v); signh=1;
! 4004: if (dx < dy)
! 4005: {
! 4006: swap(u,v); lswap(dx,dy);
! 4007: if (both_odd(dx, dy)) signh = -signh;
! 4008: }
! 4009: if (dy < 0) return gzero;
! 4010: if (dy==0) return gerepileupto(av, u_FpX_pow((GEN)v[2],dx,p));
! 4011:
! 4012: g = h = u_scalarpol(1); av2 = avma; lim = stack_lim(av2,1);
! 4013: for(;;)
! 4014: {
! 4015: r = u_FpXX_pseudorem(u,v,p); dr = lgef(r);
! 4016: if (dr == 2) { avma = av; return gzero; }
! 4017: du = degpol(u); dv = degpol(v); degq = du-dv;
! 4018: u = v; p1 = g; g = leading_term(u);
! 4019: switch(degq)
! 4020: {
! 4021: case 0: break;
! 4022: case 1:
! 4023: p1 = u_FpX_mul(h,p1, p); h = g; break;
! 4024: default:
! 4025: p1 = u_FpX_mul(u_FpX_pow(h,degq,p), p1, p);
! 4026: h = u_FpX_div(u_FpX_pow(g,degq,p), u_FpX_pow(h,degq-1,p), p);
! 4027: }
! 4028: if (both_odd(du,dv)) signh = -signh;
! 4029: v = u_FpX_divexact(r, p1, p);
! 4030: if (dr==3) break;
! 4031: if (low_stack(lim,stack_lim(av2,1)))
! 4032: {
! 4033: if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);
! 4034: gerepileall(av2,4, &u, &v, &g, &h);
! 4035: }
! 4036: }
! 4037: z = (GEN)v[2];
! 4038: if (dv > 1) z = u_FpX_div(u_FpX_pow(z,dv,p), u_FpX_pow(h,dv-1,p), p);
! 4039: if (signh < 0) z = u_FpX_neg(z,p);
! 4040: return gerepileupto(av, z);
! 4041: }
! 4042:
! 4043: /* return a t_POL (in dummy variable 0) whose coeffs are the coeffs of b,
! 4044: * in variable v. This is an incorrect PARI object if initially varn(b) < v.
! 4045: * We could return a vector of coeffs, but it is convenient to have degpol()
! 4046: * and friends available. Even in that case, it will behave nicely with all
! 4047: * functions treating a polynomial as a vector of coeffs (eg poleval).
! 4048: * FOR INTERNAL USE! */
! 4049: GEN
! 4050: swap_vars(GEN b0, long v)
! 4051: {
! 4052: long i, n = poldegree(b0, v);
! 4053: GEN b = cgetg(n+3, t_POL), x = b + 2;
! 4054: b[1] = evalsigne(1) | evallgef(n+3) | evalvarn(v);
! 4055: for (i=0; i<=n; i++) x[i] = (long)polcoeff_i(b0, i, v);
! 4056: return b;
! 4057: }
! 4058:
! 4059: /* assume varn(b) < varn(a) */
! 4060: GEN
! 4061: FpY_FpXY_resultant(GEN a, GEN b0, GEN p)
! 4062: {
! 4063: long i,n,dres,nmax, vX = varn(b0), vY = varn(a);
! 4064: GEN la,x,y,b = swap_vars(b0, vY);
! 4065:
! 4066: dres = degpol(a)*degpol(b0);
! 4067: if (OK_ULONG(p))
! 4068: {
! 4069: ulong pp = (ulong)p[2];
! 4070: long l = lgef(b);
! 4071:
! 4072: a = u_Fp_FpX(a, pp);
! 4073: for (i=2; i<l; i++)
! 4074: b[i] = (long)u_Fp_FpX((GEN)b[i], pp);
! 4075: if (dres >= pp)
! 4076: {
! 4077: l = lgef(a);
! 4078: a[0] = evaltyp(t_POL) | evallg(l);
! 4079: a[1] = evalsigne(1)|evalvarn(vY)|evallgef(l);
! 4080: for (i=2; i<l; i++)
! 4081: a[i] = (long)u_scalarpol(a[i]);
! 4082: x = u_FpYX_subres(a, b, pp);
! 4083: }
! 4084: else
! 4085: x = u_FpY_FpXY_resultant(a, b, pp, dres);
! 4086: return small_to_pol(x, vX);
! 4087: }
! 4088:
! 4089: nmax = (dres+1)>>1;
! 4090: la = leading_term(a);
! 4091: x = cgetg(dres+2, t_VEC);
! 4092: y = cgetg(dres+2, t_VEC);
! 4093: /* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X),
! 4094: * where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */
! 4095: for (i=0,n = 1; n <= nmax; n++)
! 4096: {
! 4097: i++; x[i] = lstoi(n);
! 4098: y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la);
! 4099: i++; x[i] = lsubis(p,n);
! 4100: y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la);
! 4101: }
! 4102: if (i == dres)
! 4103: {
! 4104: i++; x[i] = zero;
! 4105: y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la);
! 4106: }
! 4107: x = FpV_polint(x,y, p);
! 4108: setvarn(x, vX); return x;
! 4109: }
! 4110:
! 4111: /* check that theta(maxprime) - theta(27448) >= 2^bound */
! 4112: /* NB: theta(27449) ~ 27225.387, theta(x) > 0.98 x for x>7481
! 4113: * (Schoenfeld, 1976 for x > 1155901 + direct calculations) */
! 4114: static void
! 4115: check_theta(ulong bound)
! 4116: {
! 4117: ulong c = (ulong)ceil((bound * LOG2 + 27225.388) / 0.98);
! 4118: if (maxprime() < c)
! 4119: err(talker,"not enough precalculated primes: need primelimit ~ %lu", c);
! 4120: }
! 4121:
1.1 noro 4122: /* 0, 1, -1, 2, -2, ... */
4123: #define next_lambda(a) (a>0 ? -a : 1-a)
4124:
4125: /* If lambda = NULL, assume A in Z[Y], B in Q[Y][X], return Res_Y(A,B)
4126: * Otherwise, find a small lambda (start from *lambda, use the sequence above)
4127: * such that R(X) = Res_Y(A, B(X + lambda Y)) is squarefree, reset *lambda to
4128: * the chosen value and return R
4129: *
4130: * If LERS is non-NULL, set it to the last non-constant polynomial in the PRS */
4131: GEN
4132: ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LERS)
4133: {
1.2 ! noro 4134: int checksqfree = lambda? 1: 0, delvar = 0, first = 1, stable;
! 4135: ulong bound;
! 4136: gpmem_t av = avma, av2, lim;
1.1 noro 4137: long i,n, lb, dres = degpol(A)*degpol(B0), nmax = (dres+1)>>1;
4138: long vX = varn(B0), vY = varn(A); /* assume vX < vY */
4139: GEN x,y,dglist,cB,B,q,a,b,ev,H,H0,H1,Hp,H0p,H1p,C0,C1;
4140: byteptr d = diffptr + 3000;
4141: ulong p = 27449; /* p = prime(3000) */
4142:
4143: dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
4144: if (LERS)
4145: {
4146: if (!lambda) err(talker,"ZY_ZXY_resultant_all: LERS needs lambda");
4147: C0 = cgetg(dres+2, t_VECSMALL);
4148: C1 = cgetg(dres+2, t_VECSMALL);
4149: dglist = cgetg(dres+1, t_VECSMALL);
4150: }
4151: x = cgetg(dres+2, t_VECSMALL);
4152: y = cgetg(dres+2, t_VECSMALL);
4153: if (vY == MAXVARN)
4154: {
1.2 ! noro 4155: vY = fetch_var(); delvar = 1;
1.1 noro 4156: B0 = gsubst(B0, MAXVARN, polx[vY]);
4157: A = dummycopy(A); setvarn(A, vY);
4158: }
4159: cB = content(B0);
4160: if (typ(cB) == t_POL) cB = content(cB);
4161: if (gcmp1(cB)) cB = NULL; else B0 = gdiv(B0, cB);
4162:
4163: if (lambda)
4164: B = poleval(B0, gadd(polx[MAXVARN], gmulsg(*lambda, polx[vY])));
4165: else
4166: B = poleval(B0, polx[MAXVARN]);
4167: av2 = avma; lim = stack_lim(av,2);
4168:
4169: INIT:
4170: if (first) first = 0;
4171: else
4172: { /* lambda != NULL */
4173: avma = av2; *lambda = next_lambda(*lambda);
4174: if (DEBUGLEVEL>4) fprintferr("Starting with lambda = %ld\n",*lambda);
4175: B = poleval(B0, gadd(polx[MAXVARN], gmulsg(*lambda, polx[vY])));
4176: av2 = avma;
4177: }
4178:
4179: if (degpol(A) <= 3)
4180: { /* sub-resultant faster for small degrees */
4181: if (LERS)
4182: {
4183: H = subresall(A,B,&q);
4184: if (typ(q) != t_POL || degpol(q)!=1 || !ZX_is_squarefree(H)) goto INIT;
4185: H0 = (GEN)q[2]; if (typ(H0) == t_POL) setvarn(H0,vX);
4186: H1 = (GEN)q[3]; if (typ(H1) == t_POL) setvarn(H1,vX);
4187: }
4188: else
4189: {
4190: H = subres(A,B);
4191: if (checksqfree && !ZX_is_squarefree(H)) goto INIT;
4192: }
4193: goto END;
4194: }
4195:
1.2 ! noro 4196: /* make sure p large enough */
! 4197: while (p < (dres<<1)) NEXT_PRIME_VIADIFF(p,d);
! 4198:
1.1 noro 4199: H = H0 = H1 = NULL;
1.2 ! noro 4200: lb = lgef(B); b = u_getpol(degpol(B));
! 4201: bound = ZY_ZXY_ResBound(A, B);
1.1 noro 4202: if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound);
1.2 ! noro 4203: check_theta(bound);
1.1 noro 4204:
4205: for(;;)
4206: {
1.2 ! noro 4207: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 4208:
1.2 ! noro 4209: a = u_Fp_FpX(A, p);
1.1 noro 4210: for (i=2; i<lb; i++)
1.2 ! noro 4211: b[i] = (long)u_Fp_FpX((GEN)B[i], p);
1.1 noro 4212: if (LERS)
4213: {
4214: if (!b[lb-1] || degpol(a) < degpol(A)) continue; /* p | lc(A)lc(B) */
4215: if (checksqfree)
4216: { /* find degree list for generic Euclidean Remainder Sequence */
4217: int goal = min(degpol(a), degpol(b)); /* longest possible */
4218: for (n=1; n <= goal; n++) dglist[n] = 0;
4219: setlg(dglist, 1);
4220: for (n=0; n <= dres; n++)
4221: {
1.2 ! noro 4222: ev = u_vec_FpX_eval(b, n, p);
1.1 noro 4223: (void)u_FpX_resultant_all(a, ev, NULL, NULL, dglist, p);
4224: if (lg(dglist)-1 == goal) break;
4225: }
4226: /* last pol in ERS has degree > 1 ? */
4227: goal = lg(dglist)-1;
4228: if (degpol(B) == 1) { if (!goal) goto INIT; }
4229: else
4230: {
4231: if (goal <= 1) goto INIT;
4232: if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
4233: }
4234: if (DEBUGLEVEL>4)
4235: fprintferr("Degree list for ERS (trials: %ld) = %Z\n",n+1,dglist);
4236: }
4237:
4238: for (i=0,n = 0; i <= dres; n++)
4239: {
1.2 ! noro 4240: i++; ev = u_vec_FpX_eval(b, n, p);
1.1 noro 4241: x[i] = n;
4242: y[i] = u_FpX_resultant_all(a, ev, C0+i, C1+i, dglist, p);
4243: if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
4244: }
4245: u_FpV_polint_all(x,y,C0,C1, p);
4246: Hp = (GEN)y[1];
4247: H0p= (GEN)C0[1];
4248: H1p= (GEN)C1[1];
4249: }
4250: else
1.2 ! noro 4251: { /* cf u_FpXY_resultant */
1.1 noro 4252: ulong la = (ulong)leading_term(a);
4253: /* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X),
4254: * where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */
4255: for (i=0,n = 1; n <= nmax; n++)
4256: {
1.2 ! noro 4257: i++; x[i] = n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
! 4258: i++; x[i] = p-n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
1.1 noro 4259: }
4260: if (i == dres)
4261: {
1.2 ! noro 4262: i++; x[i] = 0; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la);
1.1 noro 4263: }
4264: Hp = u_FpV_polint(x,y, p);
4265: }
4266: if (!H && degpol(Hp) != dres) continue;
4267: if (checksqfree) {
4268: if (!u_FpX_is_squarefree(Hp, p)) goto INIT;
4269: if (DEBUGLEVEL>4) fprintferr("Final lambda = %ld\n",*lambda);
4270: checksqfree = 0;
4271: }
4272:
4273: if (!H)
4274: { /* initialize */
4275: q = utoi(p); stable = 0;
4276: H = ZX_init_CRT(Hp, p,vX);
4277: if (LERS) {
4278: H0= ZX_init_CRT(H0p, p,vX);
4279: H1= ZX_init_CRT(H1p, p,vX);
4280: }
4281: }
4282: else
4283: {
4284: GEN qp = muliu(q,p);
1.2 ! noro 4285: stable = ZX_incremental_CRT(&H, Hp, q,qp, p);
1.1 noro 4286: if (LERS) {
1.2 ! noro 4287: stable &= ZX_incremental_CRT(&H0,H0p, q,qp, p);
! 4288: stable &= ZX_incremental_CRT(&H1,H1p, q,qp, p);
1.1 noro 4289: }
4290: q = qp;
4291: }
4292: /* could make it probabilistic for H ? [e.g if stable twice, etc]
4293: * Probabilistic anyway for H0, H1 */
4294: if (DEBUGLEVEL>5)
4295: msgtimer("resultant mod %ld (bound 2^%ld, stable=%ld)", p,expi(q),stable);
4296: if (stable && (ulong)expi(q) >= bound) break; /* DONE */
4297: if (low_stack(lim, stack_lim(av,2)))
4298: {
4299: GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1;
4300: if (DEBUGMEM>1) err(warnmem,"ZY_ZXY_resultant");
1.2 ! noro 4301: gerepilemany(av2,gptr,LERS? 4: 2); b = u_getpol(degpol(B));
1.1 noro 4302: }
4303: }
4304: END:
1.2 ! noro 4305: setvarn(H, vX); if (delvar) (void)delete_var();
1.1 noro 4306: if (cB) H = gmul(H, gpowgs(cB, degpol(A)));
4307: if (LERS)
4308: {
4309: GEN *gptr[2];
4310: GEN z = cgetg(3, t_VEC);
4311: z[1] = (long)H0;
4312: z[2] = (long)H1; *LERS = z;
4313: gptr[0]=&H; gptr[1]=LERS;
4314: gerepilemany(av, gptr, 2);
4315: return H;
4316: }
4317: if (!cB) H = gcopy(H);
4318: return gerepileupto(av, H);
4319: }
4320:
4321: GEN
1.2 ! noro 4322: ZY_ZXY_resultant(GEN A, GEN B, long *lambda)
1.1 noro 4323: {
1.2 ! noro 4324: return ZY_ZXY_resultant_all(A, B, lambda, NULL);
1.1 noro 4325: }
4326:
4327: /* If lambda = NULL, return caract(Mod(B, A)), A,B in Z[X].
4328: * Otherwise find a small lambda such that caract (Mod(B + lambda X, A)) is
4329: * squarefree */
4330: GEN
4331: ZX_caract_sqf(GEN A, GEN B, long *lambda, long v)
4332: {
1.2 ! noro 4333: gpmem_t av = avma;
1.1 noro 4334: GEN B0, R, a;
4335: long dB;
1.2 ! noro 4336: int delvar;
1.1 noro 4337:
4338: if (v < 0) v = 0;
4339: switch (typ(B))
4340: {
4341: case t_POL: dB = degpol(B); if (dB > 0) break;
4342: B = dB? (GEN)B[2]: gzero; /* fall through */
4343: default:
4344: if (lambda) { B = scalarpol(B,varn(A)); dB = 0; break;}
4345: return gerepileupto(av, gpowgs(gsub(polx[v], B), degpol(A)));
4346: }
1.2 ! noro 4347: delvar = 0;
1.1 noro 4348: if (varn(A) == 0)
4349: {
1.2 ! noro 4350: long v0 = fetch_var(); delvar = 1;
1.1 noro 4351: A = dummycopy(A); setvarn(A,v0);
4352: B = dummycopy(B); setvarn(B,v0);
4353: }
4354: B0 = cgetg(4, t_POL); B0[1] = evalsigne(1)|evalvarn(0)|evallgef(4);
4355: B0[2] = (long)gneg_i(B);
4356: B0[3] = un;
4357: R = ZY_ZXY_resultant(A, B0, lambda);
1.2 ! noro 4358: if (delvar) (void)delete_var();
1.1 noro 4359: setvarn(R, v); a = leading_term(A);
4360: if (!gcmp1(a)) R = gdiv(R, gpowgs(a, dB));
4361: return gerepileupto(av, R);
4362: }
4363:
1.2 ! noro 4364:
1.1 noro 4365: GEN
4366: ZX_caract(GEN A, GEN B, long v)
4367: {
1.2 ! noro 4368: return (degpol(A) < 16) ? caractducos(A,B,v): ZX_caract_sqf(A,B, NULL, v);
! 4369: }
! 4370:
! 4371: /* assume A integral, B in Q[v] */
! 4372: GEN
! 4373: QX_caract(GEN A, GEN B, long v)
! 4374: {
! 4375: gpmem_t av = avma;
! 4376: GEN cB, B0 = Q_primitive_part(B, &cB);
! 4377: GEN ch = ZX_caract(A, B0, v);
! 4378: if (cB)
! 4379: ch = gerepilecopy(av, rescale_pol(ch, cB));
! 4380: return ch;
1.1 noro 4381: }
4382:
4383: static GEN
4384: trivial_case(GEN A, GEN B)
4385: {
1.2 ! noro 4386: long d;
1.1 noro 4387: if (typ(A) == t_INT) return gpowgs(A, degpol(B));
1.2 ! noro 4388: d = degpol(A);
! 4389: if (d == 0) return trivial_case((GEN)A[2],B);
! 4390: if (d < 0) return gzero;
1.1 noro 4391: return NULL;
4392: }
4393:
1.2 ! noro 4394: /* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */
1.1 noro 4395: GEN
1.2 ! noro 4396: ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
1.1 noro 4397: {
1.2 ! noro 4398: ulong Hp, dp;
! 4399: gpmem_t av = avma, av2, lim;
! 4400: long degA;
1.1 noro 4401: int stable;
4402: GEN q, a, b, H;
4403: byteptr d = diffptr + 3000;
4404: ulong p = 27449; /* p = prime(3000) */
4405:
4406: if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H;
4407: q = H = NULL;
4408: av2 = avma; lim = stack_lim(av,2);
1.2 ! noro 4409: degA = degpol(A);
! 4410: if (!bound)
! 4411: {
! 4412: bound = ZY_ZXY_ResBound(A, B);
! 4413: if (bound > 50000)
! 4414: {
! 4415: GEN run = realun(MEDDEFAULTPREC);
! 4416: GEN Ar = gmul(A, run), Br = gmul(B, run);
! 4417: GEN R = subres(Ar,Br);
! 4418: bound = gexpo(R) + 1;
! 4419: }
! 4420: if (dB) bound -= (long)(mylog2(dB)*degA);
! 4421: }
1.1 noro 4422: if (DEBUGLEVEL>4) fprintferr("bound for resultant: 2^%ld\n",bound);
1.2 ! noro 4423: check_theta(bound);
1.1 noro 4424:
1.2 ! noro 4425: dp = 0; /* gcc -Wall */
1.1 noro 4426: for(;;)
4427: {
1.2 ! noro 4428: NEXT_PRIME_VIADIFF_CHECK(p,d);
! 4429: if (dB) { dp = smodis(dB, p); if (!dp) continue; }
! 4430:
! 4431: a = u_Fp_FpX(A, p);
! 4432: b = u_Fp_FpX(B, p);
1.1 noro 4433: Hp= u_FpX_resultant(a, b, p);
1.2 ! noro 4434: if (dp) Hp = mulssmod(Hp, u_powmod(dp, -degA, p), p);
1.1 noro 4435: if (!H)
4436: {
4437: stable = 0; q = utoi(p);
1.2 ! noro 4438: H = stoi(u_center(Hp, p, p>>1));
1.1 noro 4439: }
4440: else /* could make it probabilistic ??? [e.g if stable twice, etc] */
4441: {
4442: GEN qp = muliu(q,p);
4443: stable = Z_incremental_CRT(&H, Hp, q,qp, p);
4444: q = qp;
4445: }
4446: if (DEBUGLEVEL>5)
4447: msgtimer("resultant mod %ld (bound 2^%ld, stable = %d)",p,expi(q),stable);
4448: if (stable && (ulong)expi(q) >= bound) break; /* DONE */
4449: if (low_stack(lim, stack_lim(av,2)))
4450: {
4451: GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q;
4452: if (DEBUGMEM>1) err(warnmem,"ZX_resultant");
4453: gerepilemany(av2,gptr, 2);
4454: }
4455: }
4456: return gerepileuptoint(av, icopy(H));
4457: }
4458:
4459: GEN
1.2 ! noro 4460: ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
! 4461:
! 4462: GEN
! 4463: ZX_QX_resultant(GEN A, GEN B)
! 4464: {
! 4465: GEN c, d, n, R;
! 4466: gpmem_t av = avma;
! 4467: B = Q_primitive_part(B, &c);
! 4468: if (!c) return ZX_resultant(A,B);
! 4469: n = numer(c);
! 4470: d = denom(c); if (is_pm1(d)) d = NULL;
! 4471: R = ZX_resultant_all(A, B, d, 0);
! 4472: if (!is_pm1(n)) R = mulii(R, gpowgs(n, degpol(A)));
! 4473: return gerepileuptoint(av, R);
! 4474: }
1.1 noro 4475:
4476: /* assume x has integral coefficients */
4477: GEN
4478: ZX_disc_all(GEN x, ulong bound)
4479: {
1.2 ! noro 4480: gpmem_t av = avma;
! 4481: GEN l, d = ZX_resultant_all(x, derivpol(x), NULL, bound);
! 4482: l = leading_term(x); if (!gcmp1(l)) d = diviiexact(d,l);
1.1 noro 4483: if (degpol(x) & 2) d = negi(d);
4484: return gerepileuptoint(av,d);
4485: }
4486: GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); }
4487:
4488: int
4489: ZX_is_squarefree(GEN x)
4490: {
1.2 ! noro 4491: gpmem_t av = avma;
1.1 noro 4492: int d = (lgef(modulargcd(x,derivpol(x))) == 3);
4493: avma = av; return d;
4494: }
4495:
1.2 ! noro 4496: static GEN
! 4497: _gcd(GEN a, GEN b)
1.1 noro 4498: {
1.2 ! noro 4499: if (!a) a = gun;
! 4500: if (!b) b = gun;
! 4501: return ggcd(a,b);
1.1 noro 4502: }
4503:
4504: /* A0 and B0 in Q[X] */
4505: GEN
4506: modulargcd(GEN A0, GEN B0)
4507: {
1.2 ! noro 4508: GEN a,b,Hp,D,A,B,q,qp,H,g;
1.1 noro 4509: long m,n;
1.2 ! noro 4510: ulong p;
! 4511: gpmem_t av2, av = avma, avlim = stack_lim(av, 1);
1.1 noro 4512: byteptr d = diffptr;
4513:
4514: if ((typ(A0) | typ(B0)) !=t_POL) err(notpoler,"modulargcd");
4515: if (!signe(A0)) return gcopy(B0);
4516: if (!signe(B0)) return gcopy(A0);
1.2 ! noro 4517: A = primitive_part(A0, &a); check_pol_int(A, "modulargcd");
! 4518: B = primitive_part(B0, &b); check_pol_int(B, "modulargcd");
! 4519: D = _gcd(a,b);
! 4520: if (varn(A) != varn(B)) err(talker,"different variables in modulargcd");
! 4521:
1.1 noro 4522: /* A, B in Z[X] */
4523: g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */
4524: if (degpol(A) < degpol(B)) swap(A, B);
4525: n = 1 + degpol(B); /* > degree(gcd) */
4526:
4527: av2 = avma; H = NULL;
4528: d += 3000; /* 27449 = prime(3000) */
1.2 ! noro 4529: for(p = 27449; ;)
1.1 noro 4530: {
4531: if (!*d) err(primer1);
1.2 ! noro 4532: if (!umodiu(g,p)) goto repeat;
1.1 noro 4533:
1.2 ! noro 4534: a = u_Fp_FpX(A, p);
! 4535: b = u_Fp_FpX(B, p); Hp = u_FpX_gcd_i(a,b, p);
1.1 noro 4536: m = degpol(Hp);
4537: if (m == 0) { H = polun[varn(A0)]; break; } /* coprime. DONE */
1.2 ! noro 4538: if (m > n) goto repeat; /* p | Res(A/G, B/G). Discard */
1.1 noro 4539:
4540: if (is_pm1(g)) /* make sure lead(H) = g mod p */
4541: Hp = u_FpX_normalize(Hp, p);
4542: else
4543: {
4544: ulong t = mulssmod(umodiu(g, p), u_invmod(Hp[m+2],p), p);
1.2 ! noro 4545: Hp = u_FpX_Fp_mul(Hp, t, p);
1.1 noro 4546: }
4547: if (m < n)
4548: { /* First time or degree drop [all previous p were as above; restart]. */
4549: H = ZX_init_CRT(Hp,p,varn(A0));
1.2 ! noro 4550: q = utoi(p); n = m; goto repeat;
1.1 noro 4551: }
4552:
4553: qp = muliu(q,p);
1.2 ! noro 4554: if (ZX_incremental_CRT(&H, Hp, q, qp, p))
1.1 noro 4555: { /* H stable: check divisibility */
1.2 ! noro 4556: if (!is_pm1(g)) H = primpart(H);
! 4557: if (gcmp0(pseudorem(A,H)) && gcmp0(pseudorem(B,H))) break; /* DONE */
1.1 noro 4558:
4559: if (DEBUGLEVEL) fprintferr("modulargcd: trial division failed");
4560: }
4561: q = qp;
4562: if (low_stack(avlim, stack_lim(av,1)))
4563: {
4564: GEN *gptr[2]; gptr[0]=&H; gptr[1]=&q;
4565: if (DEBUGMEM>1) err(warnmem,"modulargcd");
4566: gerepilemany(av2,gptr,2);
4567: }
1.2 ! noro 4568: repeat:
! 4569: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 4570: }
4571: return gerepileupto(av, gmul(D,H));
4572: }
4573:
1.2 ! noro 4574: /* lift(1 / Mod(A,B)). B0 a t_POL, A0 a scalar or a t_POL. Rational coeffs */
1.1 noro 4575: GEN
1.2 ! noro 4576: QX_invmod(GEN A0, GEN B0)
1.1 noro 4577: {
4578: GEN a,b,D,A,B,q,qp,Up,Vp,U,V,res;
4579: long stable;
1.2 ! noro 4580: ulong p;
! 4581: gpmem_t av2, av = avma, avlim = stack_lim(av, 1);
1.1 noro 4582: byteptr d = diffptr;
4583:
1.2 ! noro 4584: if (typ(B0) != t_POL) err(notpoler,"QX_invmod");
1.1 noro 4585: if (typ(A0) != t_POL)
4586: {
4587: if (is_scalar_t(typ(A0))) return ginv(A0);
1.2 ! noro 4588: err(notpoler,"QX_invmod");
1.1 noro 4589: }
1.2 ! noro 4590: if (degpol(A0) < 15) return ginvmod(A0,B0);
! 4591: A = primitive_part(A0, &D);
! 4592: B = primpart(B0);
1.1 noro 4593: /* A, B in Z[X] */
4594: av2 = avma; U = NULL;
4595: d += 3000; /* 27449 = prime(3000) */
1.2 ! noro 4596: for(p = 27449; ; )
1.1 noro 4597: {
4598: if (!*d) err(primer1);
1.2 ! noro 4599: a = u_Fp_FpX(A, p);
! 4600: b = u_Fp_FpX(B, p);
1.1 noro 4601: /* if p | Res(A/G, B/G), discard */
1.2 ! noro 4602: if (!u_FpX_extresultant(b,a,p, &Vp,&Up)) goto repeat;
1.1 noro 4603:
4604: if (!U)
4605: { /* First time */
4606: U = ZX_init_CRT(Up,p,varn(A0));
4607: V = ZX_init_CRT(Vp,p,varn(A0));
1.2 ! noro 4608: q = utoi(p); goto repeat;
1.1 noro 4609: }
1.2 ! noro 4610: if (DEBUGLEVEL>5) msgtimer("QX_invmod: mod %ld (bound 2^%ld)", p,expi(q));
1.1 noro 4611: qp = muliu(q,p);
1.2 ! noro 4612: stable = ZX_incremental_CRT(&U, Up, q,qp, p);
! 4613: stable &= ZX_incremental_CRT(&V, Vp, q,qp, p);
1.1 noro 4614: if (stable)
4615: { /* all stable: check divisibility */
4616: res = gadd(gmul(A,U), gmul(B,V));
4617: if (degpol(res) == 0) break; /* DONE */
1.2 ! noro 4618: if (DEBUGLEVEL) fprintferr("QX_invmod: char 0 check failed");
1.1 noro 4619: }
4620: q = qp;
4621: if (low_stack(avlim, stack_lim(av,1)))
4622: {
4623: GEN *gptr[3]; gptr[0]=&q; gptr[1]=&U; gptr[2]=&V;
1.2 ! noro 4624: if (DEBUGMEM>1) err(warnmem,"QX_invmod");
1.1 noro 4625: gerepilemany(av2,gptr,3);
4626: }
1.2 ! noro 4627: repeat:
! 4628: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 4629: }
1.2 ! noro 4630: D = D? gmul(D,res): res;
1.1 noro 4631: return gerepileupto(av, gdiv(U,D));
4632: }
4633:
1.2 ! noro 4634: /* irreducible (unitary) polynomial of degree n over Fp */
! 4635: GEN
! 4636: ffinit_rand(GEN p,long n)
! 4637: {
! 4638: gpmem_t av = avma;
! 4639: GEN pol;
! 4640:
! 4641: for(;; avma = av)
! 4642: {
! 4643: pol = gadd(gpowgs(polx[0],n), FpX_rand(n-1,0, p));
! 4644: if (FpX_is_irred(pol, p)) break;
! 4645: }
! 4646: return pol;
! 4647: }
! 4648:
! 4649: GEN
! 4650: FpX_direct_compositum(GEN A, GEN B, GEN p)
! 4651: {
! 4652: GEN C,a,b,x;
! 4653: a = dummycopy(A); setvarn(a, MAXVARN);
! 4654: b = dummycopy(B); setvarn(b, MAXVARN);
! 4655: x = gadd(polx[0], polx[MAXVARN]);
! 4656: C = FpY_FpXY_resultant(a, poleval(b,x),p);
! 4657: return C;
! 4658: }
! 4659:
! 4660: GEN
! 4661: FpX_compositum(GEN A, GEN B, GEN p)
! 4662: {
! 4663: GEN C, a,b;
! 4664: long k;
! 4665:
! 4666: a = dummycopy(A); setvarn(a, MAXVARN);
! 4667: b = dummycopy(B); setvarn(b, MAXVARN);
! 4668: for (k = 1;; k = next_lambda(k))
! 4669: {
! 4670: GEN x = gadd(polx[0], gmulsg(k, polx[MAXVARN]));
! 4671: C = FpY_FpXY_resultant(a, poleval(b,x),p);
! 4672: if (FpX_is_squarefree(C, p)) break;
! 4673: }
! 4674: return C;
! 4675: }
! 4676:
! 4677: /* return an extension of degree 2^l of F_2, assume l > 0
! 4678: * using Adleman-Lenstra Algorithm.
! 4679: * Not stack clean. */
! 4680: static GEN
! 4681: f2init(long l)
! 4682: {
! 4683: long i;
! 4684: GEN a, q, T, S;
! 4685:
! 4686: if (l == 1) return cyclo(3, MAXVARN);
! 4687:
! 4688: a = gun;
! 4689: S = coefs_to_pol(4, gun,gun,gzero,gzero); /* #(#^2 + #) */
! 4690: setvarn(S, MAXVARN);
! 4691: q = coefs_to_pol(3, gun,gun, S); /* X^2 + X + #(#^2+#) */
! 4692:
! 4693: /* x^4+x+1, irred over F_2, minimal polynomial of a root of q */
! 4694: T = coefs_to_pol(5, gun,gzero,gzero,gun,gun);
! 4695:
! 4696: for (i=2; i<l; i++)
! 4697: { /* q = X^2 + X + a(#) irred. over K = F2[#] / (T(#))
! 4698: * ==> X^2 + X + a(#) b irred. over K for any root b of q
! 4699: * ==> X^2 + X + (b^2+b)b */
! 4700: setvarn(T, MAXVARN);
! 4701: T = FpY_FpXY_resultant(T, q, gdeux);
! 4702: /* T = minimal polynomial of b over F2 */
! 4703: }
! 4704: return T;
! 4705: }
! 4706:
! 4707: /*Check if subcyclo(n,l,0) is irreducible modulo p*/
! 4708: static long
! 4709: fpinit_check(GEN p, long n, long l)
! 4710: {
! 4711: gpmem_t ltop=avma;
! 4712: long q,o;
! 4713: if (!isprime(stoi(n))) {avma=ltop; return 0;}
! 4714: q = smodis(p,n);
! 4715: if (!q) {avma=ltop; return 0;}
! 4716: o = itos(order(gmodulss(q,n)));
! 4717: avma = ltop;
! 4718: return ( cgcd((n-1)/o,l) == 1 );
! 4719: }
! 4720:
! 4721: /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
! 4722: * Return an irreducible polynomial of degree l over F_p.
! 4723: * This a variant of an algorithm of Adleman and Lenstra
! 4724: * "Finding irreducible polynomials over finite fields",
! 4725: * ACM, 1986 (5) 350--355
! 4726: * Not stack clean.
! 4727: */
! 4728: static GEN
! 4729: fpinit(GEN p, long l)
! 4730: {
! 4731: ulong n = 1+l, k = 1;
! 4732: while (!fpinit_check(p,n,l)) { n += l; k++; }
! 4733: if (DEBUGLEVEL>=4)
! 4734: fprintferr("FFInit: using subcyclo(%ld, %ld)\n",n,l);
! 4735: return FpX_red(subcyclo(n,l,0),p);
! 4736: }
! 4737:
! 4738: GEN
! 4739: ffinit_fact(GEN p, long n)
! 4740: {
! 4741: gpmem_t ltop=avma;
! 4742: GEN F; /* vecsmall */
! 4743: GEN P; /* pol */
! 4744: long i;
! 4745: F = decomp_primary_small(n);
! 4746: /* If n is even, then F[1] is 2^bfffo(n)*/
! 4747: if (!(n&1) && egalii(p, gdeux))
! 4748: P = f2init(vals(n));
! 4749: else
! 4750: P = fpinit(p, F[1]);
! 4751: for (i = 2; i < lg(F); ++i)
! 4752: P = FpX_direct_compositum(fpinit(p, F[i]), P, p);
! 4753: return gerepileupto(ltop,FpX(P,p));
! 4754: }
! 4755:
! 4756: GEN
! 4757: ffinit_nofact(GEN p, long n)
! 4758: {
! 4759: gpmem_t av = avma;
! 4760: GEN P,Q=NULL;
! 4761: if (lgefint(p)==3)
! 4762: {
! 4763: ulong lp=p[2], q;
! 4764: long v=svaluation(n,lp,&q);
! 4765: if (v>0)
! 4766: {
! 4767: if (lp==2)
! 4768: Q=f2init(v);
! 4769: else
! 4770: Q=fpinit(p,n/q);
! 4771: n=q;
! 4772: }
! 4773: }
! 4774: if (n==1)
! 4775: P=Q;
! 4776: else
! 4777: {
! 4778: P = fpinit(p, n);
! 4779: if (Q) P = FpX_direct_compositum(P, Q, p);
! 4780: }
! 4781: return gerepileupto(av, FpX(P,p));
! 4782: }
! 4783:
! 4784: GEN
! 4785: ffinit(GEN p, long n, long v)
! 4786: {
! 4787: gpmem_t ltop=avma;
! 4788: GEN P;
! 4789: if (n <= 0) err(talker,"non positive degree in ffinit");
! 4790: if (typ(p) != t_INT) err(typeer, "ffinit");
! 4791: if (v < 0) v = 0;
! 4792: if (n == 1) return FpX(polx[v],p);
! 4793: /*If we are in a easy case just use cyclo*/
! 4794: if (fpinit_check(p, n + 1, n))
! 4795: return gerepileupto(ltop,FpX(cyclo(n + 1, v),p));
! 4796: if (lgefint(p)-2<BITS_IN_LONG-bfffo(n))
! 4797: P=ffinit_fact(p,n);
! 4798: else
! 4799: P=ffinit_nofact(p,n);
! 4800: setvarn(P, v);
! 4801: return P;
! 4802: }
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