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