Annotation of OpenXM_contrib/pari-2.2/src/basemath/bibli2.c, Revision 1.2
1.2 ! noro 1: /* $Id: bibli2.c,v 1.38 2002/09/07 21:45:57 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: /** BIBLIOTHEQUE MATHEMATIQUE **/
19: /** (deuxieme partie) **/
20: /** **/
21: /********************************************************************/
22: #include "pari.h"
1.2 ! noro 23: extern GEN _ei(long n, long i);
1.1 noro 24:
25: /********************************************************************/
26: /** **/
27: /** DEVELOPPEMENTS LIMITES **/
28: /** **/
29: /********************************************************************/
30:
31: GEN
1.2 ! noro 32: tayl(GEN x, long v, long precS)
1.1 noro 33: {
1.2 ! noro 34: long i, vx = gvar9(x);
! 35: gpmem_t tetpil, av=avma;
1.1 noro 36: GEN p1,y;
37:
38: if (v <= vx)
39: {
1.2 ! noro 40: long p1[] = { evaltyp(t_SER)|_evallg(2), 0 };
! 41: p1[1] = evalvalp(precS) | evalvarn(v);
1.1 noro 42: return gadd(p1,x);
43: }
44: p1=cgetg(v+2,t_VEC);
45: for (i=0; i<v; i++) p1[i+1]=lpolx[i];
46: p1[vx+1]=lpolx[v]; p1[v+1]=lpolx[vx];
1.2 ! noro 47: y = tayl(changevar(x,p1), vx,precS); tetpil=avma;
1.1 noro 48: return gerepile(av,tetpil, changevar(y,p1));
49: }
50:
51: GEN
52: grando0(GEN x, long n, long do_clone)
53: {
54: long m, v, tx=typ(x);
55:
56: if (gcmp0(x)) err(talker,"zero argument in O()");
57: if (tx == t_INT)
58: {
59: if (!gcmp1(x)) /* bug 3 + O(1). We suppose x is a truc() */
60: {
1.2 ! noro 61: if (do_clone) x = gclone(x);
! 62: return padiczero(x,n);
1.1 noro 63: }
64: v=0; m=0; /* 1 = x^0 */
65: }
66: else
67: {
68: if (tx != t_POL && ! is_rfrac_t(tx))
69: err(talker,"incorrect argument in O()");
1.2 ! noro 70: v = gvar(x); if (v > MAXVARN) err(talker,"incorrect object in O()");
! 71: m = n*gval(x,v);
1.1 noro 72: }
73: return zeroser(v,m);
74: }
75:
76: /*******************************************************************/
77: /** **/
78: /** SPECIAL POLYNOMIALS **/
79: /** **/
80: /*******************************************************************/
81: #ifdef LONG_IS_64BIT
82: # define SQRTVERYBIGINT 3037000500 /* ceil(sqrt(VERYBIGINT)) */
83: #else
84: # define SQRTVERYBIGINT 46341
85: #endif
86:
87: /* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
88: * T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
89: * where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
90: * and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
91: GEN
92: tchebi(long n, long v) /* Assume 4*n < VERYBIGINT */
93: {
1.2 ! noro 94: long k, l;
! 95: gpmem_t av;
1.1 noro 96: GEN q,a,r;
97:
98: if (v<0) v = 0;
1.2 ! noro 99: if (n < 0) err(talker,"negative degree in tchebi");
1.1 noro 100: if (n==0) return polun[v];
101: if (n==1) return polx[v];
102:
103: q = cgetg(n+3, t_POL); r = q + n+2;
104: a = shifti(gun, n-1);
105: *r-- = (long)a;
106: *r-- = zero;
107: if (n < SQRTVERYBIGINT)
108: for (k=1,l=n; l>1; k++,l-=2)
109: {
110: av = avma;
111: a = divis(mulis(a, l*(l-1)), 4*k*(n-k));
112: a = gerepileuptoint(av, negi(a));
113: *r-- = (long)a;
114: *r-- = zero;
115: }
116: else
117: for (k=1,l=n; l>1; k++,l-=2)
118: {
119: av = avma;
120: a = mulis(mulis(a, l), l-1);
121: a = divis(divis(a, 4*k), n-k);
122: a = gerepileuptoint(av, negi(a));
123: *r-- = (long)a;
124: *r-- = zero;
125: }
126: q[1] = evalsigne(1) | evalvarn(v) | evallgef(n+3);
127: return q;
128: }
129:
130: GEN addshiftw(GEN x, GEN y, long d);
131: /* Legendre polynomial */
132: /* L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1) */
133: GEN
134: legendre(long n, long v)
135: {
1.2 ! noro 136: long m;
! 137: gpmem_t av, tetpil, lim;
1.1 noro 138: GEN p0,p1,p2;
139:
140: if (v<0) v = 0;
1.2 ! noro 141: if (n < 0) err(talker,"negative degree in legendre");
1.1 noro 142: if (n==0) return polun[v];
143: if (n==1) return polx[v];
144:
145: p0=polun[v]; av=avma; lim=stack_lim(av,2);
146: p1=gmul2n(polx[v],1);
147: for (m=1; m<n; m++)
148: {
149: p2 = addshiftw(gmulsg(4*m+2,p1), gmulsg(-4*m,p0), 1);
150: setvarn(p2,v);
151: p0 = p1; tetpil=avma; p1 = gdivgs(p2,m+1);
152: if (low_stack(lim, stack_lim(av,2)))
153: {
154: GEN *gptr[2];
155: if(DEBUGMEM>1) err(warnmem,"legendre");
156: p0=gcopy(p0); gptr[0]=&p0; gptr[1]=&p1;
157: gerepilemanysp(av,tetpil,gptr,2);
158: }
159: }
160: tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-n));
161: }
162:
163: /* cyclotomic polynomial */
164: GEN
165: cyclo(long n, long v)
166: {
1.2 ! noro 167: long d, q, m;
! 168: gpmem_t av=avma, tetpil;
1.1 noro 169: GEN yn,yd;
170:
171: if (n<=0) err(arither2);
172: if (v<0) v = 0;
173: yn = yd = polun[0];
174: for (d=1; d*d<=n; d++)
175: {
176: if (n%d) continue;
177: q=n/d;
178: m = mu(stoi(q));
179: if (m)
180: { /* y *= (x^d - 1) */
181: if (m>0) yn = addshiftw(yn, gneg(yn), d);
182: else yd = addshiftw(yd, gneg(yd), d);
183: }
184: if (q==d) break;
185: m = mu(stoi(d));
186: if (m)
187: { /* y *= (x^q - 1) */
188: if (m>0) yn = addshiftw(yn, gneg(yn), q);
189: else yd = addshiftw(yd, gneg(yd), q);
190: }
191: }
192: tetpil=avma; yn = gerepile(av,tetpil,gdeuc(yn,yd));
193: setvarn(yn,v); return yn;
194: }
195:
196: /* compute prod (L*x ± a[i]) */
197: GEN
198: roots_to_pol_intern(GEN L, GEN a, long v, int plus)
199: {
200: long i,k,lx = lg(a), code;
201: GEN p1,p2;
202: if (lx == 1) return polun[v];
203: p1 = cgetg(lx, t_VEC);
204: code = evalsigne(1)|evalvarn(v)|evallgef(5);
205: for (k=1,i=1; i<lx-1; i+=2)
206: {
207: p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
208: p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
209: p2[3] = ladd((GEN)a[i],(GEN)a[i+1]);
210: if (plus == 0) p2[3] = lneg((GEN)p2[3]);
211: p2[4] = (long)L; p2[1] = code;
212: }
213: if (i < lx)
214: {
215: p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
216: p2[1] = code = evalsigne(1)|evalvarn(v)|evallgef(4);
217: p2[2] = plus? a[i]: lneg((GEN)a[i]);
218: p2[3] = (long)L;
219: }
220: setlg(p1, k); return divide_conquer_prod(p1, gmul);
221: }
222:
223: GEN
224: roots_to_pol(GEN a, long v)
225: {
226: return roots_to_pol_intern(gun,a,v,0);
227: }
228:
229: /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
230: GEN
231: roots_to_pol_r1r2(GEN a, long r1, long v)
232: {
233: long i,k,lx = lg(a), code;
234: GEN p1;
235: if (lx == 1) return polun[v];
236: p1 = cgetg(lx, t_VEC);
237: code = evalsigne(1)|evalvarn(v)|evallgef(5);
238: for (k=1,i=1; i<r1; i+=2)
239: {
240: GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
241: p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
242: p2[3] = lneg(gadd((GEN)a[i],(GEN)a[i+1]));
243: p2[4] = un; p2[1] = code;
244: }
245: if (i < r1+1)
246: p1[k++] = ladd(polx[v], gneg((GEN)a[i]));
247: for (i=r1+1; i<lx; i++)
248: {
249: GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
250: p2[2] = lnorm((GEN)a[i]);
251: p2[3] = lneg(gtrace((GEN)a[i]));
252: p2[4] = un; p2[1] = code;
253: }
254: setlg(p1, k); return divide_conquer_prod(p1, gmul);
255: }
256:
257: /********************************************************************/
258: /** **/
259: /** HILBERT & PASCAL MATRICES **/
260: /** **/
261: /********************************************************************/
262: GEN
263: mathilbert(long n) /* Hilbert matrix of order n */
264: {
265: long i,j;
266: GEN a,p;
267:
268: if (n<0) n = 0;
269: p = cgetg(n+1,t_MAT);
270: for (j=1; j<=n; j++)
271: {
272: p[j]=lgetg(n+1,t_COL);
273: for (i=1+(j==1); i<=n; i++)
274: {
275: a=cgetg(3,t_FRAC); a[1]=un; a[2]=lstoi(i+j-1);
276: coeff(p,i,j)=(long)a;
277: }
278: }
1.2 ! noro 279: if ( n ) mael(p,1,1)=un;
1.1 noro 280: return p;
281: }
282:
283: /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
284: GEN
285: matqpascal(long n, GEN q)
286: {
1.2 ! noro 287: long i, j, I;
! 288: gpmem_t av = avma;
1.1 noro 289: GEN m, *qpow = NULL; /* gcc -Wall */
290:
291: if (n<0) n = -1;
292: n++; m = cgetg(n+1,t_MAT);
293: for (j=1; j<=n; j++) m[j] = lgetg(n+1,t_COL);
294: if (q)
295: {
296: I = (n+1)/2;
297: if (I > 1) { qpow = (GEN*)new_chunk(I+1); qpow[2]=q; }
298: for (j=3; j<=I; j++) qpow[j] = gmul(q, qpow[j-1]);
299: }
300: for (i=1; i<=n; i++)
301: {
302: I = (i+1)/2; coeff(m,i,1)=un;
303: if (q)
304: {
305: for (j=2; j<=I; j++)
306: coeff(m,i,j) = ladd(gmul(qpow[j],gcoeff(m,i-1,j)), gcoeff(m,i-1,j-1));
307: }
308: else
309: {
310: for (j=2; j<=I; j++)
311: coeff(m,i,j) = laddii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
312: }
313: for ( ; j<=i; j++) coeff(m,i,j) = coeff(m,i,i+1-j);
314: for ( ; j<=n; j++) coeff(m,i,j) = zero;
315: }
316: return gerepilecopy(av, m);
317: }
318:
319: /********************************************************************/
320: /** **/
321: /** LAPLACE TRANSFORM (OF A SERIES) **/
322: /** **/
323: /********************************************************************/
324:
325: GEN
326: laplace(GEN x)
327: {
1.2 ! noro 328: gpmem_t av = avma;
1.1 noro 329: long i,l,ec;
330: GEN y,p1;
331:
332: if (typ(x)!=t_SER) err(talker,"not a series in laplace");
333: if (gcmp0(x)) return gcopy(x);
334:
335: ec = valp(x);
336: if (ec<0) err(talker,"negative valuation in laplace");
337: l=lg(x); y=cgetg(l,t_SER);
338: p1=mpfact(ec); y[1]=x[1];
339: for (i=2; i<l; i++)
340: {
341: y[i]=lmul(p1,(GEN)x[i]);
342: ec++; p1=mulsi(ec,p1);
343: }
344: return gerepilecopy(av,y);
345: }
346:
347: /********************************************************************/
348: /** **/
349: /** CONVOLUTION PRODUCT (OF TWO SERIES) **/
350: /** **/
351: /********************************************************************/
352:
353: GEN
354: convol(GEN x, GEN y)
355: {
356: long l,i,j,v, vx=varn(x), lx=lg(x), ly=lg(y), ex=valp(x), ey=valp(y);
357: GEN z;
358:
359: if (typ(x) != t_SER || typ(y) != t_SER)
360: err(talker,"not a series in convol");
361: if (gcmp0(x) || gcmp0(y))
362: err(talker,"zero series in convol");
363: if (varn(y) != vx)
364: err(talker,"different variables in convol");
365: v=ex; if (ey>v) v=ey;
366: l=ex+lx; i=ey+ly; if (i<l) l=i;
367: l -= v; if (l<3) err(talker,"non significant result in convol");
368: for (i=v+2; i < v+l; i++)
369: if (!gcmp0((GEN)x[i-ex]) && !gcmp0((GEN)y[i-ey])) { i++; break; }
370: if (i == l+v) return zeroser(vx, v+l-2);
371:
372: z = cgetg(l-i+3+v,t_SER);
373: z[1] = evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
374: for (j=i-1; j<l+v; j++) z[j-i+3]=lmul((GEN)x[j-ex],(GEN)y[j-ey]);
375: return z;
376: }
377:
378: /******************************************************************/
379: /** **/
380: /** PRECISION CHANGES **/
381: /** **/
382: /******************************************************************/
383:
384: GEN
385: gprec(GEN x, long l)
386: {
387: long tx=typ(x),lx=lg(x),i,pr;
388: GEN y;
389:
390: if (l<=0) err(talker,"precision<=0 in gprec");
391: switch(tx)
392: {
393: case t_REAL:
394: pr = (long) (l*pariK1+3); y=cgetr(pr); affrr(x,y); break;
395:
396: case t_PADIC:
1.2 ! noro 397: if (!signe(x[4]))
! 398: return padiczero((GEN)x[2], l+precp(x));
1.1 noro 399: y=cgetg(lx,tx); copyifstack(x[2], y[2]);
400: y[1]=x[1]; setprecp(y,l);
401: y[3]=lpuigs((GEN)x[2],l);
402: y[4]=lmodii((GEN)x[4],(GEN)y[3]);
403: break;
404:
405: case t_SER:
406: if (gcmp0(x)) return zeroser(varn(x), l);
407: y=cgetg(l+2,t_SER); y[1]=x[1]; l++; i=l;
408: if (l>=lx)
409: for ( ; i>=lx; i--) y[i]=zero;
410: for ( ; i>=2; i--) y[i]=lcopy((GEN)x[i]);
411: break;
412:
413: case t_POL:
414: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
415: for (i=2; i<lx; i++) y[i]=lprec((GEN)x[i],l);
416: break;
417:
418: case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
419: case t_VEC: case t_COL: case t_MAT:
420: y=cgetg(lx,tx);
421: for (i=1; i<lx; i++) y[i]=lprec((GEN)x[i],l);
422: break;
423: default: y=gcopy(x);
424: }
425: return y;
426: }
427:
428: /* internal: precision given in word length (including codewords) */
429: GEN
430: gprec_w(GEN x, long pr)
431: {
432: long tx=typ(x),lx,i;
433: GEN y;
434:
435: switch(tx)
436: {
437: case t_REAL:
438: y=cgetr(pr); affrr(x,y); break;
439:
440: case t_POL:
441: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
442: for (i=2; i<lx; i++) y[i]=(long)gprec_w((GEN)x[i],pr);
443: break;
444:
445: case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
446: case t_VEC: case t_COL: case t_MAT:
447: lx=lg(x); y=cgetg(lx,tx);
448: for (i=1; i<lx; i++) y[i]=(long)gprec_w((GEN)x[i],pr);
449: break;
450: default: y=gprec(x,pr);
451: }
452: return y;
453: }
454:
455: /*******************************************************************/
456: /** **/
457: /** RECIPROCAL POLYNOMIAL **/
458: /** **/
459: /*******************************************************************/
460:
461: GEN
462: polrecip(GEN x)
463: {
464: long lx=lgef(x),i,j;
465: GEN y;
466:
467: if (typ(x) != t_POL) err(typeer,"polrecip");
468: y=cgetg(lx,t_POL); y[1]=x[1];
469: for (i=2,j=lx-1; i<lx; i++,j--) y[i]=lcopy((GEN)x[j]);
470: return normalizepol_i(y,lx);
471: }
472:
473: /* as above. Internal (don't copy or normalize) */
474: GEN
475: polrecip_i(GEN x)
476: {
477: long lx=lgef(x),i,j;
478: GEN y;
479:
480: y=cgetg(lx,t_POL); y[1]=x[1];
481: for (i=2,j=lx-1; i<lx; i++,j--) y[i]=x[j];
482: return y;
483: }
484:
485: /*******************************************************************/
486: /** **/
487: /** BINOMIAL COEFFICIENTS **/
488: /** **/
489: /*******************************************************************/
490:
491: GEN
492: binome(GEN n, long k)
493: {
1.2 ! noro 494: long i;
! 495: gpmem_t av;
1.1 noro 496: GEN y;
497:
498: if (k <= 1)
499: {
500: if (is_noncalc_t(typ(n))) err(typeer,"binomial");
501: if (k < 0) return gzero;
502: if (k == 0) return gun;
503: return gcopy(n);
504: }
505: av = avma; y = n;
506: if (typ(n) == t_INT)
507: {
508: if (signe(n) > 0)
509: {
510: GEN z = subis(n,k);
511: if (cmpis(z,k) < 0) k = itos(z);
512: avma = av;
513: if (k <= 1)
514: {
515: if (k < 0) return gzero;
516: if (k == 0) return gun;
517: return gcopy(n);
518: }
519: }
520: for (i=2; i<=k; i++)
521: y = gdivgs(gmul(y,addis(n,i-1-k)), i);
522: }
523: else
524: {
525: for (i=2; i<=k; i++)
526: y = gdivgs(gmul(y,gaddgs(n,i-1-k)), i);
527: }
528: return gerepileupto(av, y);
529: }
530:
1.2 ! noro 531: /* Assume n >= 1, return bin, bin[k+1] = binomial(n, k) */
! 532: GEN
! 533: vecbinome(long n)
! 534: {
! 535: long d = (n + 1)/2, k;
! 536: GEN bin = cgetg(n+2, t_VEC), *C;
! 537: C = (GEN*)(bin + 1); /* C[k] = binomial(n, k) */
! 538: C[0] = gun;
! 539: for (k=1; k <= d; k++)
! 540: {
! 541: gpmem_t av = avma;
! 542: C[k] = gerepileuptoint(av, diviiexact(mulsi(n-k+1, C[k-1]), stoi(k)));
! 543: }
! 544: for ( ; k <= n; k++) C[k] = C[n - k];
! 545: return bin;
! 546: }
! 547:
1.1 noro 548: /********************************************************************/
549: /** **/
550: /** POLYNOMIAL INTERPOLATION **/
551: /** **/
552: /********************************************************************/
553: /* assume n > 1 */
554: GEN
555: polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy)
556: {
1.2 ! noro 557: long i, m, ns=0, tx=typ(x);
! 558: gpmem_t av = avma, tetpil;
1.1 noro 559: GEN den,ho,hp,w,y,c,d,dy;
560:
561: if (!xa)
562: {
563: xa = cgetg(n+1, t_VEC);
564: for (i=1; i<=n; i++) xa[i] = lstoi(i);
565: xa++;
566: }
567: if (is_scalar_t(tx) && tx != t_INTMOD && tx != t_PADIC && tx != t_POLMOD)
568: {
569: GEN dif = NULL, dift;
570: for (i=0; i<n; i++)
571: {
572: dift = gabs(gsub(x,(GEN)xa[i]), MEDDEFAULTPREC);
573: if (!dif || gcmp(dift,dif)<0) { ns=i; dif=dift; }
574: }
575: }
576: c=new_chunk(n);
577: d=new_chunk(n); for (i=0; i<n; i++) c[i] = d[i] = ya[i];
578: y=(GEN)d[ns--];
579: dy = NULL; tetpil = 0; /* gcc -Wall */
580: for (m=1; m<n; m++)
581: {
582: for (i=0; i<n-m; i++)
583: {
584: ho = gsub((GEN)xa[i],x);
585: hp = gsub((GEN)xa[i+m],x); den = gsub(ho,hp);
586: if (gcmp0(den)) err(talker,"two abcissas are equal in polint");
587: w=gsub((GEN)c[i+1],(GEN)d[i]); den = gdiv(w,den);
588: c[i]=lmul(ho,den);
589: d[i]=lmul(hp,den);
590: }
591: dy = (2*(ns+1) < n-m)? (GEN)c[ns+1]: (GEN)d[ns--];
592: tetpil=avma; y=gadd(y,dy);
593: }
594: if (!ptdy) y = gerepile(av,tetpil,y);
595: else
596: {
597: GEN *gptr[2];
598: *ptdy=gcopy(dy); gptr[0]=&y; gptr[1]=ptdy;
599: gerepilemanysp(av,tetpil,gptr,2);
600: }
601: return y;
602: }
603:
604: GEN
605: polint(GEN xa, GEN ya, GEN x, GEN *ptdy)
606: {
607: long tx=typ(xa), ty, lx=lg(xa);
1.2 ! noro 608:
1.1 noro 609: if (ya) ty = typ(ya); else { ya = xa; ty = tx; xa = NULL; }
610:
611: if (! is_vec_t(tx) || ! is_vec_t(ty))
612: err(talker,"not vectors in polinterpolate");
613: if (lx != lg(ya))
614: err(talker,"different lengths in polinterpolate");
615: if (lx <= 2)
616: {
617: if (lx == 1) err(talker,"no data in polinterpolate");
618: ya=gcopy((GEN)ya[1]); if (ptdy) *ptdy = ya;
619: return ya;
620: }
621: if (!x) x = polx[0];
622: return polint_i(xa? xa+1: xa,ya+1,x,lx-1,ptdy);
623: }
624:
625: /***********************************************************************/
626: /* */
627: /* SET OPERATIONS */
628: /* */
629: /***********************************************************************/
630:
631: static GEN
632: gtostr(GEN x)
633: {
634: char *s=GENtostr(x);
635: x = strtoGENstr(s,0); free(s); return x;
636: }
637:
638: GEN
639: gtoset(GEN x)
640: {
1.2 ! noro 641: gpmem_t av;
1.1 noro 642: long i,c,tx,lx;
643: GEN y;
644:
645: if (!x) return cgetg(1, t_VEC);
646: tx = typ(x); lx = lg(x);
647: if (!is_vec_t(tx))
648: {
649: if (tx != t_LIST)
650: { y=cgetg(2,t_VEC); y[1]=(long)gtostr(x); return y; }
651: lx = lgef(x)-1; x++;
652: }
653: if (lx==1) return cgetg(1,t_VEC);
654: av=avma; y=cgetg(lx,t_VEC);
655: for (i=1; i<lx; i++) y[i]=(long)gtostr((GEN)x[i]);
656: y = sort(y);
657: c=1;
658: for (i=2; i<lx; i++)
659: if (!gegal((GEN)y[i], (GEN)y[c])) y[++c] = y[i];
660: setlg(y,c+1); return gerepilecopy(av,y);
661: }
662:
663: long
664: setisset(GEN x)
665: {
666: long lx,i;
667:
668: if (typ(x)!=t_VEC) return 0;
669: lx=lg(x)-1; if (!lx) return 1;
670: for (i=1; i<lx; i++)
671: if (typ(x[i]) != t_STR || gcmp((GEN)x[i+1],(GEN)x[i])<=0) return 0;
672: return typ(x[i]) == t_STR;
673: }
674:
675: /* looks if y belongs to the set x and returns the index if yes, 0 if no */
676: long
1.2 ! noro 677: gen_search(GEN x, GEN y, int flag, int (*cmp)(GEN,GEN))
1.1 noro 678: {
1.2 ! noro 679: long lx,j,li,ri,fl, tx = typ(x);
1.1 noro 680:
681: if (tx==t_VEC) lx = lg(x);
682: else
683: {
684: if (tx!=t_LIST) err(talker,"not a set in setsearch");
685: lx=lgef(x)-1; x++;
686: }
687: if (lx==1) return flag? 1: 0;
688:
689: li=1; ri=lx-1;
690: do
691: {
1.2 ! noro 692: j = (ri+li)>>1; fl = cmp((GEN)x[j],y);
! 693: if (!fl) return flag? 0: j;
1.1 noro 694: if (fl<0) li=j+1; else ri=j-1;
695: } while (ri>=li);
1.2 ! noro 696: if (!flag) return 0;
1.1 noro 697: return (fl<0)? j+1: j;
698: }
699:
1.2 ! noro 700:
! 701: long
! 702: setsearch(GEN x, GEN y, long flag)
! 703: {
! 704: gpmem_t av = avma;
! 705: long res;
! 706: if (typ(y) != t_STR) y = gtostr(y);
! 707: res=gen_search(x,y,flag,gcmp);
! 708: avma=av;
! 709: return res;
! 710: }
! 711:
! 712: #if 0
! 713: GEN
! 714: gen_union(GEN x, GEN y, int (*cmp)(GEN,GEN))
! 715: {
! 716: if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
! 717:
! 718: }
! 719: #endif
! 720:
1.1 noro 721: GEN
722: setunion(GEN x, GEN y)
723: {
1.2 ! noro 724: gpmem_t av=avma, tetpil;
1.1 noro 725: GEN z;
726:
727: if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
728: z=concatsp(x,y); tetpil=avma; return gerepile(av,tetpil,gtoset(z));
729: }
730:
731: GEN
732: setintersect(GEN x, GEN y)
733: {
1.2 ! noro 734: long i, lx, c;
! 735: gpmem_t av=avma;
1.1 noro 736: GEN z;
737:
738: if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect");
739: lx=lg(x); z=cgetg(lx,t_VEC); c=1;
740: for (i=1; i<lx; i++)
741: if (setsearch(y, (GEN)x[i], 0)) z[c++] = x[i];
742: setlg(z,c); return gerepilecopy(av,z);
743: }
744:
745: GEN
1.2 ! noro 746: gen_setminus(GEN set1, GEN set2, int (*cmp)(GEN,GEN))
! 747: {
! 748: gpmem_t ltop=avma;
! 749: long find;
! 750: long i,j,k;
! 751: GEN diff=cgetg(lg(set1),t_VEC);
! 752: for(i=1,j=1,k=1; i < lg(set1); i++)
! 753: {
! 754: for(find=0; j < lg(set2); j++)
! 755: {
! 756: int s=cmp((GEN)set1[i],(GEN)set2[j]);
! 757: if (s<0) break ;
! 758: if (s>0) continue;
! 759: find=1;
! 760: }
! 761: if (!find)
! 762: diff[k++]=set1[i];
! 763: }
! 764: setlg(diff,k);
! 765: return gerepilecopy(ltop,diff);
! 766: }
! 767:
! 768: GEN
1.1 noro 769: setminus(GEN x, GEN y)
770: {
771: if (!setisset(x) || !setisset(y)) err(talker,"not a set in setminus");
1.2 ! noro 772: return gen_setminus(x,y,gcmp);
1.1 noro 773: }
774:
775: /***********************************************************************/
776: /* */
777: /* OPERATIONS ON DIRICHLET SERIES */
778: /* */
779: /***********************************************************************/
780:
781: /* Addition, subtraction and scalar multiplication of Dirichlet series
782: are done on the corresponding vectors */
783:
784: static long
785: dirval(GEN x)
786: {
787: long i=1,lx=lg(x);
788: while (i<lx && gcmp0((GEN)x[i])) i++;
789: return i;
790: }
791:
792: GEN
793: dirmul(GEN x, GEN y)
794: {
1.2 ! noro 795: gpmem_t av = avma, lim = stack_lim(av, 1);
1.1 noro 796: long lx,ly,lz,dx,dy,i,j,k;
797: GEN z,p1;
798:
799: if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
800: dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
801: if (ly-dy<lx-dx) { z=y; y=x; x=z; lz=ly; ly=lx; lx=lz; lz=dy; dy=dx; dx=lz; }
802: lz=min(lx*dy,ly*dx);
803: z=cgetg(lz,t_VEC); for (i=1; i<lz; i++) z[i]=zero;
804: for (j=dx; j<lx; j++)
805: {
806: p1=(GEN)x[j];
807: if (!gcmp0(p1))
808: {
809: if (gcmp1(p1))
810: for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=ladd((GEN)z[i],(GEN)y[k]);
811: else
812: {
813: if (gcmp_1(p1))
814: for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=lsub((GEN)z[i],(GEN)y[k]);
815: else
816: for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=ladd((GEN)z[i],gmul(p1,(GEN)y[k]));
817: }
818: }
819: if (low_stack(lim, stack_lim(av,1)))
820: {
821: if (DEBUGLEVEL) fprintferr("doubling stack in dirmul\n");
822: z = gerepilecopy(av,z);
823: }
824: }
825: return gerepilecopy(av,z);
826: }
827:
828: GEN
829: dirdiv(GEN x, GEN y)
830: {
1.2 ! noro 831: gpmem_t av = avma;
1.1 noro 832: long lx,ly,lz,dx,dy,i,j;
833: GEN z,p1;
834:
835: if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
836: dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
837: if (dy!=1) err(talker,"not an invertible dirseries in dirdiv");
838: lz=min(lx,ly*dx); p1=(GEN)y[1];
839: if (!gcmp1(p1)) { y=gdiv(y,p1); x=gdiv(x,p1); }
840: else x=gcopy(x);
841: z=cgetg(lz,t_VEC); for (i=1; i<dx; i++) z[i]=zero;
842: for (j=dx; j<lz; j++)
843: {
844: p1=(GEN)x[j]; z[j]=(long)p1;
845: if (!gcmp0(p1))
846: {
847: if (gcmp1(p1))
848: for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],(GEN)y[i/j]);
849: else
850: {
851: if (gcmp_1(p1))
852: for (i=j+j; i<lz; i+=j) x[i]=ladd((GEN)x[i],(GEN)y[i/j]);
853: else
854: for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],gmul(p1,(GEN)y[i/j]));
855: }
856: }
857: }
858: return gerepilecopy(av,z);
859: }
860:
861: /*************************************************************************/
862: /** **/
863: /** RANDOM **/
864: /** **/
865: /*************************************************************************/
866: static long pari_randseed = 1;
867:
868: /* BSD rand gives this: seed = 1103515245*seed + 12345 */
869: long
870: mymyrand(void)
871: {
872: #if BITS_IN_RANDOM == 64
873: pari_randseed = (1000000000000654397*pari_randseed + 12347) & ~HIGHBIT;
874: #else
875: pari_randseed = (1000276549*pari_randseed + 12347) & 0x7fffffff;
876: #endif
877: return pari_randseed;
878: }
879:
880: GEN muluu(ulong x, ulong y);
881:
882: static ulong
883: gp_rand(void)
884: {
885: #define GLUE2(hi, lo) (((hi) << BITS_IN_HALFULONG) | (lo))
886: #if !defined(LONG_IS_64BIT) || BITS_IN_RANDOM == 64
887: return GLUE2((mymyrand()>>12)&LOWMASK,
888: (mymyrand()>>12)&LOWMASK);
889: #else
890: #define GLUE4(hi1,hi2, lo1,lo2) GLUE2(((hi1)<<16)|(hi2), ((lo1)<<16)|(lo2))
891: # define LOWMASK2 0xffffUL
892: return GLUE4((mymyrand()>>12)&LOWMASK2,
893: (mymyrand()>>12)&LOWMASK2,
894: (mymyrand()>>12)&LOWMASK2,
895: (mymyrand()>>12)&LOWMASK2);
896: #endif
897: }
898:
899: GEN
900: genrand(GEN N)
901: {
902: long lx,i,nz;
903: GEN x, p1;
904:
905: if (!N) return stoi(mymyrand());
906: if (typ(N)!=t_INT || signe(N)<=0) err(talker,"invalid bound in random");
907:
908: lx = lgefint(N); x = new_chunk(lx);
909: nz = lx-1; while (!N[nz]) nz--; /* nz = index of last non-zero word */
910: for (i=2; i<lx; i++)
911: {
912: ulong n = N[i], r;
913: if (n == 0) r = 0;
914: else
1.2 ! noro 915: {
! 916: gpmem_t av = avma;
1.1 noro 917: if (i < nz) n++; /* allow for equality if we can go down later */
918: p1 = muluu(n, gp_rand()); /* < n * 2^32, so 0 <= first word < n */
919: r = (lgefint(p1)<=3)? 0: p1[2]; avma = av;
920: }
921: x[i] = r;
922: if (r < (ulong)N[i]) break;
923: }
924: for (i++; i<lx; i++) x[i] = gp_rand();
925: i=2; while (i<lx && !x[i]) i++;
926: i -= 2; x += i; lx -= i;
927: x[1] = evalsigne(lx>2) | evallgefint(lx);
928: x[0] = evaltyp(t_INT) | evallg(lx);
1.2 ! noro 929: avma = (gpmem_t)x; return x;
1.1 noro 930: }
931:
932: long
933: setrand(long seed) { return (pari_randseed = seed); }
934:
935: long
936: getrand(void) { return pari_randseed; }
937:
938: long
939: getstack(void) { return top-avma; }
940:
941: long
942: gettime(void) { return timer(); }
943:
944: /***********************************************************************/
945: /** **/
946: /** PERMUTATIONS **/
947: /** **/
948: /***********************************************************************/
949:
950: GEN
951: numtoperm(long n, GEN x)
952: {
1.2 ! noro 953: gpmem_t av;
1.1 noro 954: long i,a,r;
955: GEN v,w;
956:
957: if (n < 1) err(talker,"n too small (%ld) in numtoperm",n);
958: if (typ(x) != t_INT) err(arither1);
959: v = cgetg(n+1, t_VEC);
960: v[1]=1; av = avma;
961: if (signe(x) <= 0) x = modii(x, mpfact(n));
962: for (r=2; r<=n; r++)
963: {
964: x = dvmdis(x,r,&w); a = itos(w);
965: for (i=r; i>=a+2; i--) v[i] = v[i-1];
966: v[i]=r;
967: }
968: avma = av;
969: for (i=1; i<=n; i++) v[i] = lstoi(v[i]);
970: return v;
971: }
972:
973: GEN
974: permtonum(GEN x)
975: {
1.2 ! noro 976: long lx=lg(x)-1, n=lx, last, ind, tx = typ(x);
! 977: gpmem_t av=avma;
1.1 noro 978: GEN ary,res;
979:
980: if (!is_vec_t(tx)) err(talker,"not a vector in permtonum");
981: ary = cgetg(lx+1,t_VECSMALL);
982: for (ind=1; ind<=lx; ind++)
983: {
984: res = (GEN)*++x;
985: if (typ(res) != t_INT) err(typeer,"permtonum");
986: ary[ind] = itos(res);
987: }
988: ary++; res = gzero;
989: for (last=lx; last>0; last--)
990: {
991: lx--; ind = lx;
992: while (ind>0 && ary[ind] != last) ind--;
993: res = addis(mulis(res,last), ind);
994: while (ind++<lx) ary[ind-1] = ary[ind];
995: }
996: if (!signe(res)) res = mpfact(n);
997: return gerepileuptoint(av, res);
998: }
999:
1000: /********************************************************************/
1001: /** **/
1002: /** MODREVERSE **/
1003: /** **/
1004: /********************************************************************/
1.2 ! noro 1005: /* evaluate f(x) mod T */
! 1006: GEN
! 1007: RX_RXQ_compo(GEN f, GEN x, GEN T)
! 1008: {
! 1009: gpmem_t av = avma;
! 1010: long l;
! 1011: GEN y;
! 1012:
! 1013: if (typ(f) != t_POL) return gcopy(f);
! 1014: l = lgef(f)-1; y = (GEN)f[l];
! 1015: for (l--; l>=2; l--)
! 1016: y = gres(gadd(gmul(y,x), (GEN)f[l]), T);
! 1017: return gerepileupto(av, y);
! 1018: }
! 1019:
! 1020: /* return (1,...a^l) mod T. Not memory clean */
! 1021: GEN
! 1022: RXQ_powers(GEN a, GEN T, long l)
! 1023: {
! 1024: long i;
! 1025: GEN v;
! 1026:
! 1027: if (typ(a) != t_POL) err(typeer,"RXQ_powers");
! 1028: l += 2;
! 1029: v = cgetg(l,t_VEC);
! 1030: v[1] = un; if (l == 2) return v;
! 1031:
! 1032: if (degpol(a) >= degpol(T)) a = gres(a, T);
! 1033: v[2] = (long)a;
! 1034: for (i=3; i<l; i++) v[i] = lres(gmul((GEN)v[i-1], a), T);
! 1035: return v;
! 1036: }
! 1037:
! 1038: /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
! 1039: GEN
! 1040: modreverse_i(GEN a, GEN T)
! 1041: {
! 1042: gpmem_t av = avma;
! 1043: long n = degpol(T);
! 1044: GEN y;
! 1045:
! 1046: if (n <= 0) return gcopy(a);
! 1047: if (n == 1)
! 1048: return gerepileupto(av, gneg(gdiv((GEN)T[2], (GEN)T[3])));
! 1049: if (gcmp0(a) || typ(a) != t_POL) err(talker,"reverse polmod does not exist");
! 1050:
! 1051: y = vecpol_to_mat(RXQ_powers(a,T,n-1), n);
! 1052: y = gauss(y, _ei(n, 2));
! 1053: return gerepilecopy(av, vec_to_pol(y, varn(T)));
! 1054: }
1.1 noro 1055:
1056: GEN
1057: polymodrecip(GEN x)
1058: {
1.2 ! noro 1059: long v, n;
! 1060: GEN T, a, y;
1.1 noro 1061:
1.2 ! noro 1062: if (typ(x)!=t_POLMOD) err(talker,"not a polmod in modreverse");
! 1063: T = (GEN)x[1];
! 1064: a = (GEN)x[2];
! 1065: n = degpol(T); if (n <= 0) return gcopy(x);
! 1066: v = varn(T);
! 1067: y = cgetg(3,t_POLMOD);
! 1068: y[1] = (n==1)? lsub(polx[v], a): (long)caract2(T, a, v);
! 1069: y[2] = (long)modreverse_i(a, T); return y;
1.1 noro 1070: }
1071:
1072: /********************************************************************/
1073: /** **/
1074: /** HEAPSORT **/
1075: /** **/
1076: /********************************************************************/
1077: static GEN vcmp_k;
1078: static int vcmp_lk;
1079: static int (*vcmp_cmp)(GEN,GEN);
1080:
1081: int
1082: pari_compare_int(int *a,int *b)
1083: {
1084: return *a - *b;
1085: }
1086:
1087: int
1088: pari_compare_long(long *a,long *b)
1089: {
1090: return *a - *b;
1091: }
1092:
1093: static int
1094: veccmp(GEN x, GEN y)
1095: {
1096: int i,s;
1097:
1098: for (i=1; i<vcmp_lk; i++)
1099: {
1100: s = vcmp_cmp((GEN) x[vcmp_k[i]], (GEN) y[vcmp_k[i]]);
1101: if (s) return s;
1102: }
1103: return 0;
1104: }
1105:
1106: static int
1107: longcmp(GEN x, GEN y)
1108: {
1109: return ((long)x > (long)y)? 1: ((x == y)? 0: -1);
1110: }
1111:
1112: /* Sort x = vector of elts, using cmp to compare them.
1113: * flag & cmp_IND: indirect sort: return permutation that would sort x
1114: * flag & cmp_C : as cmp_IND, but return permutation as vector of C-longs
1115: */
1116: GEN
1117: gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN))
1118: {
1119: long i,j,indxt,ir,l,tx=typ(x),lx=lg(x);
1120: GEN q,y,indx;
1121:
1122: if (!is_matvec_t(tx) && tx != t_VECSMALL) err(typeer,"gen_sort");
1123: if (flag & cmp_C) tx = t_VECSMALL;
1124: else if (flag & cmp_IND) tx = t_VEC;
1125: y = cgetg(lx, tx);
1126: if (lx==1) return y;
1127: if (lx==2)
1128: {
1129: if (flag & cmp_C)
1130: y[1] = 1;
1131: else if (flag & cmp_IND)
1132: y[1] = un;
1133: else
1134: y[1] = lcopy((GEN)x[1]);
1135: return y;
1136: }
1137: if (!cmp) cmp = &longcmp;
1138: indx = (GEN) gpmalloc(lx*sizeof(long));
1.2 ! noro 1139: for (j=1; j<lx; j++) indx[j] = j;
1.1 noro 1140:
1.2 ! noro 1141: ir = lx-1; l = (ir>>1)+1;
1.1 noro 1142: for(;;)
1143: {
1.2 ! noro 1144: if (l > 1) indxt = indx[--l];
1.1 noro 1145: else
1146: {
1.2 ! noro 1147: indxt = indx[ir]; indx[ir] = indx[1];
! 1148: if (--ir == 1) break; /* DONE */
1.1 noro 1149: }
1.2 ! noro 1150: q = (GEN)x[indxt]; i = l;
! 1151: for (j=i<<1; j<=ir; j<<=1)
! 1152: {
! 1153: if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) < 0) j++;
! 1154: if (cmp(q,(GEN)x[indx[j]]) >= 0) break;
1.1 noro 1155:
1.2 ! noro 1156: indx[i] = indx[j]; i = j;
! 1157: }
! 1158: indx[i] = indxt;
! 1159: }
! 1160: indx[1] = indxt;
1.1 noro 1161:
1.2 ! noro 1162: if (flag & cmp_REV)
! 1163: { /* reverse order */
! 1164: GEN indy = (GEN) gpmalloc(lx*sizeof(long));
! 1165: for (j=1; j<lx; j++) indy[j] = indx[lx-j];
! 1166: free(indx); indx = indy;
! 1167: }
! 1168: if (flag & cmp_C)
! 1169: for (i=1; i<lx; i++) y[i] = indx[i];
! 1170: else if (flag & cmp_IND)
! 1171: for (i=1; i<lx; i++) y[i] = lstoi(indx[i]);
! 1172: else
! 1173: for (i=1; i<lx; i++) y[i] = lcopy((GEN)x[indx[i]]);
! 1174: free(indx); return y;
1.1 noro 1175: }
1176:
1177: #define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp)
1178:
1179: GEN
1180: gen_vecsort(GEN x, GEN k, long flag)
1181: {
1182: long i,j,l,t, lx = lg(x), tmp[2];
1183:
1184: if (lx<=2) return gen_sort(x,flag,sort_fun(flag));
1185: t = typ(k); vcmp_cmp = sort_fun(flag);
1186: if (t==t_INT)
1187: {
1188: tmp[1] = (long)k; k = tmp;
1189: vcmp_lk = 2;
1190: }
1191: else
1192: {
1193: if (! is_vec_t(t)) err(talker,"incorrect lextype in vecsort");
1194: vcmp_lk = lg(k);
1195: }
1196: l = 0;
1197: vcmp_k = (GEN)gpmalloc(vcmp_lk * sizeof(long));
1198: for (i=1; i<vcmp_lk; i++)
1199: {
1200: j = itos((GEN)k[i]);
1201: if (j<=0) err(talker,"negative index in vecsort");
1202: vcmp_k[i]=j; if (j>l) l=j;
1203: }
1204: t = typ(x);
1205: if (! is_matvec_t(t)) err(typeer,"vecsort");
1206: for (j=1; j<lx; j++)
1207: {
1208: t = typ(x[j]);
1209: if (! is_vec_t(t)) err(typeer,"vecsort");
1210: if (lg((GEN)x[j]) <= l) err(talker,"index too large in vecsort");
1211: }
1212: x = gen_sort(x, flag, veccmp);
1213: free(vcmp_k); return x;
1214: }
1215:
1216: GEN
1217: vecsort0(GEN x, GEN k, long flag)
1218: {
1219: if (flag < 0 || flag >= cmp_C) err(flagerr,"vecsort");
1220: return k? gen_vecsort(x,k,flag): gen_sort(x,flag, sort_fun(flag));
1221: }
1222:
1223: GEN
1224: vecsort(GEN x, GEN k)
1225: {
1226: return gen_vecsort(x,k, 0);
1227: }
1228:
1229: GEN
1230: sindexsort(GEN x)
1231: {
1232: return gen_sort(x, cmp_IND | cmp_C, gcmp);
1233: }
1234:
1235: GEN
1236: sindexlexsort(GEN x)
1237: {
1238: return gen_sort(x, cmp_IND | cmp_C, lexcmp);
1239: }
1240:
1241: GEN
1242: indexsort(GEN x)
1243: {
1244: return gen_sort(x, cmp_IND, gcmp);
1245: }
1246:
1247: GEN
1248: indexlexsort(GEN x)
1249: {
1250: return gen_sort(x, cmp_IND, lexcmp);
1251: }
1252:
1253: GEN
1254: sort(GEN x)
1255: {
1256: return gen_sort(x, 0, gcmp);
1257: }
1258:
1259: GEN
1260: lexsort(GEN x)
1261: {
1262: return gen_sort(x, 0, lexcmp);
1263: }
1264:
1265: /* index of x in table T, 0 otherwise */
1266: long
1267: tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
1268: {
1269: long l=1,u=lg(T)-1,i,s;
1270:
1271: while (u>=l)
1272: {
1273: i = (l+u)>>1; s = cmp(x,(GEN)T[i]);
1274: if (!s) return i;
1275: if (s<0) u=i-1; else l=i+1;
1276: }
1277: return 0;
1278: }
1279:
1280: /* assume lg(x) = lg(y), x,y in Z^n */
1281: int
1282: cmp_vecint(GEN x, GEN y)
1283: {
1284: long fl,i, lx = lg(x);
1285: for (i=1; i<lx; i++)
1286: if (( fl = cmpii((GEN)x[i], (GEN)y[i]) )) return fl;
1287: return 0;
1288: }
1289:
1290: /* assume x and y come from the same primedec call (uniformizer unique) */
1291: int
1292: cmp_prime_over_p(GEN x, GEN y)
1293: {
1294: int k = mael(x,4,2) - mael(y,4,2); /* diff. between residue degree */
1295: return k? ((k > 0)? 1: -1)
1296: : cmp_vecint((GEN)x[2], (GEN)y[2]);
1297: }
1298:
1299: int
1300: cmp_prime_ideal(GEN x, GEN y)
1301: {
1302: int k = cmpii((GEN)x[1], (GEN)y[1]);
1303: return k? k: cmp_prime_over_p(x,y);
1304: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bibli2.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