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