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