Annotation of OpenXM_contrib/pari-2.2/src/basemath/gen2.c, Revision 1.1
1.1 ! noro 1: /* $Id: gen2.c,v 1.33 2001/10/01 12:11:30 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: /** GENERIC OPERATIONS **/
! 19: /** (second part) **/
! 20: /** **/
! 21: /********************************************************************/
! 22: #include "pari.h"
! 23:
! 24: /*******************************************************************/
! 25: /* */
! 26: /* OPERATIONS BY VALUE */
! 27: /* f is a pointer to the function called. */
! 28: /* result is gaffect-ed to last parameter */
! 29: /* */
! 30: /*******************************************************************/
! 31:
! 32: void
! 33: gop1z(GEN (*f)(GEN), GEN x, GEN y)
! 34: {
! 35: long av=avma; gaffect(f(x),y); avma=av;
! 36: }
! 37:
! 38: void
! 39: gop2z(GEN (*f)(GEN, GEN), GEN x, GEN y, GEN z)
! 40: {
! 41: long av=avma; gaffect(f(x,y),z); avma=av;
! 42: }
! 43:
! 44: void
! 45: gops2gsz(GEN (*f)(GEN, long), GEN x, long s, GEN z)
! 46: {
! 47: long av=avma; gaffect(f(x,s),z); avma=av;
! 48: }
! 49:
! 50: void
! 51: gops2sgz(GEN (*f)(long, GEN), long s, GEN y, GEN z)
! 52: {
! 53: long av=avma; gaffect(f(s,y),z); avma=av;
! 54: }
! 55:
! 56: void
! 57: gops2ssz(GEN (*f)(long, long), long s, long y, GEN z)
! 58: {
! 59: long av=avma; gaffect(f(s,y),z); avma=av;
! 60: }
! 61:
! 62: /*******************************************************************/
! 63: /* */
! 64: /* OPERATIONS USING SMALL INTEGERS */
! 65: /* */
! 66: /*******************************************************************/
! 67:
! 68: /* small int prototype */
! 69: static long court_p[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
! 70:
! 71: GEN
! 72: gopsg2(GEN (*f)(GEN, GEN), long s, GEN y)
! 73: {
! 74: affsi(s,court_p); return f(court_p,y);
! 75: }
! 76:
! 77: GEN
! 78: gopgs2(GEN (*f)(GEN, GEN), GEN y, long s)
! 79: {
! 80: affsi(s,court_p); return f(y,court_p);
! 81: }
! 82:
! 83: long
! 84: opgs2(int (*f)(GEN, GEN), GEN y, long s)
! 85: {
! 86: affsi(s,court_p); return f(y,court_p);
! 87: }
! 88:
! 89: void
! 90: gopsg2z(GEN (*f)(GEN, GEN), long s, GEN y, GEN z)
! 91: {
! 92: long av=avma;
! 93: affsi(s,court_p); gaffect(f(court_p,y),z); avma=av;
! 94: }
! 95:
! 96: void
! 97: gopgs2z(GEN (*f)(GEN, GEN), GEN y, long s, GEN z)
! 98: {
! 99: long av=avma;
! 100: affsi(s,court_p); gaffect(f(y,court_p),z); avma=av;
! 101: }
! 102:
! 103: /*******************************************************************/
! 104: /* */
! 105: /* CREATION OF A P-ADIC GEN */
! 106: /* */
! 107: /*******************************************************************/
! 108:
! 109: /* y[4] not filled */
! 110: static GEN
! 111: cgetp2(GEN x, long v)
! 112: {
! 113: GEN y = cgetg(5,t_PADIC);
! 114: y[1] = evalprecp(precp(x)) | evalvalp(v);
! 115: icopyifstack(x[2], y[2]);
! 116: y[3] = licopy((GEN)x[3]); return y;
! 117: }
! 118:
! 119: GEN
! 120: cgetp(GEN x)
! 121: {
! 122: GEN y = cgetp2(x, 0);
! 123: y[4] = lgeti(lgefint(x[3])); return y;
! 124: }
! 125:
! 126: GEN
! 127: pureimag(GEN x)
! 128: {
! 129: GEN y = cgetg(3,t_COMPLEX);
! 130: y[1] = zero;
! 131: y[2] = (long)x; return y;
! 132: }
! 133:
! 134: /*******************************************************************/
! 135: /* */
! 136: /* SIZES */
! 137: /* */
! 138: /*******************************************************************/
! 139:
! 140: long
! 141: glength(GEN x)
! 142: {
! 143: switch(typ(x))
! 144: {
! 145: case t_INT: return lgefint(x)-2;
! 146: case t_POL: case t_LIST: return lgef(x)-2;
! 147: case t_REAL: return signe(x)? lg(x)-2: 0;
! 148: case t_STR: return strlen(GSTR(x));
! 149: }
! 150: return lg(x)-lontyp[typ(x)];
! 151: }
! 152:
! 153: GEN
! 154: matsize(GEN x)
! 155: {
! 156: GEN y=cgetg(3,t_VEC);
! 157:
! 158: switch(typ(x))
! 159: {
! 160: case t_VEC:
! 161: y[1]=un; y[2]=lstoi(lg(x)-1); break;
! 162: case t_COL:
! 163: y[1]=lstoi(lg(x)-1); y[2]=un; break;
! 164: case t_MAT:
! 165: y[2]=lstoi(lg(x)-1);
! 166: y[1]=(lg(x)==1)? zero: lstoi(lg(x[1])-1); break;
! 167: default: err(typeer,"matsize");
! 168: }
! 169: return y;
! 170: }
! 171:
! 172: long
! 173: taille(GEN x)
! 174: {
! 175: long i,n,lx, tx = typ(x);
! 176: n = lg(x);
! 177: if (is_recursive_t(tx))
! 178: {
! 179: lx = (tx==t_POL || tx==t_LIST)? lgef(x): n;
! 180: for (i=lontyp[tx]; i<lx; i++) n += taille((GEN)x[i]);
! 181: }
! 182: return n;
! 183: }
! 184:
! 185: long
! 186: taille2(GEN x) { return taille(x)<<TWOPOTBYTES_IN_LONG; }
! 187:
! 188: /*******************************************************************/
! 189: /* */
! 190: /* GREFFE */
! 191: /* Greffe d'une serie sur un polynome */
! 192: /* */
! 193: /*******************************************************************/
! 194:
! 195: GEN
! 196: greffe(GEN x, long l, long use_stack)
! 197: {
! 198: long i,e,k,vx;
! 199: GEN y;
! 200:
! 201: if (typ(x)!=t_POL) err(notpoler,"greffe");
! 202: if (use_stack) y = cgetg(l,t_SER);
! 203: else
! 204: {
! 205: y = (GEN) gpmalloc(l*sizeof(long));
! 206: y[0] = evaltyp(t_SER) | evallg(l);
! 207: }
! 208: if (gcmp0(x))
! 209: {
! 210: y[1]=evalvalp(l-2) | evalvarn(varn(x));
! 211: i=2; while(i<l) { y[i]=x[2]; i++; }
! 212: return y;
! 213: }
! 214: vx=varn(x); e=gval(x,vx);
! 215: y[1]=evalsigne(1) | evalvalp(e) | evalvarn(vx);
! 216: k=lgef(x)-e-1; i=l-1;
! 217: if (k<l)
! 218: while (i>k) { y[i]=zero; i--; }
! 219: while (i>=2) { y[i]=x[i+e]; i--; }
! 220: return y;
! 221: }
! 222:
! 223: /*******************************************************************/
! 224: /* */
! 225: /* CONVERSION GEN --> long */
! 226: /* */
! 227: /*******************************************************************/
! 228:
! 229: long
! 230: gtolong(GEN x)
! 231: {
! 232: long y,tx=typ(x),av=avma;
! 233:
! 234: switch(tx)
! 235: {
! 236: case t_INT:
! 237: return itos(x);
! 238: case t_REAL: case t_FRAC: case t_FRACN:
! 239: y=itos(ground(x)); avma=av; return y;
! 240: case t_COMPLEX:
! 241: if (gcmp0((GEN)x[2])) return gtolong((GEN)x[1]); break;
! 242: case t_QUAD:
! 243: if (gcmp0((GEN)x[3])) return gtolong((GEN)x[2]); break;
! 244: }
! 245: err(typeer,"gtolong");
! 246: return 0; /* not reached */
! 247: }
! 248:
! 249: /*******************************************************************/
! 250: /* */
! 251: /* COMPARISONS */
! 252: /* */
! 253: /*******************************************************************/
! 254:
! 255: /* returns 1 whenever x = 0, and 0 otherwise */
! 256: int
! 257: gcmp0(GEN x)
! 258: {
! 259: switch(typ(x))
! 260: {
! 261: case t_INT: case t_REAL: case t_POL: case t_SER:
! 262: return !signe(x);
! 263:
! 264: case t_INTMOD: case t_POLMOD:
! 265: return gcmp0((GEN)x[2]);
! 266:
! 267: case t_FRAC: case t_FRACN:
! 268: return !signe(x[1]);
! 269:
! 270: case t_COMPLEX:
! 271: /* is 0 iff norm(x) would be 0 (can happen with Re(x) and Im(x) != 0
! 272: * only if Re(x) and Im(x) are of type t_REAL). See mp.c:addrr().
! 273: */
! 274: if (gcmp0((GEN)x[1]))
! 275: {
! 276: if (gcmp0((GEN)x[2])) return 1;
! 277: if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
! 278: return (expo(x[1])>expo(x[2]));
! 279: }
! 280: if (gcmp0((GEN)x[2]))
! 281: {
! 282: if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
! 283: return (expo(x[2])>expo(x[1]));
! 284: }
! 285: return 0;
! 286:
! 287: case t_PADIC:
! 288: return !signe(x[4]);
! 289:
! 290: case t_QUAD:
! 291: return gcmp0((GEN)x[2]) && gcmp0((GEN)x[3]);
! 292:
! 293: case t_RFRAC: case t_RFRACN:
! 294: return gcmp0((GEN)x[1]);
! 295:
! 296: case t_VEC: case t_COL: case t_MAT:
! 297: {
! 298: long i;
! 299: for (i=lg(x)-1; i; i--)
! 300: if (!gcmp0((GEN)x[i])) return 0;
! 301: return 1;
! 302: }
! 303: }
! 304: return 0;
! 305: }
! 306:
! 307: /* returns 1 whenever x = 1, 0 otherwise */
! 308: int
! 309: gcmp1(GEN x)
! 310: {
! 311: switch(typ(x))
! 312: {
! 313: case t_INT:
! 314: return is_pm1(x) && signe(x)==1;
! 315:
! 316: case t_REAL:
! 317: if (signe(x) > 0 && expo(x)==0 && (ulong)x[2]==HIGHBIT)
! 318: {
! 319: long i,lx = lg(x);
! 320: for (i=3; i<lx; i++)
! 321: if (x[i]) return 0;
! 322: return 1;
! 323: }
! 324: return 0;
! 325:
! 326: case t_INTMOD: case t_POLMOD:
! 327: return gcmp1((GEN)x[2]);
! 328:
! 329: case t_FRAC:
! 330: return gcmp1((GEN)x[1]) && gcmp1((GEN)x[2]);
! 331:
! 332: case t_FRACN:
! 333: return egalii((GEN)x[1],(GEN)x[2]);
! 334:
! 335: case t_COMPLEX:
! 336: return gcmp1((GEN)x[1]) && gcmp0((GEN)x[2]);
! 337:
! 338: case t_PADIC:
! 339: return !valp(x) && gcmp1((GEN)x[4]);
! 340:
! 341: case t_QUAD:
! 342: return gcmp1((GEN)x[2]) && gcmp0((GEN)x[3]);
! 343:
! 344: case t_POL:
! 345: return lgef(x)==3 && gcmp1((GEN)x[2]);
! 346: }
! 347: return 0;
! 348: }
! 349:
! 350: /* returns 1 whenever the x = -1, 0 otherwise */
! 351: int
! 352: gcmp_1(GEN x)
! 353: {
! 354: long l,y;
! 355: GEN p1;
! 356:
! 357: switch(typ(x))
! 358: {
! 359: case t_INT:
! 360: return is_pm1(x) && signe(x)== -1;
! 361:
! 362: case t_REAL:
! 363: if (signe(x) < 0 && expo(x)==0 && (ulong)x[2]==HIGHBIT)
! 364: {
! 365: long i,lx = lg(x);
! 366: for (i=3; i<lx; i++)
! 367: if (x[i]) return 0;
! 368: return 1;
! 369: }
! 370: return 0;
! 371:
! 372: case t_INTMOD:
! 373: l=avma; y=egalii(addsi(1,(GEN)x[2]), (GEN)x[1]); avma=l; return y;
! 374:
! 375: case t_FRAC: case t_FRACN:
! 376: l = signe(x[1]);
! 377: return l && l == -signe(x[2]) && !absi_cmp((GEN)x[1],(GEN)x[2]);
! 378:
! 379: case t_COMPLEX:
! 380: return gcmp_1((GEN)x[1]) && gcmp0((GEN)x[2]);
! 381:
! 382: case t_QUAD:
! 383: return gcmp_1((GEN)x[2]) && gcmp0((GEN)x[3]);
! 384:
! 385: case t_PADIC:
! 386: l=avma; y=gegal(addsi(1,(GEN)x[4]), (GEN)x[3]); avma=l; return y;
! 387:
! 388: case t_POLMOD:
! 389: l=avma; p1=gadd(gun,(GEN)x[2]);
! 390: y = signe(p1) && !gegal(p1,(GEN)x[1]); avma=l; return !y;
! 391:
! 392: case t_POL:
! 393: return lgef(x)==3 && gcmp_1((GEN)x[2]);
! 394: }
! 395: return 0;
! 396: }
! 397:
! 398: /* returns the sign of x - y when it makes sense. 0 otherwise */
! 399: int
! 400: gcmp(GEN x, GEN y)
! 401: {
! 402: long tx,ty,f,av;
! 403:
! 404: tx=typ(x); ty=typ(y);
! 405: if (is_intreal_t(tx))
! 406: { if (is_intreal_t(ty)) return mpcmp(x,y); }
! 407: else
! 408: {
! 409: if (tx==t_STR)
! 410: {
! 411: if (ty != t_STR) return 1;
! 412: return strcmp(GSTR(x),GSTR(y));
! 413: }
! 414: if (!is_frac_t(tx)) err(typeer,"comparison");
! 415: }
! 416: if (ty == t_STR) return -1;
! 417: if (!is_intreal_t(ty) && !is_frac_t(ty)) err(typeer,"comparison");
! 418: av=avma; y=gneg_i(y); f=gsigne(gadd(x,y)); avma=av; return f;
! 419: }
! 420:
! 421: /* as gcmp for vector/matrices, using lexicographic ordering on components */
! 422: int
! 423: lexcmp(GEN x, GEN y)
! 424: {
! 425: const long tx=typ(x), ty=typ(y);
! 426: long lx,ly,l,fl,i;
! 427:
! 428: ly=lg(y);
! 429: if (!is_matvec_t(tx))
! 430: {
! 431: if (!is_matvec_t(ty)) return gcmp(x,y);
! 432: if (ly==1) return 1;
! 433: fl = lexcmp(x,(GEN)y[1]);
! 434: if (fl) return fl;
! 435: return (ly>2)? -1:0;
! 436: }
! 437:
! 438: lx=lg(x);
! 439: if (!is_matvec_t(ty))
! 440: {
! 441: if (lx==1) return -1;
! 442: fl = lexcmp(y,(GEN)x[1]);
! 443: if (fl) return -fl;
! 444: return (lx>2)? 1:0;
! 445: }
! 446:
! 447: /* x and y are matvec_t */
! 448: if (ly==1) return (lx==1)?0:1;
! 449: if (lx==1) return -1;
! 450: if (ty==t_MAT)
! 451: {
! 452: if (tx != t_MAT)
! 453: {
! 454: fl = lexcmp(x,(GEN)y[1]);
! 455: if (fl) return fl;
! 456: return (ly>2)?-1:0;
! 457: }
! 458: }
! 459: else if (tx==t_MAT)
! 460: {
! 461: fl = lexcmp(y,(GEN)x[1]);
! 462: if (fl) return -fl;
! 463: return (ly>2)?1:0;
! 464: }
! 465:
! 466: /* tx = ty = t_MAT, or x and y are both vect_t */
! 467: l=min(lx,ly);
! 468: for (i=1; i<l; i++)
! 469: {
! 470: fl = lexcmp((GEN)x[i],(GEN)y[i]);
! 471: if (fl) return fl;
! 472: }
! 473: return (ly != lx)? -1 : 0;
! 474: }
! 475:
! 476: /*****************************************************************/
! 477: /* */
! 478: /* EQUALITY */
! 479: /* returns 1 if x == y, 0 otherwise */
! 480: /* */
! 481: /*****************************************************************/
! 482:
! 483: /* x,y t_POL */
! 484: static int
! 485: polegal(GEN x, GEN y)
! 486: {
! 487: long i, lx = lgef(x);
! 488:
! 489: if (x[1] != y[1] && (lgef(y) != lx || lx > 3)) return 0;
! 490: for (i = 2; i < lx; i++)
! 491: if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
! 492: return 1;
! 493: }
! 494:
! 495: #define MASK(x) (((ulong)(x)) & (TYPBITS | LGBITS))
! 496: /* typ(x) = typ(y) = t_VEC/COL/MAT */
! 497: static int
! 498: vecegal(GEN x, GEN y)
! 499: {
! 500: long i;
! 501:
! 502: if (MASK(x[0]) != MASK(y[0])) return 0;
! 503:
! 504: i = lg(x)-1;
! 505: if (typ(x) != t_MAT)
! 506: {
! 507: for ( ; i; i--)
! 508: if (! gegal((GEN)x[i],(GEN)y[i]) ) return 0;
! 509: return 1;
! 510: }
! 511: for ( ; i; i--)
! 512: if (! vecegal((GEN)x[i],(GEN)y[i]) ) return 0;
! 513: return 1;
! 514: }
! 515:
! 516: int
! 517: gegal(GEN x, GEN y)
! 518: {
! 519: ulong av,tx;
! 520: long i;
! 521:
! 522: if (x == y) return 1;
! 523: tx = typ(x);
! 524: if (tx==typ(y))
! 525: switch(tx)
! 526: {
! 527: case t_INT:
! 528: return egalii(x,y);
! 529:
! 530: case t_POL:
! 531: return polegal(x,y);
! 532:
! 533: case t_COMPLEX:
! 534: return gegal((GEN)x[1],(GEN)y[1]) && gegal((GEN)x[2],(GEN)y[2]);
! 535:
! 536: case t_INTMOD: case t_POLMOD:
! 537: return gegal((GEN)x[2],(GEN)y[2])
! 538: && (x[1]==y[1] || gegal((GEN)x[1],(GEN)y[1]));
! 539:
! 540: case t_QFR:
! 541: if (!gegal((GEN)x[4],(GEN)y[4])) return 0; /* fall through */
! 542: case t_QUAD: case t_QFI:
! 543: return gegal((GEN)x[1],(GEN)y[1])
! 544: && gegal((GEN)x[2],(GEN)y[2])
! 545: && gegal((GEN)x[3],(GEN)y[3]);
! 546:
! 547: case t_FRAC:
! 548: return gegal((GEN)x[1], (GEN)y[1])
! 549: && gegal((GEN)x[2], (GEN)y[2]);
! 550:
! 551: case t_FRACN: case t_RFRAC: case t_RFRACN:
! 552: av=avma; i=gegal(gmul((GEN)x[1],(GEN)y[2]),gmul((GEN)x[2],(GEN)y[1]));
! 553: avma=av; return i;
! 554:
! 555: case t_STR:
! 556: return !strcmp(GSTR(x),GSTR(y));
! 557:
! 558: case t_VEC: case t_COL: case t_MAT:
! 559: return vecegal(x,y);
! 560: case t_VECSMALL:
! 561: if (MASK(x[0]) != MASK(y[0])) return 0;
! 562: for (i = lg(x)-1; i; i--)
! 563: if (x[i] != y[i]) return 0;
! 564: return 1;
! 565: }
! 566: {
! 567: jmp_buf env;
! 568: VOLATILE long av = avma;
! 569: void *c;
! 570: if (setjmp(env)) i = 0;
! 571: else
! 572: {
! 573: c = err_catch(-1, env, NULL);
! 574: i = gcmp0(gadd(x, gneg_i(y)));
! 575: }
! 576: err_leave(&c);
! 577: avma = av;
! 578: }
! 579: return i;
! 580: }
! 581: #undef MASK
! 582:
! 583: /*******************************************************************/
! 584: /* */
! 585: /* VALUATION */
! 586: /* p is either an int or a polynomial. */
! 587: /* returns the largest exponent of p dividing x when this makes */
! 588: /* sense : error for types real, integermod and polymod if p does */
! 589: /* not divide the modulus, q-adic if q!=p. */
! 590: /* */
! 591: /*******************************************************************/
! 592:
! 593: static long
! 594: minval(GEN x, GEN p, long first, long lx)
! 595: {
! 596: long i,k, val = VERYBIGINT;
! 597: for (i=first; i<lx; i++)
! 598: {
! 599: k = ggval((GEN)x[i],p);
! 600: if (k < val) val = k;
! 601: }
! 602: return val;
! 603: }
! 604:
! 605: long
! 606: ggval(GEN x, GEN p)
! 607: {
! 608: long tx=typ(x), tp=typ(p), av, limit,vx,v,i,val;
! 609: GEN p1,p2;
! 610:
! 611: if (isexactzero(x)) return VERYBIGINT;
! 612: if (is_const_t(tx) && tp==t_POL) return 0;
! 613: switch(tx)
! 614: {
! 615: case t_INT:
! 616: if (tp != t_INT) break;
! 617: av=avma; val = pvaluation(x,p, &p1);
! 618: avma=av; return val;
! 619:
! 620: case t_INTMOD:
! 621: av=avma;
! 622: p1=cgeti(lgefint(x[1]));
! 623: p2=cgeti(lgefint(x[2]));
! 624: if (tp!=t_INT || !mpdivis((GEN)x[1],p,p1)) break;
! 625: if (!mpdivis((GEN)x[2],p,p2)) { avma=av; return 0; }
! 626: val=1; while (mpdivis(p1,p,p1) && mpdivis(p2,p,p2)) val++;
! 627: avma=av; return val;
! 628:
! 629: case t_PADIC:
! 630: if (tp!=t_INT || !gegal(p,(GEN)x[2])) break;
! 631: return valp(x);
! 632:
! 633: case t_POLMOD:
! 634: if (tp==t_INT) return ggval((GEN)x[2],p);
! 635: av = avma;
! 636: if (tp!=t_POL || !poldivis((GEN)x[1],p,&p1)) break;
! 637: if (!poldivis((GEN)x[2],p,&p2)) { avma=av; return 0; }
! 638: val=1; while (poldivis(p1,p,&p1)&&poldivis(p2,p,&p2)) val++;
! 639: avma = av; return val;
! 640:
! 641: case t_POL:
! 642: if (tp==t_POL)
! 643: {
! 644: v=varn(p); vx=varn(x);
! 645: if (vx == v)
! 646: {
! 647: if ((p>=(GEN)polx && p <= (GEN)(polx+MAXVARN)) || ismonome(p))
! 648: {
! 649: i=2; while (isexactzero((GEN)x[i])) i++;
! 650: return i-2;
! 651: }
! 652: av = avma; limit=stack_lim(av,1);
! 653: for (val=0; ; val++)
! 654: {
! 655: if (!poldivis(x,p,&x)) break;
! 656: if (low_stack(limit, stack_lim(av,1)))
! 657: {
! 658: if(DEBUGMEM>1) err(warnmem,"ggval");
! 659: x = gerepilecopy(av, x);
! 660: }
! 661: }
! 662: avma = av; return val;
! 663: }
! 664: if (vx > v) return 0;
! 665: }
! 666: else if (tp!=t_INT) break;
! 667: i=2; while (isexactzero((GEN)x[i])) i++;
! 668: return minval(x,p,i,lgef(x));
! 669:
! 670: case t_SER:
! 671: if (tp!=t_POL && tp!=t_SER && tp!=t_INT) break;
! 672: v=gvar(p); vx=varn(x);
! 673: if (vx==v) return (long)(valp(x)/ggval(p,polx[v]));
! 674: if (vx>v) return 0;
! 675: return minval(x,p,2,lg(x));
! 676:
! 677: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
! 678: return ggval((GEN)x[1],p) - ggval((GEN)x[2],p);
! 679:
! 680: case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
! 681: return minval(x,p,1,lg(x));
! 682: }
! 683: err(talker,"forbidden or conflicting type in gval");
! 684: return 0; /* not reached */
! 685: }
! 686:
! 687: /* x is non-zero */
! 688: long
! 689: svaluation(ulong x, ulong p, long *py)
! 690: {
! 691: ulong v = 0;
! 692: for(;;)
! 693: {
! 694: if (x%p) { *py = x; return v; }
! 695: v++; x/=p;
! 696: }
! 697: }
! 698:
! 699: /* x is a non-zero integer */
! 700: long
! 701: pvaluation(GEN x, GEN p, GEN *py)
! 702: {
! 703: long av,v;
! 704: GEN p1,p2;
! 705:
! 706: if (egalii(p,gdeux))
! 707: {
! 708: v = vali(x);
! 709: if (py) *py = shifti(x, -v);
! 710: return v;
! 711: }
! 712: if (is_pm1(p))
! 713: {
! 714: v = (signe(p) < 0 && signe(x) < 0);
! 715: if (py) { *py = v? negi(x): icopy(x); }
! 716: return v;
! 717: }
! 718: if (!is_bigint(x))
! 719: {
! 720: long y;
! 721: if (!is_bigint(p))
! 722: {
! 723: v = svaluation(x[2],p[2], &y);
! 724: if (signe(x) < 0) y = -y;
! 725: if (py) *py = stoi(y);
! 726: }
! 727: else
! 728: {
! 729: v = 0;
! 730: if (py) *py = icopy(x);
! 731: }
! 732: return v;
! 733: }
! 734: av = avma; v = 0; (void)new_chunk(lgefint(x));
! 735: for(;;)
! 736: {
! 737: p1 = dvmdii(x,p,&p2);
! 738: if (p2 != gzero) { avma=av; if (py) *py = icopy(x); return v; }
! 739: v++; x = p1;
! 740: }
! 741: }
! 742:
! 743: /*******************************************************************/
! 744: /* */
! 745: /* NEGATION: Create -x */
! 746: /* */
! 747: /*******************************************************************/
! 748:
! 749: GEN
! 750: gneg(GEN x)
! 751: {
! 752: long tx=typ(x),lx,i;
! 753: GEN y;
! 754:
! 755: if (gcmp0(x)) return gcopy(x);
! 756: switch(tx)
! 757: {
! 758: case t_INT: case t_REAL:
! 759: return mpneg(x);
! 760:
! 761: case t_INTMOD: y=cgetg(3,t_INTMOD);
! 762: icopyifstack(x[1],y[1]);
! 763: y[2] = lsubii((GEN)y[1],(GEN)x[2]);
! 764: break;
! 765:
! 766: case t_POLMOD: y=cgetg(3,t_POLMOD);
! 767: copyifstack(x[1],y[1]);
! 768: y[2]=lneg((GEN)x[2]); break;
! 769:
! 770: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
! 771: y=cgetg(3,tx);
! 772: y[1]=lneg((GEN)x[1]);
! 773: y[2]=lcopy((GEN)x[2]); break;
! 774:
! 775: case t_PADIC:
! 776: y=cgetp2(x,valp(x));
! 777: y[4]=lsubii((GEN)x[3],(GEN)x[4]);
! 778: break;
! 779:
! 780: case t_QUAD:
! 781: y=cgetg(4,t_QUAD); copyifstack(x[1],y[1]);
! 782: y[2]=lneg((GEN)x[2]);
! 783: y[3]=lneg((GEN)x[3]); break;
! 784:
! 785: case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
! 786: lx=lg(x); y=cgetg(lx,tx);
! 787: for (i=1; i<lx; i++) y[i]=lneg((GEN)x[i]);
! 788: break;
! 789:
! 790: case t_POL:
! 791: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
! 792: for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
! 793: break;
! 794:
! 795: case t_SER:
! 796: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 797: for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
! 798: break;
! 799:
! 800: default:
! 801: err(typeer,"negation");
! 802: return NULL; /* not reached */
! 803: }
! 804: return y;
! 805: }
! 806:
! 807: GEN
! 808: gneg_i(GEN x)
! 809: {
! 810: long tx=typ(x),lx,i;
! 811: GEN y;
! 812:
! 813: if (gcmp0(x)) return x;
! 814: switch(tx)
! 815: {
! 816: case t_INT: case t_REAL:
! 817: return mpneg(x);
! 818:
! 819: case t_INTMOD: y=cgetg(3,t_INTMOD); y[1]=x[1];
! 820: y[2] = lsubii((GEN)y[1],(GEN)x[2]);
! 821: break;
! 822:
! 823: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
! 824: y=cgetg(3,tx); y[2]=x[2];
! 825: y[1]=(long)gneg_i((GEN)x[1]); break;
! 826:
! 827: case t_PADIC: y = cgetg(5,t_PADIC); y[2]=x[2]; y[3]=x[3];
! 828: y[1] = evalprecp(precp(x)) | evalvalp(valp(x));
! 829: y[4]=lsubii((GEN)x[3],(GEN)x[4]); break;
! 830:
! 831: case t_POLMOD: y=cgetg(3,t_POLMOD); y[1]=x[1];
! 832: y[2]=(long)gneg_i((GEN)x[2]); break;
! 833:
! 834: case t_QUAD: y=cgetg(4,t_QUAD); y[1]=x[1];
! 835: y[2]=(long)gneg_i((GEN)x[2]);
! 836: y[3]=(long)gneg_i((GEN)x[3]); break;
! 837:
! 838: case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
! 839: lx=lg(x); y=cgetg(lx,tx);
! 840: for (i=1; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
! 841: break;
! 842:
! 843: case t_POL: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
! 844: for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
! 845: break;
! 846:
! 847: case t_SER: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 848: for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
! 849: break;
! 850:
! 851: default:
! 852: err(typeer,"negation");
! 853: return NULL; /* not reached */
! 854: }
! 855: return y;
! 856: }
! 857:
! 858: /******************************************************************/
! 859: /* */
! 860: /* ABSOLUTE VALUE */
! 861: /* Create abs(x) if x is integer, real, fraction or complex. */
! 862: /* Error otherwise. */
! 863: /* */
! 864: /******************************************************************/
! 865:
! 866: GEN
! 867: gabs(GEN x, long prec)
! 868: {
! 869: long tx=typ(x),lx,i,l,tetpil;
! 870: GEN y,p1;
! 871:
! 872: switch(tx)
! 873: {
! 874: case t_INT: case t_REAL:
! 875: return mpabs(x);
! 876:
! 877: case t_FRAC: case t_FRACN: y=cgetg(lg(x),tx);
! 878: y[1]=labsi((GEN)x[1]);
! 879: y[2]=labsi((GEN)x[2]); return y;
! 880:
! 881: case t_COMPLEX:
! 882: l=avma; p1=gnorm(x);
! 883: switch(typ(p1))
! 884: {
! 885: case t_INT:
! 886: if (!carrecomplet(p1, &y)) break;
! 887: return gerepileupto(l, y);
! 888: case t_FRAC:
! 889: case t_FRACN:
! 890: {
! 891: GEN a,b;
! 892: if (!carrecomplet((GEN)p1[1], &a)) break;
! 893: if (!carrecomplet((GEN)p1[2], &b)) break;
! 894: return gerepileupto(l, gdiv(a,b));
! 895: }
! 896: }
! 897: tetpil=avma;
! 898: return gerepile(l,tetpil,gsqrt(p1,prec));
! 899:
! 900: case t_QUAD:
! 901: l=avma; p1=gmul(x, realun(prec)); tetpil=avma;
! 902: return gerepile(l,tetpil,gabs(p1,prec));
! 903:
! 904: case t_POL:
! 905: lx=lgef(x); if (lx<=2) return gcopy(x);
! 906: p1=(GEN)x[lx-1];
! 907: switch(typ(p1))
! 908: {
! 909: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 910: if (gsigne(p1)<0) return gneg(x);
! 911: }
! 912: return gcopy(x);
! 913:
! 914: case t_SER:
! 915: if (valp(x) || !signe(x) || lg(x)<3)
! 916: err(talker,"abs is not analytic at 0");
! 917: if (gsigne((GEN) x[2])<0) return gneg(x);
! 918: return gcopy(x);
! 919:
! 920: case t_VEC: case t_COL: case t_MAT:
! 921: lx=lg(x); y=cgetg(lx,tx);
! 922: for (i=1; i<lx; i++)
! 923: y[i]=(long)gabs((GEN)x[i],prec);
! 924: return y;
! 925: }
! 926: err(typeer,"gabs");
! 927: return NULL; /* not reached */
! 928: }
! 929:
! 930: GEN
! 931: gmax(GEN x, GEN y)
! 932: {
! 933: if (gcmp(x,y)>0) y=x; return gcopy(y);
! 934: }
! 935:
! 936: GEN
! 937: gmin(GEN x, GEN y)
! 938: {
! 939: if (gcmp(x,y)<0) y=x; return gcopy(y);
! 940: }
! 941:
! 942: GEN
! 943: vecmax(GEN x)
! 944: {
! 945: long tx=typ(x),lx,lx2,i,j;
! 946: GEN *p1,s;
! 947:
! 948: if (!is_matvec_t(tx)) return gcopy(x);
! 949: lx=lg(x); if (lx==1) return stoi(-VERYBIGINT);
! 950: if (tx!=t_MAT)
! 951: {
! 952: s=(GEN)x[1];
! 953: for (i=2; i<lx; i++)
! 954: if (gcmp((GEN)x[i],s) > 0) s=(GEN)x[i];
! 955: }
! 956: else
! 957: {
! 958: lx2 = lg(x[1]);
! 959: if (lx2==1) return stoi(-VERYBIGINT);
! 960: s=gmael(x,1,1); i=2;
! 961: for (j=1; j<lx; j++)
! 962: {
! 963: p1 = (GEN *) x[j];
! 964: for (; i<lx2; i++)
! 965: if (gcmp(p1[i],s) > 0) s=p1[i];
! 966: i=1;
! 967: }
! 968: }
! 969: return gcopy(s);
! 970: }
! 971:
! 972: GEN
! 973: vecmin(GEN x)
! 974: {
! 975: long tx=typ(x),lx,lx2,i,j;
! 976: GEN *p1,s;
! 977:
! 978: if (!is_matvec_t(tx)) return gcopy(x);
! 979: lx=lg(x); if (lx==1) return stoi(VERYBIGINT);
! 980: if (tx!=t_MAT)
! 981: {
! 982: s=(GEN)x[1];
! 983: for (i=2; i<lx; i++)
! 984: if (gcmp((GEN)x[i],s) < 0) s=(GEN)x[i];
! 985: }
! 986: else
! 987: {
! 988: lx2 = lg(x[1]);
! 989: if (lx2==1) return stoi(VERYBIGINT);
! 990: s=gmael(x,1,1); i=2;
! 991: for (j=1; j<lx; j++)
! 992: {
! 993: p1 = (GEN *) x[j];
! 994: for (; i<lx2; i++)
! 995: if (gcmp(p1[i],s) < 0) s=p1[i];
! 996: i=1;
! 997: }
! 998: }
! 999: return gcopy(s);
! 1000: }
! 1001:
! 1002: /*******************************************************************/
! 1003: /* */
! 1004: /* AFFECT long --> GEN */
! 1005: /* affect long s to GEN x. Useful for initialization. */
! 1006: /* */
! 1007: /*******************************************************************/
! 1008:
! 1009: static void
! 1010: padicaff0(GEN x)
! 1011: {
! 1012: if (signe(x[4]))
! 1013: {
! 1014: setvalp(x, valp(x)|precp(x));
! 1015: affsi(0,(GEN)x[4]);
! 1016: }
! 1017: }
! 1018:
! 1019: void
! 1020: gaffsg(long s, GEN x)
! 1021: {
! 1022: long l,i,v;
! 1023: GEN p;
! 1024:
! 1025: switch(typ(x))
! 1026: {
! 1027: case t_INT:
! 1028: affsi(s,x); break;
! 1029:
! 1030: case t_REAL:
! 1031: affsr(s,x); break;
! 1032:
! 1033: case t_INTMOD:
! 1034: modsiz(s,(GEN)x[1],(GEN)x[2]); break;
! 1035:
! 1036: case t_FRAC: case t_FRACN:
! 1037: affsi(s,(GEN)x[1]); affsi(1,(GEN)x[2]); break;
! 1038:
! 1039: case t_COMPLEX:
! 1040: gaffsg(s,(GEN)x[1]); gaffsg(0,(GEN)x[2]); break;
! 1041:
! 1042: case t_PADIC:
! 1043: if (!s) { padicaff0(x); break; }
! 1044: p = (GEN)x[2];
! 1045: v = (is_bigint(p))? 0: svaluation(s,p[2], &i);
! 1046: setvalp(x,v); modsiz(i,(GEN)x[3],(GEN)x[4]);
! 1047: break;
! 1048:
! 1049: case t_QUAD:
! 1050: gaffsg(s,(GEN)x[2]); gaffsg(0,(GEN)x[3]); break;
! 1051:
! 1052: case t_POLMOD:
! 1053: gaffsg(s,(GEN)x[2]); break;
! 1054:
! 1055: case t_POL:
! 1056: v=varn(x);
! 1057: if (!s) x[1]=evallgef(2) | evalvarn(v);
! 1058: else
! 1059: {
! 1060: x[1]=evalsigne(1) | evallgef(3) | evalvarn(v);
! 1061: gaffsg(s,(GEN)x[2]);
! 1062: }
! 1063: break;
! 1064:
! 1065: case t_SER:
! 1066: v=varn(x); gaffsg(s,(GEN)x[2]); l=lg(x);
! 1067: if (!s)
! 1068: x[1] = evalvalp(l-2) | evalvarn(v);
! 1069: else
! 1070: x[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
! 1071: for (i=3; i<l; i++) gaffsg(0,(GEN)x[i]);
! 1072: break;
! 1073:
! 1074: case t_RFRAC: case t_RFRACN:
! 1075: gaffsg(s,(GEN)x[1]); gaffsg(1,(GEN)x[2]); break;
! 1076:
! 1077: case t_VEC: case t_COL: case t_MAT:
! 1078: if (lg(x)!=2) err(operi,"",t_INT,typ(x));
! 1079: gaffsg(s,(GEN)x[1]); break;
! 1080:
! 1081: default: err(operf,"",t_INT,typ(x));
! 1082: }
! 1083: }
! 1084:
! 1085: /*******************************************************************/
! 1086: /* */
! 1087: /* GENERIC AFFECTATION */
! 1088: /* Affect the content of x to y, whenever possible */
! 1089: /* */
! 1090: /*******************************************************************/
! 1091: void
! 1092: gaffect(GEN x, GEN y)
! 1093: {
! 1094: long i,j,k,l,v,vy,lx,ly,tx,ty,av;
! 1095: GEN p1,num,den;
! 1096:
! 1097:
! 1098: tx=typ(x); ty=typ(y); lx=lg(x); ly=lg(y);
! 1099: if (is_scalar_t(tx))
! 1100: {
! 1101: if (is_scalar_t(ty))
! 1102: {
! 1103: switch(tx)
! 1104: {
! 1105: case t_INT:
! 1106: switch(ty)
! 1107: {
! 1108: case t_INT:
! 1109: /* y = gzero, gnil, gun or gdeux */
! 1110: if (is_universal_constant(y))
! 1111: {
! 1112: if (y==gzero) err(overwriter,"gaffect (gzero)");
! 1113: if (y==gun) err(overwriter,"gaffect (gun)");
! 1114: if (y==gdeux) err(overwriter,"gaffect (gdeux)");
! 1115: err(overwriter,"gaffect (gnil)");
! 1116: }
! 1117: affii(x,y); break;
! 1118:
! 1119: case t_REAL:
! 1120: if (y == gpi) err(overwriter,"gaffect (gpi)");
! 1121: if (y==geuler) err(overwriter,"gaffect (geuler)");
! 1122: affir(x,y); break;
! 1123:
! 1124: case t_INTMOD:
! 1125: modiiz(x,(GEN)y[1],(GEN)y[2]); break;
! 1126:
! 1127: case t_FRAC:
! 1128: if (y == ghalf) err(overwriter,"gaffect (ghalf)");
! 1129: case t_FRACN:
! 1130: affii(x,(GEN)y[1]); affsi(1,(GEN)y[2]); break;
! 1131:
! 1132: case t_COMPLEX:
! 1133: if (y == gi) err(overwriter,"gaffect (gi)");
! 1134: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
! 1135:
! 1136: case t_PADIC:
! 1137: if (!signe(x)) { padicaff0(y); break; }
! 1138: l=avma;
! 1139: setvalp(y, pvaluation(x,(GEN)y[2],&p1));
! 1140: modiiz(p1,(GEN)y[3],(GEN)y[4]);
! 1141: avma=l; break;
! 1142:
! 1143: case t_QUAD:
! 1144: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
! 1145:
! 1146: case t_POLMOD:
! 1147: gaffect(x,(GEN)y[2]); break;
! 1148:
! 1149: default: err(operf,"",tx,ty);
! 1150: }
! 1151: break;
! 1152:
! 1153: case t_REAL:
! 1154: switch(ty)
! 1155: {
! 1156: case t_REAL:
! 1157: affrr(x,y); break;
! 1158:
! 1159: case t_COMPLEX:
! 1160: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
! 1161:
! 1162: case t_POLMOD:
! 1163: gaffect(x,(GEN)y[2]); break;
! 1164:
! 1165: case t_INT: case t_INTMOD: case t_FRAC:
! 1166: case t_FRACN: case t_PADIC: case t_QUAD:
! 1167: err(operf,"",tx,ty);
! 1168:
! 1169: default: err(operf,"",tx,ty);
! 1170: }
! 1171: break;
! 1172:
! 1173: case t_INTMOD:
! 1174: switch(ty)
! 1175: {
! 1176: case t_INTMOD:
! 1177: if (!divise((GEN)x[1],(GEN)y[1]))
! 1178: err(operi,"",tx,ty);
! 1179: modiiz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;
! 1180:
! 1181: case t_POLMOD:
! 1182: gaffect(x,(GEN)y[2]); break;
! 1183:
! 1184: default: err(operf,"",tx,ty);
! 1185: }
! 1186: break;
! 1187:
! 1188: case t_FRAC:
! 1189: switch(ty)
! 1190: {
! 1191: case t_INT:
! 1192: if (! mpdivis((GEN)x[1],(GEN)x[2],y))
! 1193: err(operi,"",tx,ty);
! 1194: break;
! 1195:
! 1196: case t_REAL:
! 1197: diviiz((GEN)x[1],(GEN)x[2],y); break;
! 1198:
! 1199: case t_INTMOD:
! 1200: av=avma; p1=mpinvmod((GEN)x[2],(GEN)y[1]);
! 1201: modiiz(mulii((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
! 1202: avma=av; break;
! 1203:
! 1204: case t_FRAC: case t_FRACN:
! 1205: affii((GEN)x[1],(GEN)y[1]);
! 1206: affii((GEN)x[2],(GEN)y[2]);
! 1207: break;
! 1208:
! 1209: case t_COMPLEX:
! 1210: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
! 1211:
! 1212: case t_PADIC:
! 1213: if (!signe(x[1])) { padicaff0(y); break; }
! 1214: av=avma; num=(GEN)x[1]; den=(GEN)x[2];
! 1215: v = pvaluation(num, (GEN) y[2], &num);
! 1216: if (!v) v = -pvaluation(den,(GEN)y[2],&den);
! 1217: p1=mulii(num,mpinvmod(den,(GEN)y[3]));
! 1218: modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
! 1219: setvalp(y,v); break;
! 1220:
! 1221: case t_QUAD:
! 1222: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
! 1223: case t_POLMOD:
! 1224: gaffect(x,(GEN)y[2]); break;
! 1225: default: err(operf,"",tx,ty);
! 1226: }
! 1227: break;
! 1228:
! 1229: case t_FRACN:
! 1230: switch(ty)
! 1231: {
! 1232: case t_INT:
! 1233: if (! mpdivis((GEN)x[1],(GEN)x[2],y)) err(operi,"",tx,ty);
! 1234: break;
! 1235:
! 1236: case t_REAL:
! 1237: diviiz((GEN)x[1],(GEN)x[2],y); break;
! 1238: case t_INTMOD:
! 1239: av=avma; gaffect(gred(x),y); avma=av; break;
! 1240: case t_FRAC:
! 1241: gredz(x,y); break;
! 1242: case t_FRACN:
! 1243: affii((GEN)x[1],(GEN)y[1]);
! 1244: affii((GEN)x[2],(GEN)y[2]); break;
! 1245: case t_COMPLEX:
! 1246: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
! 1247: case t_PADIC:
! 1248: if (!signe(x[1])) { padicaff0(y); break; }
! 1249: av=avma; num=(GEN)x[1]; den=(GEN)x[2];
! 1250: v = pvaluation(num,(GEN)y[2],&num)
! 1251: - pvaluation(den,(GEN)y[2],&den);
! 1252: p1=mulii(num,mpinvmod(den,(GEN)y[3]));
! 1253: modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
! 1254: setvalp(y,v); break;
! 1255:
! 1256: case t_QUAD:
! 1257: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
! 1258: case t_POLMOD:
! 1259: gaffect(x,(GEN)y[2]); break;
! 1260: default: err(operf,"",tx,ty);
! 1261: }
! 1262: break;
! 1263:
! 1264: case t_COMPLEX:
! 1265: switch(ty)
! 1266: {
! 1267: case t_INT: case t_REAL: case t_INTMOD:
! 1268: case t_FRAC: case t_FRACN: case t_PADIC: case t_QUAD:
! 1269: if (!gcmp0((GEN)x[2])) err(operi,"",tx,ty);
! 1270: gaffect((GEN)x[1],y); break;
! 1271:
! 1272: case t_COMPLEX:
! 1273: gaffect((GEN)x[1],(GEN)y[1]);
! 1274: gaffect((GEN)x[2],(GEN)y[2]);
! 1275: break;
! 1276:
! 1277: case t_POLMOD:
! 1278: gaffect(x,(GEN)y[2]); break;
! 1279:
! 1280: default: err(operf,"",tx,ty);
! 1281: }
! 1282: break;
! 1283:
! 1284: case t_PADIC:
! 1285: switch(ty)
! 1286: {
! 1287: case t_INTMOD:
! 1288: if (valp(x)<0) err(operi,"",tx,ty);
! 1289: av=avma;
! 1290: v = pvaluation((GEN)y[1],(GEN)x[2],&p1);
! 1291: k = signe(x[4])? (precp(x)+valp(x)): valp(x)+1;
! 1292: if (!gcmp1(p1) || v > k) err(operi,"",tx,ty);
! 1293: modiiz(gtrunc(x),(GEN)y[1],(GEN)y[2]); avma=av; break;
! 1294:
! 1295: case t_PADIC:
! 1296: if (!egalii((GEN)x[2],(GEN)y[2])) err(operi,"",tx,ty);
! 1297: modiiz((GEN)x[4],(GEN)y[3],(GEN)y[4]);
! 1298: setvalp(y,valp(x)); break;
! 1299:
! 1300: case t_POLMOD:
! 1301: gaffect(x,(GEN)y[2]); break;
! 1302:
! 1303: default: err(operf,"",tx,ty);
! 1304: }
! 1305: break;
! 1306:
! 1307: case t_QUAD:
! 1308: switch(ty)
! 1309: {
! 1310: case t_INT: case t_INTMOD: case t_FRAC:
! 1311: case t_FRACN: case t_PADIC:
! 1312: if (!gcmp0((GEN)x[3])) err(operi,"",tx,ty);
! 1313: gaffect((GEN)x[2],y); break;
! 1314:
! 1315: case t_REAL:
! 1316: av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av; break;
! 1317: case t_COMPLEX:
! 1318: ly=precision(y);
! 1319: if (ly)
! 1320: {
! 1321: av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av;
! 1322: }
! 1323: else
! 1324: {
! 1325: if (!gcmp0((GEN)x[3])) err(operi,"",tx,ty);
! 1326: gaffect((GEN)x[2],y);
! 1327: }
! 1328: break;
! 1329:
! 1330: case t_QUAD:
! 1331: if (! gegal((GEN)x[1],(GEN)y[1])) err(operi,"",tx,ty);
! 1332: affii((GEN)x[2],(GEN)y[2]);
! 1333: affii((GEN)x[3],(GEN)y[3]);
! 1334: break;
! 1335: case t_POLMOD:
! 1336: gaffect(x,(GEN)y[2]); break;
! 1337: default: err(operf,"",tx,ty);
! 1338: }
! 1339: break;
! 1340:
! 1341: case t_POLMOD:
! 1342: if (ty!=t_POLMOD) err(operf,"",tx,ty);
! 1343: if (! gdivise((GEN)x[1],(GEN)y[1])) err(operi,"",tx,ty);
! 1344: gmodz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;
! 1345:
! 1346: default: err(operf,"",tx,ty);
! 1347: }
! 1348: return;
! 1349: }
! 1350:
! 1351: /* here y is not scalar */
! 1352: switch(ty)
! 1353: {
! 1354: case t_POL:
! 1355: v=varn(y);
! 1356: if (y==polun[v] || y==polx[v])
! 1357: err(talker,"trying to overwrite a universal polynomial");
! 1358: gaffect(x,(GEN)y[2]);
! 1359: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1360: if (gcmp0(x)) y[1]=evallgef(2)+evalvarn(v);
! 1361: else y[1]=evalsigne(1)+evallgef(3)+evalvarn(v);
! 1362: break;
! 1363:
! 1364: case t_SER:
! 1365: v=varn(y); gaffect(x,(GEN)y[2]);
! 1366: if (gcmp0(x))
! 1367: y[1] = evalvalp(ly-2) | evalvarn(v);
! 1368: else
! 1369: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
! 1370: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1371: break;
! 1372:
! 1373: case t_RFRAC: case t_RFRACN:
! 1374: gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;
! 1375:
! 1376: default: err(operf,"",tx,ty);
! 1377: }
! 1378: return;
! 1379: }
! 1380:
! 1381: if (is_const_t(ty)) {
! 1382: #define initial_value(ep) ((ep)+1)
! 1383: if (tx == t_POL) {
! 1384: entree *var = varentries[ordvar[varn(x)]];
! 1385: /* Is a bound expression, thus should not loop: */
! 1386: if (var && var->value != (void*)initial_value(var)) {
! 1387: gaffect(geval(x),y);
! 1388: return;
! 1389: }
! 1390: } else if (is_rfrac_t(tx)) {
! 1391: GEN num = (GEN)x[1], den = (GEN)x[2];
! 1392: entree *varnum = varentries[ordvar[varn(num)]];
! 1393: entree *varden = varentries[ordvar[varn(den)]];
! 1394:
! 1395: /* Are bound expressions, thus should not loop: */
! 1396: if (varnum && varnum->value != (void*)initial_value(varnum)
! 1397: && varden && varden->value != (void*)initial_value(varden)) {
! 1398: gaffect(geval(x),y);
! 1399: return;
! 1400: }
! 1401: }
! 1402: #undef initial_value
! 1403: err(operf,"",tx,ty);
! 1404: }
! 1405:
! 1406: switch(tx)
! 1407: {
! 1408: case t_POL:
! 1409: v=varn(x);
! 1410: switch(ty)
! 1411: {
! 1412: case t_POL:
! 1413: vy=varn(y); if (vy>v) err(operf,"",tx,ty);
! 1414: if (vy==v)
! 1415: {
! 1416: l=lgef(x); if (l>ly) err(operi,"",tx,ty);
! 1417: y[1]=x[1]; for (i=2; i<l; i++) gaffect((GEN)x[i],(GEN)y[i]);
! 1418: }
! 1419: else
! 1420: {
! 1421: gaffect(x,(GEN)y[2]);
! 1422: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1423: if (signe(x)) y[1]=evalsigne(1)+evallgef(3)+evalvarn(vy);
! 1424: else y[1]=evallgef(2)+evalvarn(vy);
! 1425: }
! 1426: break;
! 1427:
! 1428: case t_SER:
! 1429: vy=varn(y); if (vy>v) err(operf,"",tx,ty);
! 1430: if (!signe(x)) { gaffsg(0,y); return; }
! 1431: if (vy==v)
! 1432: {
! 1433: i=gval(x,v); y[1]=evalvarn(v) | evalvalp(i) | evalsigne(1);
! 1434: k=lgef(x)-i; if (k>ly) k=ly;
! 1435: for (j=2; j<k; j++) gaffect((GEN)x[i+j],(GEN)y[j]);
! 1436: for ( ; j<ly; j++) gaffsg(0,(GEN)y[j]);
! 1437: }
! 1438: else
! 1439: {
! 1440: gaffect(x,(GEN)y[2]);
! 1441: if (!signe(x))
! 1442: y[1] = evalvalp(ly-2) | evalvarn(vy);
! 1443: else
! 1444: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
! 1445: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1446: }
! 1447: break;
! 1448:
! 1449: case t_POLMOD:
! 1450: gmodz(x,(GEN)y[1],(GEN)y[2]); break;
! 1451:
! 1452: case t_RFRAC: case t_RFRACN:
! 1453: gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;
! 1454:
! 1455: default: err(operf,"",tx,ty);
! 1456: }
! 1457: break;
! 1458:
! 1459: case t_SER:
! 1460: if (ty!=t_SER) err(operf,"",tx,ty);
! 1461: v=varn(x); vy=varn(y); if (vy>v) err(operf,"",tx,ty);
! 1462: if (vy==v)
! 1463: {
! 1464: y[1]=x[1]; k=lx; if (k>ly) k=ly;
! 1465: for (i=2; i< k; i++) gaffect((GEN)x[i],(GEN)y[i]);
! 1466: for ( ; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1467: }
! 1468: else
! 1469: {
! 1470: gaffect(x,(GEN)y[2]);
! 1471: if (!signe(x))
! 1472: y[1] = evalvalp(ly-2) | evalvarn(vy);
! 1473: else
! 1474: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
! 1475: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
! 1476: }
! 1477: break;
! 1478:
! 1479: case t_RFRAC: case t_RFRACN:
! 1480: switch(ty)
! 1481: {
! 1482: case t_POL: case t_VEC: case t_COL: case t_MAT:
! 1483: err(operf,"",tx,ty);
! 1484:
! 1485: case t_POLMOD:
! 1486: av=avma; p1=ginvmod((GEN)x[2],(GEN)y[1]);
! 1487: gmodz(gmul((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
! 1488: avma=av; break;
! 1489:
! 1490: case t_SER:
! 1491: gdivz((GEN)x[1],(GEN)x[2],y); break;
! 1492:
! 1493: case t_RFRAC:
! 1494: if (tx==t_RFRACN) gredz(x,y);
! 1495: else
! 1496: {
! 1497: gaffect((GEN)x[1],(GEN)y[1]);
! 1498: gaffect((GEN)x[2],(GEN)y[2]);
! 1499: }
! 1500: break;
! 1501:
! 1502: case t_RFRACN:
! 1503: gaffect((GEN)x[1],(GEN)y[1]);
! 1504: gaffect((GEN)x[2],(GEN)y[2]); break;
! 1505:
! 1506: default: err(operf,"",tx,ty);
! 1507: }
! 1508: break;
! 1509:
! 1510: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
! 1511: if (ty != tx || ly != lx) err(operi,"",tx,ty);
! 1512: for (i=1; i<lx; i++) gaffect((GEN)x[i],(GEN)y[i]);
! 1513: break;
! 1514:
! 1515: default: err(operf,"",tx,ty);
! 1516: }
! 1517: }
! 1518:
! 1519: /*******************************************************************/
! 1520: /* */
! 1521: /* CONVERSION QUAD --> REAL, COMPLEX OR P-ADIC */
! 1522: /* */
! 1523: /*******************************************************************/
! 1524:
! 1525: GEN
! 1526: co8(GEN x, long prec)
! 1527: {
! 1528: long av=avma,tetpil;
! 1529: GEN p1, p = (GEN) x[1];
! 1530:
! 1531: p1 = subii(sqri((GEN)p[3]), shifti((GEN)p[2],2));
! 1532: if (signe(p1) > 0)
! 1533: {
! 1534: p1 = subri(gsqrt(p1,prec), (GEN)p[3]);
! 1535: setexpo(p1, expo(p1)-1);
! 1536: }
! 1537: else
! 1538: {
! 1539: p1 = gsub(gsqrt(p1,prec), (GEN)p[3]);
! 1540: p1[1] = lmul2n((GEN)p1[1],-1);
! 1541: setexpo(p1[2], expo(p1[2])-1);
! 1542: }/* p1 = (-b + sqrt(D)) / 2 */
! 1543: p1 = gmul((GEN)x[3],p1); tetpil=avma;
! 1544: return gerepile(av,tetpil,gadd((GEN)x[2],p1));
! 1545: }
! 1546:
! 1547: GEN
! 1548: cvtop(GEN x, GEN p, long l)
! 1549: {
! 1550: GEN p1,p2,p3;
! 1551: long av,tetpil,n;
! 1552:
! 1553: if (typ(p)!=t_INT)
! 1554: err(talker,"not an integer modulus in cvtop or gcvtop");
! 1555: if (gcmp0(x)) return ggrandocp(p,l);
! 1556: switch(typ(x))
! 1557: {
! 1558: case t_INT:
! 1559: return gadd(x,ggrandocp(p,ggval(x,p)+l));
! 1560:
! 1561: case t_INTMOD:
! 1562: n=ggval((GEN)x[1],p); if (n>l) n=l;
! 1563: return gadd((GEN)x[2],ggrandocp(p,n));
! 1564:
! 1565: case t_FRAC: case t_FRACN:
! 1566: n = ggval((GEN)x[1],p) - ggval((GEN)x[2],p);
! 1567: return gadd(x,ggrandocp(p,n+l));
! 1568:
! 1569: case t_COMPLEX:
! 1570: av=avma; p1=gsqrt(gaddgs(ggrandocp(p,l),-1),0);
! 1571: p1=gmul(p1,(GEN)x[2]); tetpil=avma;
! 1572: return gerepile(av,tetpil,gadd(p1,(GEN)x[1]));
! 1573:
! 1574: case t_PADIC:
! 1575: return gprec(x,l);
! 1576:
! 1577: case t_QUAD:
! 1578: av=avma; p1=(GEN)x[1]; p3=gmul2n((GEN)p1[3],-1);
! 1579: p2=gsub(gsqr(p3),(GEN)p1[2]);
! 1580: switch(typ(p2))
! 1581: {
! 1582: case t_INT:
! 1583: n=ggval(p2,p); p2=gadd(p2,ggrandocp(p,n+l)); break;
! 1584:
! 1585: case t_INTMOD: case t_PADIC:
! 1586: break;
! 1587:
! 1588: case t_FRAC: case t_FRACN:
! 1589: n = ggval((GEN)p2[1],p) - ggval((GEN)p2[2],p);
! 1590: p2=gadd(p2,ggrandocp(p,n+l)); break;
! 1591:
! 1592: default: err(operi,"",t_QUAD,t_QUAD);
! 1593: }
! 1594: p2=gsqrt(p2,0); p1=gmul((GEN)x[3],gsub(p2,p3)); tetpil=avma;
! 1595: return gerepile(av,tetpil,gadd((GEN)x[2],p1));
! 1596:
! 1597: }
! 1598: err(typeer,"cvtop");
! 1599: return NULL; /* not reached */
! 1600: }
! 1601:
! 1602: GEN
! 1603: gcvtop(GEN x, GEN p, long r)
! 1604: {
! 1605: long i,lx, tx=typ(x);
! 1606: GEN y;
! 1607:
! 1608: if (is_const_t(tx)) return cvtop(x,p,r);
! 1609: switch(tx)
! 1610: {
! 1611: case t_POL:
! 1612: lx=lgef(x); y=cgetg(lx,t_POL); y[1]=x[1];
! 1613: for (i=2; i<lx; i++)
! 1614: y[i]=(long)gcvtop((GEN)x[i],p,r);
! 1615: break;
! 1616:
! 1617: case t_SER:
! 1618: lx=lg(x); y=cgetg(lx,t_SER); y[1]=x[1];
! 1619: for (i=2; i<lx; i++)
! 1620: y[i]=(long)gcvtop((GEN)x[i],p,r);
! 1621: break;
! 1622:
! 1623: case t_POLMOD: case t_RFRAC: case t_RFRACN:
! 1624: case t_VEC: case t_COL: case t_MAT:
! 1625: lx=lg(x); y=cgetg(lx,tx);
! 1626: for (i=1; i<lx; i++)
! 1627: y[i]=(long)gcvtop((GEN)x[i],p,r);
! 1628: break;
! 1629:
! 1630: default: err(typeer,"gcvtop");
! 1631: return NULL; /* not reached */
! 1632: }
! 1633: return y;
! 1634: }
! 1635:
! 1636: long
! 1637: gexpo(GEN x)
! 1638: {
! 1639: long tx=typ(x),lx,e,i,y,av;
! 1640:
! 1641: switch(tx)
! 1642: {
! 1643: case t_INT:
! 1644: return expi(x);
! 1645:
! 1646: case t_FRAC: case t_FRACN:
! 1647: if (!signe(x[1])) return -HIGHEXPOBIT;
! 1648: return expi((GEN)x[1]) - expi((GEN)x[2]);
! 1649:
! 1650: case t_REAL:
! 1651: return expo(x);
! 1652:
! 1653: case t_COMPLEX:
! 1654: e = gexpo((GEN)x[1]);
! 1655: y = gexpo((GEN)x[2]); return max(e,y);
! 1656:
! 1657: case t_QUAD:
! 1658: av=avma; e = gexpo(co8(x,3)); avma=av; return e;
! 1659:
! 1660: case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
! 1661: lx=(tx==t_POL)? lgef(x): lg(x);
! 1662: y = -HIGHEXPOBIT;
! 1663: for (i=lontyp[tx]; i<lx; i++) { e=gexpo((GEN)x[i]); if (e>y) y=e; }
! 1664: return y;
! 1665: }
! 1666: err(typeer,"gexpo");
! 1667: return 0; /* not reached */
! 1668: }
! 1669:
! 1670: long
! 1671: sizedigit(GEN x)
! 1672: {
! 1673: return gcmp0(x)? 0: (long) ((gexpo(x)+1) * L2SL10) + 1;
! 1674: }
! 1675:
! 1676: /* Normalize series x in place.
! 1677: * Assumption: x,x[2],...,x[lg(x)-1] have been created in that order.
! 1678: * All intermediate objects will be destroyed.
! 1679: */
! 1680: GEN
! 1681: normalize(GEN x)
! 1682: {
! 1683: long i,j, lx = lg(x);
! 1684:
! 1685: if (typ(x)!=t_SER) err(typeer,"normalize");
! 1686: if (lx==2) { setsigne(x,0); avma = (long) x; return x; }
! 1687: if (! isexactzero((GEN)x[2])) { setsigne(x,1); return x; }
! 1688:
! 1689: for (i=3; i<lx; i++)
! 1690: if (! isexactzero((GEN)x[i]))
! 1691: {
! 1692: long tetpil = avma;
! 1693: GEN p1 = cgetg(lx-i+2,t_SER);
! 1694: p1[1] = evalsigne(1) | evalvalp(valp(x)+i-2) | evalvarn(varn(x));
! 1695: j=i; i=2; while (j<lx) p1[i++] = lcopy((GEN)x[j++]);
! 1696: return gerepile((long) (x+lx),tetpil,p1);
! 1697: }
! 1698: avma = (long) (x+lx); return zeroser(varn(x),lx-2);
! 1699: }
! 1700:
! 1701: GEN
! 1702: normalizepol_i(GEN x, long lx)
! 1703: {
! 1704: long i;
! 1705: for (i=lx-1; i>1; i--)
! 1706: if (! isexactzero((GEN)x[i])) break;
! 1707: setlgef(x,i+1);
! 1708: for (; i>1; i--)
! 1709: if (! gcmp0((GEN)x[i]) ) { setsigne(x,1); return x; }
! 1710: setsigne(x,0); return x;
! 1711: }
! 1712:
! 1713: /* Normalize polynomial x in place. See preceding comment */
! 1714: GEN
! 1715: normalizepol(GEN x)
! 1716: {
! 1717: if (typ(x)!=t_POL) err(typeer,"normalizepol");
! 1718: return normalizepol_i(x, lgef(x));
! 1719: }
! 1720:
! 1721: int
! 1722: gsigne(GEN x)
! 1723: {
! 1724: switch(typ(x))
! 1725: {
! 1726: case t_INT: case t_REAL:
! 1727: return signe(x);
! 1728:
! 1729: case t_FRAC: case t_FRACN:
! 1730: return (signe(x[2])>0) ? signe(x[1]) : -signe(x[1]);
! 1731: }
! 1732: err(typeer,"gsigne");
! 1733: return 0; /* not reached */
! 1734: }
! 1735:
! 1736: /*******************************************************************/
! 1737: /* */
! 1738: /* LISTS */
! 1739: /* */
! 1740: /*******************************************************************/
! 1741:
! 1742: GEN
! 1743: listcreate(long n)
! 1744: {
! 1745: GEN list;
! 1746:
! 1747: if (n<0) err(talker,"negative length in listcreate");
! 1748: n += 2;
! 1749: if (n & ~LGEFBITS) err(talker,"list too long (max = %ld)",LGEFBITS-2);
! 1750: list=cgetg(n,t_LIST); list[1]=evallgef(2);
! 1751: return list;
! 1752: }
! 1753:
! 1754: static void
! 1755: listaffect(GEN list, long index, GEN object)
! 1756: {
! 1757: if (index<lgef(list) && isclone(list[index]))
! 1758: gunclone((GEN)list[index]);
! 1759: list[index]=lclone(object);
! 1760: }
! 1761:
! 1762: void
! 1763: listkill(GEN list)
! 1764: {
! 1765: long lx,i;
! 1766:
! 1767: if (typ(list)!=t_LIST) err(typeer,"listkill");
! 1768: lx=lgef(list);
! 1769: for (i=2; i<lx; i++)
! 1770: if (isclone(list[i])) gunclone((GEN)list[i]);
! 1771: list[1]=evallgef(2); return;
! 1772: }
! 1773:
! 1774: GEN
! 1775: listput(GEN list, GEN object, long index)
! 1776: {
! 1777: long lx=lgef(list);
! 1778:
! 1779: if (typ(list)!=t_LIST) err(typeer,"listput");
! 1780: if (index < 0) err(talker,"negative index (%ld) in listput", index);
! 1781: if (!index || index >= lx-1)
! 1782: {
! 1783: index = lx-1; lx++;
! 1784: if (lx > lg(list))
! 1785: err(talker,"no more room in this list (size %ld)", lg(list)-2);
! 1786: }
! 1787: listaffect(list,index+1,object);
! 1788: list[1]=evallgef(lx);
! 1789: return (GEN)list[index+1];
! 1790: }
! 1791:
! 1792: GEN
! 1793: listinsert(GEN list, GEN object, long index)
! 1794: {
! 1795: long lx = lgef(list), i;
! 1796:
! 1797: if (typ(list)!=t_LIST) err(typeer,"listinsert");
! 1798: if (index <= 0 || index > lx-1) err(talker,"bad index in listinsert");
! 1799: lx++; if (lx > lg(list)) err(talker,"no more room in this list");
! 1800:
! 1801: for (i=lx-2; i > index; i--) list[i+1]=list[i];
! 1802: list[index+1] = lclone(object);
! 1803: list[1] = evallgef(lx); return (GEN)list[index+1];
! 1804: }
! 1805:
! 1806: GEN
! 1807: gtolist(GEN x)
! 1808: {
! 1809: long tx,lx,i;
! 1810: GEN list;
! 1811:
! 1812: if (!x) { list = cgetg(2, t_LIST); list[1] = evallgef(2); return list; }
! 1813:
! 1814: tx = typ(x);
! 1815: lx = (tx==t_LIST)? lgef(x): lg(x);
! 1816: switch(tx)
! 1817: {
! 1818: case t_VEC: case t_COL:
! 1819: lx++; x--; /* fall through */
! 1820: case t_LIST:
! 1821: list = cgetg(lx,t_LIST);
! 1822: for (i=2; i<lx; i++) list[i] = lclone((GEN)x[i]);
! 1823: break;
! 1824: default: err(typeer,"gtolist");
! 1825: return NULL; /* not reached */
! 1826: }
! 1827: list[1] = evallgef(lx); return list;
! 1828: }
! 1829:
! 1830: GEN
! 1831: listconcat(GEN list1, GEN list2)
! 1832: {
! 1833: long i,l1,lx;
! 1834: GEN list;
! 1835:
! 1836: if (typ(list1)!=t_LIST || typ(list2)!=t_LIST)
! 1837: err(typeer,"listconcat");
! 1838: l1=lgef(list1)-2; lx=l1+lgef(list2);
! 1839: list = listcreate(lx-2);
! 1840: for (i=2; i<=l1+1; i++) listaffect(list,i,(GEN)list1[i]);
! 1841: for ( ; i<lx; i++) listaffect(list,i,(GEN)list2[i-l1]);
! 1842: list[1]=evallgef(lx); return list;
! 1843: }
! 1844:
! 1845: GEN
! 1846: listsort(GEN list, long flag)
! 1847: {
! 1848: long i,av=avma, c=list[1], lx = lgef(list)-1;
! 1849: GEN perm,vec,l;
! 1850:
! 1851: if (typ(list) != t_LIST) err(typeer,"listsort");
! 1852: vec = list+1;
! 1853: vec[0] = evaltyp(t_VEC) | evallg(lx);
! 1854: perm = sindexsort(vec);
! 1855: list[1] = c; l = cgetg(lx,t_VEC);
! 1856: for (i=1; i<lx; i++) l[i] = vec[perm[i]];
! 1857: if (flag)
! 1858: {
! 1859: c=1; vec[1] = l[1];
! 1860: for (i=2; i<lx; i++)
! 1861: if (!gegal((GEN)l[i], (GEN)vec[c]))
! 1862: vec[++c] = l[i];
! 1863: else
! 1864: if (isclone(l[i])) gunclone((GEN)l[i]);
! 1865: setlgef(list,c+2);
! 1866: }
! 1867: else
! 1868: for (i=1; i<lx; i++) vec[i] = l[i];
! 1869:
! 1870: avma=av; return list;
! 1871: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>