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