Annotation of OpenXM_contrib/pari/src/basemath/gen3.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** GENERIC OPERATIONS **/
! 4: /** (third part) **/
! 5: /** **/
! 6: /********************************************************************/
! 7: /* $Id: gen3.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: /********************************************************************/
! 11: /** **/
! 12: /** PRINCIPAL VARIABLE NUMBER **/
! 13: /** **/
! 14: /********************************************************************/
! 15:
! 16: int
! 17: gvar(GEN x)
! 18: {
! 19: long tx=typ(x),i,v,w;
! 20:
! 21: if (is_polser_t(tx)) return varn(x);
! 22: if (tx==t_POLMOD) return varn(x[1]);
! 23: if (is_const_t(tx) || is_qf_t(tx) || tx == t_STR || tx == t_LIST)
! 24: return BIGINT;
! 25: v=BIGINT;
! 26: for (i=1; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
! 27: return v;
! 28: }
! 29:
! 30: int
! 31: gvar2(GEN x)
! 32: {
! 33: long tx=typ(x),i,v,w;
! 34:
! 35: if (is_const_t(tx) || is_qf_t(tx)) return BIGINT;
! 36: v = BIGINT;
! 37: switch(tx)
! 38: {
! 39: case t_POLMOD:
! 40: v=gvar2((GEN)x[1]);
! 41: w=gvar2((GEN)x[2]); if (w<v) v=w;
! 42: return v;
! 43:
! 44: case t_SER:
! 45: for (i=2; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
! 46: return v;
! 47:
! 48: case t_POL:
! 49: for (i=2; i<lgef(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
! 50: return v;
! 51: }
! 52: for (i=1; i<lg(x); i++) { w=gvar2((GEN)x[i]); if (w<v) v=w; }
! 53: return v;
! 54: }
! 55:
! 56: GEN
! 57: gpolvar(GEN x)
! 58: {
! 59: if (typ(x)==t_PADIC) x = (GEN)x[2];
! 60: else
! 61: {
! 62: long v=gvar(x);
! 63: if (v==BIGINT) err(typeer,"polvar");
! 64: x = polx[v];
! 65: }
! 66: return gcopy(x);
! 67: }
! 68:
! 69: /*******************************************************************/
! 70: /* */
! 71: /* PRECISION OF SCALAR OBJECTS */
! 72: /* */
! 73: /*******************************************************************/
! 74:
! 75: long
! 76: precision(GEN x)
! 77: {
! 78: long tx=typ(x),k,l;
! 79:
! 80: if (tx==t_REAL)
! 81: {
! 82: k=2-(expo(x)>>TWOPOTBITS_IN_LONG);
! 83: l=lg(x); if (l>k) k=l;
! 84: return k;
! 85: }
! 86: if (tx==t_COMPLEX)
! 87: {
! 88: k=precision((GEN)x[1]);
! 89: l=precision((GEN)x[2]); if (l && l<k) k=l;
! 90: return k;
! 91: }
! 92: return 0;
! 93: }
! 94:
! 95: long
! 96: gprecision(GEN x)
! 97: {
! 98: long tx=typ(x),lx=lg(x),i,k,l;
! 99:
! 100: if (is_scalar_t(tx)) return precision(x);
! 101: switch(tx)
! 102: {
! 103: case t_POL: lx=lgef(x);
! 104: case t_VEC: case t_COL: case t_MAT:
! 105: k=VERYBIGINT;
! 106: for (i=lontyp[tx]; i<lx; i++)
! 107: {
! 108: l = gprecision((GEN)x[i]);
! 109: if (l && l<k) k = l;
! 110: }
! 111: return (k==VERYBIGINT)? 0: k;
! 112:
! 113: case t_RFRAC: case t_RFRACN:
! 114: {
! 115: k=gprecision((GEN)x[1]);
! 116: l=gprecision((GEN)x[2]); if (l && l<k) k=l;
! 117: return k;
! 118: }
! 119: case t_QFR:
! 120: return gprecision((GEN)x[4]);
! 121: }
! 122: return 0;
! 123: }
! 124:
! 125: GEN
! 126: ggprecision(GEN x)
! 127: {
! 128: long a=gprecision(x);
! 129: return stoi(a ? (long) ((a-2)*pariK): VERYBIGINT);
! 130: }
! 131:
! 132: GEN
! 133: precision0(GEN x, long n)
! 134: {
! 135: if (n) return gprec(x,n);
! 136: return ggprecision(x);
! 137: }
! 138:
! 139: /* attention: precision p-adique absolue */
! 140: long
! 141: padicprec(GEN x, GEN p)
! 142: {
! 143: long i,s,t,lx=lg(x),tx=typ(x);
! 144:
! 145: switch(tx)
! 146: {
! 147: case t_INT: case t_FRAC: case t_FRACN:
! 148: return VERYBIGINT;
! 149:
! 150: case t_INTMOD:
! 151: return ggval((GEN)x[1],p);
! 152:
! 153: case t_PADIC:
! 154: if (!gegal((GEN)x[2],p))
! 155: err(talker,"not the same prime in padicprec");
! 156: return precp(x)+valp(x);
! 157:
! 158: case t_POL:
! 159: lx=lgef(x);
! 160:
! 161: case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_SER: case t_RFRAC:
! 162: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 163: for (s=VERYBIGINT, i=lontyp[tx]; i<lx; i++)
! 164: {
! 165: t = padicprec((GEN)x[i],p); if (t<s) s = t;
! 166: }
! 167: return s;
! 168: }
! 169: err(typeer,"padicprec");
! 170: return 0; /* not reached */
! 171: }
! 172:
! 173: /* Degre de x par rapport a la variable v si v>=0, par rapport a la variable
! 174: * principale si v<0. On suppose x fraction rationnelle ou polynome.
! 175: * Convention deg(0)=-1.
! 176: */
! 177: long
! 178: poldegree(GEN x, long v)
! 179: {
! 180: long tx=typ(x), av, w, d;
! 181:
! 182: if (is_scalar_t(tx)) return gcmp0(x)? -1: 0;
! 183: switch(tx)
! 184: {
! 185: case t_POL:
! 186: w = varn(x);
! 187: if (v < 0 || v == w) return lgef(x)-3;
! 188: if (v < w) return signe(x)? 0: -1;
! 189: av = avma; x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
! 190: if (gvar(x)) { d = gcmp0(x)? -1: 0; } else d = lgef(x)-3;
! 191: avma = av; return d;
! 192:
! 193: case t_RFRAC: case t_RFRACN:
! 194: if (gcmp0((GEN)x[1])) return -1;
! 195: return poldegree((GEN)x[1],v) - poldegree((GEN)x[2],v);
! 196: }
! 197: err(typeer,"degree");
! 198: return 0; /* not reached */
! 199: }
! 200:
! 201: long
! 202: degree(GEN x)
! 203: {
! 204: return poldegree(x,-1);
! 205: }
! 206:
! 207: /* si v<0, par rapport a la variable principale, sinon par rapport a v.
! 208: * On suppose que x est un polynome ou une serie.
! 209: */
! 210: GEN
! 211: pollead(GEN x, long v)
! 212: {
! 213: long l,tx = typ(x),av,tetpil,w;
! 214: GEN xinit;
! 215:
! 216: if (is_scalar_t(tx)) return gcopy(x);
! 217: w = varn(x);
! 218: switch(tx)
! 219: {
! 220: case t_POL:
! 221: if (v < 0 || v == w)
! 222: {
! 223: l=lgef(x);
! 224: return (l==2)? gzero: gcopy((GEN)x[l-1]);
! 225: }
! 226: if (v < w) return gcopy(x);
! 227: av = avma; xinit = x;
! 228: x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
! 229: if (gvar(x)) { avma = av; return gcopy(xinit); }
! 230: l = lgef(x); if (l == 2) { avma = av; return gzero; }
! 231: tetpil = avma; x = gsubst((GEN)x[l-1],MAXVARN,polx[w]);
! 232: return gerepile(av,tetpil,x);
! 233:
! 234: case t_SER:
! 235: if (v < 0 || v == w) return (!signe(x))? gzero: gcopy((GEN)x[2]);
! 236: if (v < w) return gcopy(x);
! 237: av = avma; xinit = x;
! 238: x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
! 239: if (gvar(x)) { avma = av; return gcopy(xinit);}
! 240: if (!signe(x)) { avma = av; return gzero;}
! 241: tetpil = avma; x = gsubst((GEN)x[2],MAXVARN,polx[w]);
! 242: return gerepile(av,tetpil,x);
! 243: }
! 244: err(typeer,"pollead");
! 245: return NULL; /* not reached */
! 246: }
! 247:
! 248: /* returns 1 if there's a real component in the structure, 0 otherwise */
! 249: int
! 250: isinexactreal(GEN x)
! 251: {
! 252: long tx=typ(x),lx,i;
! 253:
! 254: if (is_scalar_t(tx))
! 255: {
! 256: switch(tx)
! 257: {
! 258: case t_REAL:
! 259: return 1;
! 260:
! 261: case t_COMPLEX: case t_QUAD:
! 262: return (typ(x[1])==t_REAL || typ(x[2])==t_REAL);
! 263: }
! 264: return 0;
! 265: }
! 266: switch(tx)
! 267: {
! 268: case t_QFR: case t_QFI:
! 269: return 0;
! 270:
! 271: case t_RFRAC: case t_RFRACN:
! 272: return isinexactreal((GEN)x[1]) || isinexactreal((GEN)x[2]);
! 273: }
! 274: lx = (tx==t_POL)? lgef(x): lg(x);
! 275: for (i=lontyp[tx]; i<lx; i++)
! 276: if (isinexactreal((GEN)x[i])) return 1;
! 277: return 0;
! 278: }
! 279:
! 280: int
! 281: isexactzero(GEN g)
! 282: {
! 283: long i;
! 284: switch (typ(g))
! 285: {
! 286: case t_INT:
! 287: return !signe(g);
! 288: case t_REAL: case t_PADIC: case t_SER:
! 289: return 0;
! 290: case t_INTMOD: case t_POLMOD:
! 291: return isexactzero((GEN)g[2]);
! 292: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
! 293: return isexactzero((GEN)g[1]);
! 294: case t_COMPLEX:
! 295: return isexactzero((GEN)g[1]) && isexactzero((GEN)g[2]);
! 296: case t_QUAD:
! 297: return isexactzero((GEN)g[2]) && isexactzero((GEN)g[3]);
! 298:
! 299: case t_POL:
! 300: for (i=lgef(g)-1; i>1; i--)
! 301: if (!isexactzero((GEN)g[i])) return 0;
! 302: return 1;
! 303: case t_VEC: case t_COL: case t_MAT:
! 304: for (i=lg(g)-1; i; i--)
! 305: if (!isexactzero((GEN)g[i])) return 0;
! 306: return 1;
! 307: }
! 308: return 0;
! 309: }
! 310:
! 311: int
! 312: iscomplex(GEN x)
! 313: {
! 314: switch(typ(x))
! 315: {
! 316: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 317: return 0;
! 318: case t_COMPLEX:
! 319: return !gcmp0((GEN)x[2]);
! 320: case t_QUAD:
! 321: return signe(mael(x,1,2)) > 0;
! 322: }
! 323: err(typeer,"iscomplex");
! 324: return 0; /* not reached */
! 325: }
! 326:
! 327: int
! 328: ismonome(GEN x)
! 329: {
! 330: long i;
! 331: if (typ(x)!=t_POL || !signe(x)) return 0;
! 332: for (i=lgef(x)-2; i>1; i--)
! 333: if (!isexactzero((GEN)x[i])) return 0;
! 334: return 1;
! 335: }
! 336:
! 337: /********************************************************************/
! 338: /** **/
! 339: /** MULTIPLICATION SIMPLE **/
! 340: /** **/
! 341: /********************************************************************/
! 342: #define fix_frac(z) if (signe(z[2])<0)\
! 343: {\
! 344: setsigne(z[1],-signe(z[1]));\
! 345: setsigne(z[2],1);\
! 346: }
! 347:
! 348: /* assume z[1] was created last */
! 349: #define fix_frac_if_int(z) if (is_pm1(z[2]))\
! 350: z = gerepileupto((long)(z+3), (GEN)z[1]);
! 351:
! 352: GEN
! 353: gmulsg(long s, GEN y)
! 354: {
! 355: long ty=typ(y),ly=lg(y),i,av,tetpil;
! 356: GEN z,p1,p2;
! 357:
! 358: switch(ty)
! 359: {
! 360: case t_INT:
! 361: return mulsi(s,y);
! 362:
! 363: case t_REAL:
! 364: return mulsr(s,y);
! 365:
! 366: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
! 367: (void)new_chunk(lgefint(p2)<<2);
! 368: p1=mulsi(s,(GEN)y[2]); avma=(long)z;
! 369: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 370:
! 371: case t_FRAC:
! 372: if (!s) return gzero;
! 373: z = cgetg(3,t_FRAC);
! 374: i = cgcd(s, smodis((GEN)y[2], s));
! 375: if (i == 1)
! 376: {
! 377: z[2] = licopy((GEN)y[2]);
! 378: z[1] = lmulis((GEN)y[1], s);
! 379: }
! 380: else
! 381: {
! 382: z[2] = ldivis((GEN)y[2], i);
! 383: z[1] = lmulis((GEN)y[1], s/i);
! 384: fix_frac_if_int(z);
! 385: }
! 386: return z;
! 387:
! 388: case t_FRACN: z=cgetg(3,ty);
! 389: z[1]=lmulsi(s,(GEN)y[1]);
! 390: z[2]=licopy((GEN)y[2]);
! 391: return z;
! 392:
! 393: case t_COMPLEX: z=cgetg(ly,ty);
! 394: z[1]=lmulsg(s,(GEN)y[1]);
! 395: z[2]=lmulsg(s,(GEN)y[2]); return z;
! 396:
! 397: case t_PADIC:
! 398: if (!s) return gzero;
! 399: av=avma; p1=cgetp(y); gaffsg(s,p1); tetpil=avma;
! 400: return gerepile(av,tetpil,gmul(p1,y));
! 401:
! 402: case t_QUAD: z=cgetg(ly,ty);
! 403: copyifstack(y[1],z[1]);
! 404: z[2]=lmulsg(s,(GEN)y[2]);
! 405: z[3]=lmulsg(s,(GEN)y[3]); return z;
! 406:
! 407: case t_POLMOD: z=cgetg(ly,ty);
! 408: z[2]=lmulsg(s,(GEN)y[2]);
! 409: copyifstack(y[1],z[1]); return z;
! 410:
! 411: case t_POL:
! 412: if (!s || !signe(y)) return zeropol(varn(y));
! 413: ly=lgef(y); z=cgetg(ly,t_POL); z[1]=y[1];
! 414: for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
! 415: return normalizepol_i(z, ly);
! 416:
! 417: case t_SER:
! 418: if (!s) return zeropol(varn(y));
! 419: if (gcmp0(y)) return gcopy(y);
! 420: z=cgetg(ly,ty);
! 421: for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
! 422: z[1]=y[1]; return normalize(z);
! 423:
! 424: case t_RFRAC:
! 425: if (!s) return zeropol(gvar(y));
! 426: z = cgetg(3, t_RFRAC);
! 427: i = ggcd(stoi(s),(GEN)y[2])[2];
! 428: avma = (long)z;
! 429: if (i == 1)
! 430: {
! 431: z[1]=lmulgs((GEN)y[1], s);
! 432: z[2]= lcopy((GEN)y[2]);
! 433: }
! 434: else
! 435: {
! 436: z[1] = lmulgs((GEN)y[1], s/i);
! 437: z[2] = ldivgs((GEN)y[2], i);
! 438: }
! 439: return z;
! 440:
! 441: case t_RFRACN:
! 442: if (!s) return zeropol(gvar(y));
! 443: z=cgetg(ly,t_RFRACN);
! 444: z[1]=lmulsg(s,(GEN)y[1]);
! 445: z[2]=lcopy((GEN)y[2]); return z;
! 446:
! 447: case t_VEC: case t_COL: case t_MAT:
! 448: z=cgetg(ly,ty);
! 449: for (i=1; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
! 450: return z;
! 451: }
! 452: err(typeer,"gmulsg");
! 453: return NULL; /* not reached */
! 454: }
! 455:
! 456: /********************************************************************/
! 457: /** **/
! 458: /** DIVISION SIMPLE **/
! 459: /** **/
! 460: /********************************************************************/
! 461:
! 462: GEN
! 463: gdivgs(GEN x, long s)
! 464: {
! 465: static long court[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
! 466: long tx=typ(x),lx=lg(x),i,av;
! 467: GEN z,p1;
! 468:
! 469: if (!s) err(gdiver2);
! 470: switch(tx)
! 471: {
! 472: case t_INT:
! 473: av=avma; z=dvmdis(x,s,&p1);
! 474: if (p1==gzero) return z;
! 475:
! 476: i = cgcd(s, smodis(x,s));
! 477: avma=av; z=cgetg(3,t_FRAC);
! 478: if (i == 1)
! 479: x = icopy(x);
! 480: else
! 481: {
! 482: s /= i;
! 483: x = divis(x, i);
! 484: }
! 485: z[1]=(long)x; z[2]=lstoi(s);
! 486: fix_frac(z); return z;
! 487:
! 488: case t_REAL:
! 489: return divrs(x,s);
! 490:
! 491: case t_FRAC: z=cgetg(3,tx);
! 492: i = cgcd(s, smodis((GEN)x[1], s));
! 493: if (i == 1)
! 494: {
! 495: z[2] = lmulsi(s, (GEN)x[2]);
! 496: z[1] = licopy((GEN)x[1]);
! 497: }
! 498: else
! 499: {
! 500: z[2] = lmulsi(s/i, (GEN)x[2]);
! 501: z[1] = ldivis((GEN)x[1], i);
! 502: }
! 503: fix_frac(z);
! 504: fix_frac_if_int(z); return z;
! 505:
! 506: case t_FRACN: z = cgetg(3,tx);
! 507: z[1]=licopy((GEN)x[1]);
! 508: z[2]=lmulsi(s,(GEN)x[2]);
! 509: fix_frac(z); return z;
! 510:
! 511: case t_COMPLEX: z=cgetg(lx,tx);
! 512: z[1]=ldivgs((GEN)x[1],s);
! 513: z[2]=ldivgs((GEN)x[2],s);
! 514: return z;
! 515:
! 516: case t_QUAD: z=cgetg(lx,tx);
! 517: copyifstack(x[1],z[1]);
! 518: for (i=2; i<4; i++) z[i]=ldivgs((GEN)x[i],s);
! 519: return z;
! 520:
! 521: case t_POLMOD: z=cgetg(lx,tx);
! 522: copyifstack(x[1],z[1]);
! 523: z[2]=ldivgs((GEN)x[2],s);
! 524: return z;
! 525:
! 526: case t_POL: lx=lgef(x); z=cgetg(lx,tx);
! 527: for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
! 528: z[1]=x[1]; return z;
! 529:
! 530: case t_SER: z=cgetg(lx,tx);
! 531: for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
! 532: z[1]=x[1]; return z;
! 533:
! 534: case t_RFRAC:
! 535: z = cgetg(3, t_RFRAC);
! 536: i = ggcd(stoi(s),(GEN)x[1])[2];
! 537: avma = (long)z;
! 538: if (i == 1)
! 539: {
! 540: z[2]=lmulsg(s,(GEN)x[2]);
! 541: z[1]=lcopy((GEN)x[1]); return z;
! 542: }
! 543: z[1] = ldivgs((GEN)x[1], i);
! 544: z[2] = lmulgs((GEN)x[2], s/i); return z;
! 545:
! 546: case t_RFRACN: z=cgetg(3,t_RFRACN);
! 547: z[2]=lmulsg(s,(GEN)x[2]);
! 548: z[1]=lcopy((GEN)x[1]); return z;
! 549:
! 550: case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
! 551: for (i=1; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
! 552: return z;
! 553: }
! 554: affsi(s,court); return gdiv(x,court);
! 555: }
! 556:
! 557: /*******************************************************************/
! 558: /* */
! 559: /* MODULO GENERAL */
! 560: /* */
! 561: /*******************************************************************/
! 562:
! 563: GEN
! 564: gmod(GEN x, GEN y)
! 565: {
! 566: long av,tetpil,i, tx=typ(x);
! 567: GEN z,p1;
! 568:
! 569: if (is_matvec_t(tx))
! 570: {
! 571: i=lg(x); z=cgetg(i,tx);
! 572: for (i--; i; i--) z[i]=lmod((GEN)x[i],y);
! 573: return z;
! 574: }
! 575: switch(typ(y))
! 576: {
! 577: case t_INT:
! 578: switch(tx)
! 579: {
! 580: case t_INT:
! 581: return modii(x,y);
! 582:
! 583: case t_INTMOD: z=cgetg(3,tx);
! 584: z[1]=lmppgcd((GEN)x[1],y);
! 585: z[2]=lmodii((GEN)x[2],(GEN)z[1]); return z;
! 586:
! 587: case t_FRAC: case t_FRACN:
! 588: av=avma; if (tx==t_FRACN) x=gred(x);
! 589: p1=mulii((GEN)x[1],mpinvmod((GEN)x[2],y));
! 590: tetpil=avma; return gerepile(av,tetpil,modii(p1,y));
! 591:
! 592: case t_QUAD: z=cgetg(4,tx);
! 593: copyifstack(x[1],z[1]);
! 594: z[2]=lmod((GEN)x[2],y);
! 595: z[3]=lmod((GEN)x[3],y); return z;
! 596:
! 597: case t_PADIC:
! 598: {
! 599: long p1[] = {evaltyp(t_INTMOD)|m_evallg(3),0,0};
! 600: p1[1] = (long)y;
! 601: p1[2] = lgeti(lgefint(y));
! 602: gaffect(x,p1); return (GEN)p1[2];
! 603: }
! 604: case t_POLMOD: case t_POL:
! 605: return gzero;
! 606:
! 607: default: err(gmoder1);
! 608: }
! 609:
! 610: case t_REAL: case t_FRAC: case t_FRACN:
! 611: switch(tx)
! 612: {
! 613: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 614: av=avma; p1 = gfloor(gdiv(x,y)); p1 = gneg_i(gmul(p1,y));
! 615: tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
! 616:
! 617: case t_POLMOD: case t_POL:
! 618: return gzero;
! 619:
! 620: default: err(gmoder1);
! 621: }
! 622:
! 623: case t_POL:
! 624: if (is_scalar_t(tx))
! 625: {
! 626: if (tx!=t_POLMOD || varn(x[1]) > varn(y))
! 627: return (lgef(y) < 4)? gzero: gcopy(x);
! 628: if (varn(x[1])!=varn(y)) return gzero;
! 629: z=cgetg(3,t_POLMOD);
! 630: z[1]=lgcd((GEN)x[1],y);
! 631: z[2]=lres((GEN)x[2],(GEN)z[1]); return z;
! 632: }
! 633: switch(tx)
! 634: {
! 635: case t_POL:
! 636: return gres(x,y);
! 637:
! 638: case t_RFRAC: case t_RFRACN:
! 639: av=avma; if (tx==t_RFRACN) x=gred_rfrac(x);
! 640: p1=gmul((GEN)x[1],ginvmod((GEN)x[2],y)); tetpil=avma;
! 641: return gerepile(av,tetpil,gres(p1,y));
! 642:
! 643: default: err(talker,"type mod polynomial forbidden in gmod");
! 644: }
! 645: }
! 646: err(talker,"modulus type forbidden in gmod");
! 647: return NULL; /* not reached */
! 648: }
! 649:
! 650: GEN
! 651: gmodulsg(long x, GEN y)
! 652: {
! 653: GEN z;
! 654:
! 655: switch(typ(y))
! 656: {
! 657: case t_INT: z=cgetg(3,t_INTMOD);
! 658: z[1]=labsi(y); z[2]=lmodsi(x,y); return z;
! 659:
! 660: case t_POL: z=cgetg(3,t_POLMOD);
! 661: z[1]=lcopy(y); z[2]=lstoi(x); return z;
! 662: }
! 663: err(gmoder1); return NULL; /* not reached */
! 664: }
! 665:
! 666: GEN
! 667: gmodulss(long x, long y)
! 668: {
! 669: GEN z=cgetg(3,t_INTMOD);
! 670:
! 671: y=labs(y); z[1]=lstoi(y); z[2]=lstoi(x % y); return z;
! 672: }
! 673:
! 674: GEN
! 675: gmodulo(GEN x,GEN y)
! 676: {
! 677: long tx=typ(x),l,i;
! 678: GEN z;
! 679:
! 680: if (is_matvec_t(tx))
! 681: {
! 682: l=lg(x); z=cgetg(l,tx);
! 683: for (i=1; i<l; i++) z[i] = (long) gmodulo((GEN)x[i],y);
! 684: return z;
! 685: }
! 686: switch(typ(y))
! 687: {
! 688: case t_INT:
! 689: if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
! 690: z=cgetg(3,t_INTMOD);
! 691: if (!signe(y)) err(talker,"zero modulus in gmodulo");
! 692: y = gclone(y); setsigne(y,1);
! 693: z[1]=(long)y;
! 694: z[2]=lmod(x,y); return z;
! 695:
! 696: case t_POL: z=cgetg(3,t_POLMOD);
! 697: z[1] = lclone(y);
! 698: if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
! 699: if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
! 700: }
! 701: err(gmoder1); return NULL; /* not reached */
! 702: }
! 703:
! 704:
! 705: GEN
! 706: gmodulcp(GEN x,GEN y)
! 707: {
! 708: long tx=typ(x),l,i;
! 709: GEN z;
! 710:
! 711: if (is_matvec_t(tx))
! 712: {
! 713: l=lg(x); z=cgetg(l,tx);
! 714: for (i=1; i<l; i++) z[i] = lmodulcp((GEN)x[i],y);
! 715: return z;
! 716: }
! 717: switch(typ(y))
! 718: {
! 719: case t_INT:
! 720: if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
! 721: z=cgetg(3,t_INTMOD);
! 722: z[1]=labsi(y);
! 723: z[2]=lmod(x,y); return z;
! 724:
! 725: case t_POL: z=cgetg(3,t_POLMOD);
! 726: z[1]=lcopy(y);
! 727: if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
! 728: if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
! 729: }
! 730: err(gmoder1); return NULL; /* not reached */
! 731: }
! 732:
! 733: GEN
! 734: Mod0(GEN x,GEN y,long flag)
! 735: {
! 736: switch(flag)
! 737: {
! 738: case 0: return gmodulcp(x,y);
! 739: case 1: return gmodulo(x,y);
! 740: default: err(flagerr,"Mod");
! 741: }
! 742: return NULL; /* not reached */
! 743: }
! 744:
! 745: /*******************************************************************/
! 746: /* */
! 747: /* DIVISION ENTIERE GENERALE */
! 748: /* DIVISION ENTIERE AVEC RESTE GENERALE */
! 749: /* */
! 750: /*******************************************************************/
! 751:
! 752: GEN
! 753: gdivent(GEN x, GEN y)
! 754: {
! 755: long tx=typ(x),ty=typ(y);
! 756:
! 757: if (tx == t_INT)
! 758: {
! 759: if (ty == t_INT)
! 760: {
! 761: GEN z,p1;
! 762: long av,tetpil;
! 763:
! 764: if (signe(x)>=0) return divii(x,y);
! 765: av=avma; z=dvmdii(x,y,&p1);
! 766: if (p1==gzero) return z;
! 767: tetpil=avma; return gerepile(av,tetpil,addsi(-signe(y),z));
! 768: }
! 769: if (ty!=t_POL) err(typeer,"gdivent");
! 770: return gzero;
! 771: }
! 772: if (ty != t_POL) err(typeer,"gdivent");
! 773: if (tx == t_POL) return gdeuc(x,y);
! 774: if (! is_scalar_t(tx)) err(typeer,"gdivent");
! 775: return gzero;
! 776: }
! 777:
! 778: GEN
! 779: gdiventres(GEN x, GEN y)
! 780: {
! 781: long tx=typ(x),ty=typ(y);
! 782: GEN z = cgetg(3,t_COL);
! 783:
! 784: if (tx==t_INT)
! 785: {
! 786: if (ty==t_INT) z[1]=(long)truedvmdii(x,y,(GEN*)(z+2));
! 787: else
! 788: {
! 789: if (ty!=t_POL) err(typeer,"gdiventres");
! 790: z[1]=zero; z[2]=licopy(x);
! 791: }
! 792: return z;
! 793: }
! 794: if (ty != t_POL) err(typeer,"gdiventres");
! 795: if (tx == t_POL)
! 796: {
! 797: z[1]=ldivres(x,y,(GEN *)(z+2));
! 798: return z;
! 799: }
! 800: if (! is_scalar_t(tx)) err(typeer,"gdiventres");
! 801: z[1]=zero; z[2]=lcopy(x); return z;
! 802: }
! 803:
! 804: GEN
! 805: gdivmod(GEN x, GEN y, GEN *pr)
! 806: {
! 807: long ty,tx=typ(x);
! 808:
! 809: if (tx==t_INT)
! 810: {
! 811: ty=typ(y);
! 812: if (ty==t_INT) return dvmdii(x,y,pr);
! 813: if (ty==t_POL) { *pr=gcopy(x); return gzero; }
! 814: err(typeer,"gdivmod");
! 815: }
! 816: if (tx!=t_POL) err(typeer,"gdivmod");
! 817: return poldivres(x,y,pr);
! 818: }
! 819:
! 820: /* When x and y are integers, compute the quotient x/y, rounded to the
! 821: * nearest integer. If there is a tie, the quotient is rounded towards zero.
! 822: * If x and y are not both integers, same as gdivent.
! 823: */
! 824: GEN
! 825: gdivround(GEN x, GEN y)
! 826: {
! 827: long tx=typ(x),ty=typ(y);
! 828:
! 829: if (tx==t_INT)
! 830: {
! 831: if (ty==t_INT)
! 832: {
! 833: long av = avma, av1,fl;
! 834: GEN r, q=dvmdii(x,y,&r); /* q = x/y rounded towards 0, sqn(r)=sgn(x) */
! 835:
! 836: if (r==gzero) return q;
! 837: av1 = avma;
! 838: fl = absi_cmp(shifti(r,1),y);
! 839: avma = av1; cgiv(r);
! 840: if (fl >= 0) /* If 2*|r| >= |y| */
! 841: {
! 842: long sz = signe(x)*signe(y);
! 843: if (fl || sz > 0)
! 844: { av1=avma; q = gerepile(av,av1,addis(q,sz)); }
! 845: }
! 846: return q;
! 847: }
! 848: if (ty!=t_POL) err(typeer,"gdivround");
! 849: return gzero;
! 850: }
! 851: if (ty != t_POL) err(typeer,"gdivround");
! 852: if (tx == t_POL) return gdeuc(x,y);
! 853: if (! is_scalar_t(tx)) err(typeer,"gdivround");
! 854: return gzero;
! 855: }
! 856:
! 857: /*******************************************************************/
! 858: /* */
! 859: /* SHIFT D'UN GEN */
! 860: /* */
! 861: /*******************************************************************/
! 862:
! 863: /* Shift tronque si n<0 (multiplication tronquee par 2^n) */
! 864: GEN
! 865: gshift(GEN x, long n)
! 866: {
! 867: long i,l,lx, tx = typ(x);
! 868: GEN y;
! 869:
! 870: switch(tx)
! 871: {
! 872: case t_INT:
! 873: return shifti(x,n);
! 874: case t_REAL:
! 875: return shiftr(x,n);
! 876:
! 877: case t_VEC: case t_COL: case t_MAT:
! 878: lx=lg(x); y=cgetg(lx,tx); l=lontyp[tx];
! 879: for (i=1; i<l ; i++) y[i]=x[i];
! 880: for ( ; i<lx; i++) y[i]=lshift((GEN)x[i],n);
! 881: return y;
! 882: }
! 883: return gmul2n(x,n);
! 884: }
! 885:
! 886: GEN mulscalrfrac(GEN x, GEN y);
! 887:
! 888: /* Shift vrai (multiplication exacte par 2^n) */
! 889: GEN
! 890: gmul2n(GEN x, long n)
! 891: {
! 892: long tx,lx,i,k,l,av,tetpil;
! 893: GEN p2,p1,y;
! 894:
! 895: tx=typ(x);
! 896: switch(tx)
! 897: {
! 898: case t_INT:
! 899: if (n>=0) return shifti(x,n);
! 900: if (!signe(x)) return gzero;
! 901: l = vali(x); n = -n;
! 902: if (n<=l) return shifti(x,-n);
! 903: y=cgetg(3,t_FRAC);
! 904: y[1]=lshifti(x,-l);
! 905: y[2]=lshifti(gun,n-l); return y;
! 906:
! 907: case t_REAL:
! 908: return shiftr(x,n);
! 909:
! 910: case t_INTMOD:
! 911: if (n > 0)
! 912: {
! 913: y=cgetg(3,t_INTMOD); p2=(GEN)x[1];
! 914: av=avma; new_chunk(lgefint(p2) * (3 + (n>>TWOPOTBITS_IN_LONG)));
! 915: p1 = shifti((GEN)x[2],n); avma=av;
! 916: y[2]=lmodii(p1,p2); icopyifstack(p2,y[1]); return y;
! 917: }
! 918: l=avma; y=gmul2n(gun,n); tetpil=avma;
! 919: return gerepile(l,tetpil,gmul(y,x));
! 920:
! 921: case t_FRAC: case t_FRACN:
! 922: l = vali((GEN)x[1]);
! 923: k = vali((GEN)x[2]);
! 924: if (n+l-k>=0)
! 925: {
! 926: if (expi((GEN)x[2]) == k) /* x[2] power of 2 */
! 927: return shifti((GEN)x[1],n-k);
! 928: l = n-k; k = -k;
! 929: }
! 930: else
! 931: {
! 932: k = -l-n; l = -l;
! 933: }
! 934: y=cgetg(3,t_FRAC);
! 935: y[1]=lshifti((GEN)x[1],l);
! 936: y[2]=lshifti((GEN)x[2],k); return y;
! 937:
! 938: case t_QUAD: y=cgetg(4,t_QUAD);
! 939: copyifstack(x[1],y[1]);
! 940: for (i=2; i<4; i++) y[i]=lmul2n((GEN)x[i],n);
! 941: return y;
! 942:
! 943: case t_POLMOD: y=cgetg(3,t_POLMOD);
! 944: copyifstack(x[1],y[1]);
! 945: y[2]=lmul2n((GEN)x[2],n); return y;
! 946:
! 947: case t_POL: case t_COMPLEX: case t_SER:
! 948: case t_VEC: case t_COL: case t_MAT:
! 949: lx = (tx==t_POL)? lgef(x): lg(x);
! 950: y=cgetg(lx,tx); l=lontyp[tx];
! 951: for (i=1; i<l ; i++) y[i]=x[i];
! 952: for ( ; i<lx; i++) y[i]=lmul2n((GEN)x[i],n);
! 953: return y;
! 954:
! 955: case t_RFRAC: av=avma; p1 = gmul2n(gun,n); tetpil = avma;
! 956: return gerepile(av,tetpil, mulscalrfrac(p1,x));
! 957:
! 958: case t_RFRACN: y=cgetg(3,tx);
! 959: if (n>=0)
! 960: {
! 961: y[1]=lmul2n((GEN)x[1],n); y[2]=lcopy((GEN)x[2]);
! 962: }
! 963: else
! 964: {
! 965: y[2]=lmul2n((GEN)x[2],-n); y[1]=lcopy((GEN)x[1]);
! 966: }
! 967: return y;
! 968:
! 969: case t_PADIC:
! 970: l=avma; y=gmul2n(gun,n); tetpil=avma;
! 971: return gerepile(l,tetpil,gmul(y,x));
! 972: }
! 973: err(typeer,"gmul2n");
! 974: return NULL; /* not reached */
! 975: }
! 976:
! 977: /*******************************************************************/
! 978: /* */
! 979: /* INVERSE D'UN GEN */
! 980: /* */
! 981: /*******************************************************************/
! 982: GEN fix_rfrac_if_pol(GEN x, GEN y);
! 983:
! 984: GEN
! 985: ginv(GEN x)
! 986: {
! 987: long tx=typ(x),av,tetpil,s;
! 988: GEN z,y,p1,p2;
! 989:
! 990: switch(tx)
! 991: {
! 992: case t_INT:
! 993: if (is_pm1(x)) return icopy(x);
! 994: if (!signe(x)) err(gdiver2);
! 995: z=cgetg(3,t_FRAC);
! 996: z[1] = (signe(x)<0)? lnegi(gun): un;
! 997: z[2]=labsi(x); return z;
! 998:
! 999: case t_REAL:
! 1000: return divsr(1,x);
! 1001:
! 1002: case t_INTMOD: z=cgetg(3,t_INTMOD);
! 1003: icopyifstack(x[1],z[1]);
! 1004: z[2]=lmpinvmod((GEN)x[2],(GEN)x[1]); return z;
! 1005:
! 1006: case t_FRAC: case t_FRACN: z=cgetg(3,tx);
! 1007: s = signe(x[1]);
! 1008: if (!s) err(gdiver2);
! 1009: if (is_pm1(x[1]))
! 1010: return s>0? icopy((GEN)x[2]): negi((GEN)x[2]);
! 1011: z[1]=licopy((GEN)x[2]);
! 1012: z[2]=licopy((GEN)x[1]);
! 1013: if (s < 0)
! 1014: {
! 1015: setsigne(z[1],-signe(z[1]));
! 1016: setsigne(z[2],1);
! 1017: }
! 1018: return z;
! 1019:
! 1020: case t_COMPLEX: case t_QUAD:
! 1021: av=avma; p1=gnorm(x); p2=gconj(x); tetpil=avma;
! 1022: return gerepile(av,tetpil,gdiv(p2,p1));
! 1023:
! 1024: case t_PADIC: z = cgetg(5,t_PADIC);
! 1025: if (!signe(x[4])) err(gdiver2);
! 1026: z[1] = evalprecp(precp(x)) | evalvalp(-valp(x));
! 1027: icopyifstack(x[2], z[2]);
! 1028: z[3] = licopy((GEN)x[3]);
! 1029: z[4] = lmpinvmod((GEN)x[4],(GEN)z[3]); return z;
! 1030:
! 1031: case t_POLMOD: z=cgetg(3,t_POLMOD);
! 1032: copyifstack(x[1],z[1]);
! 1033: z[2]=linvmod((GEN)x[2],(GEN)x[1]); return z;
! 1034:
! 1035: case t_POL: case t_SER:
! 1036: return gdiv(gun,x);
! 1037:
! 1038: case t_RFRAC: case t_RFRACN:
! 1039: if (gcmp0((GEN) x[1])) err(gdiver2);
! 1040: p1 = fix_rfrac_if_pol((GEN)x[2],(GEN)x[1]);
! 1041: if (p1) return p1;
! 1042: z=cgetg(3,tx);
! 1043: z[1]=lcopy((GEN)x[2]);
! 1044: z[2]=lcopy((GEN)x[1]); return z;
! 1045:
! 1046: case t_QFR:
! 1047: {
! 1048: long k,l;
! 1049: l=signe(x[2]); setsigne(x[2],-l);
! 1050: k=signe(x[4]); setsigne(x[4],-k); z=redreal(x);
! 1051: setsigne(x[2],l); setsigne(x[4],k); return z;
! 1052: }
! 1053: case t_QFI:
! 1054: y=gcopy(x);
! 1055: if (!egalii((GEN)x[1],(GEN)x[2]) && !egalii((GEN)x[1],(GEN)x[3]))
! 1056: setsigne(y[2],-signe(y[2]));
! 1057: return y;
! 1058: case t_MAT:
! 1059: return (lg(x)==1)? cgetg(1,t_MAT): invmat(x);
! 1060: }
! 1061: err(typeer,"inverse");
! 1062: return NULL; /* not reached */
! 1063: }
! 1064:
! 1065: /*******************************************************************/
! 1066: /* */
! 1067: /* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
! 1068: /* */
! 1069: /*******************************************************************/
! 1070:
! 1071: /* Convert t_SER --> t_POL / t_RFRAC */
! 1072: static GEN
! 1073: gconvsp(GEN x, int flpile)
! 1074: {
! 1075: long v = varn(x), av,tetpil,i;
! 1076: GEN p1,y;
! 1077:
! 1078: if (gcmp0(x)) return zeropol(v);
! 1079: av=avma; y=dummycopy(x); settyp(y,t_POL);
! 1080: i=lg(x)-1; while (i>1 && gcmp0((GEN)y[i])) i--;
! 1081: setlgef(y,i+1);
! 1082: p1=gpuigs(polx[v],valp(x));
! 1083: tetpil=avma; p1=gmul(p1,y);
! 1084: return flpile? gerepile(av,tetpil,p1): p1;
! 1085: }
! 1086:
! 1087: GEN
! 1088: gsubst(GEN x, long v, GEN y)
! 1089: {
! 1090: long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
! 1091: long l,vx,vy,e,ex,ey,tetpil,av,i,j,k,jb,limite;
! 1092: GEN t,p1,p2,z;
! 1093:
! 1094: if (ty==t_MAT)
! 1095: {
! 1096: if (ly==1) return cgetg(1,t_MAT);
! 1097: if (ly != lg(y[1]))
! 1098: err(talker,"forbidden substitution by a non square matrix");
! 1099: } else if (is_graphicvec_t(ty))
! 1100: err(talker,"forbidden substitution by a vector");
! 1101:
! 1102: if (is_scalar_t(tx))
! 1103: {
! 1104: if (tx!=t_POLMOD || v <= varn(x[1]))
! 1105: {
! 1106: if (ty==t_MAT) return gscalmat(x,ly-1);
! 1107: return gcopy(x);
! 1108: }
! 1109: av=avma;
! 1110: p1=gsubst((GEN)x[1],v,y); vx=varn(p1);
! 1111: p2=gsubst((GEN)x[2],v,y); vy=gvar(p2);
! 1112: if (typ(p1)!=t_POL)
! 1113: err(talker,"forbidden substitution in a scalar type");
! 1114: if (vy>=vx)
! 1115: {
! 1116: tetpil=avma;
! 1117: return gerepile(av,tetpil,gmodulcp(p2,p1));
! 1118: }
! 1119: lx=lgef(p2); tetpil=avma; z=cgetg(lx,t_POL); z[1]=p2[1];
! 1120: for (i=2; i<lx; i++) z[i]=lmodulcp((GEN)p2[i],p1);
! 1121: return gerepile(av,tetpil, normalizepol_i(z,lx));
! 1122: }
! 1123:
! 1124: switch(tx)
! 1125: {
! 1126: case t_POL:
! 1127: l=lgef(x);
! 1128: if (l==2)
! 1129: return (ty==t_MAT)? gscalmat(gzero,ly-1): gzero;
! 1130:
! 1131: vx=varn(x);
! 1132: if (vx<v)
! 1133: {
! 1134: av=avma; p1=polx[vx]; z= gsubst((GEN)x[l-1],v,y);
! 1135: for (i=l-1; i>2; i--) z=gadd(gmul(z,p1),gsubst((GEN)x[i-1],v,y));
! 1136: return gerepileupto(av,z);
! 1137: }
! 1138: if (ty!=t_MAT)
! 1139: return (vx>v)? gcopy(x): poleval(x,y);
! 1140:
! 1141: if (vx>v) return gscalmat(x,ly-1);
! 1142: if (l==3) return gscalmat((GEN)x[2],ly-1);
! 1143: av=avma; z=(GEN)x[l-1];
! 1144: for (i=l-1; i>2; i--) z=gaddmat((GEN)x[i-1],gmul(z,y));
! 1145: return gerepileupto(av,z);
! 1146:
! 1147: case t_SER:
! 1148: vx=varn(x);
! 1149: if (vx > v)
! 1150: return (ty==t_MAT)? gscalmat(x,ly-1): gcopy(x);
! 1151: ex = valp(x);
! 1152: if (vx < v)
! 1153: {
! 1154: if (!signe(x)) return gcopy(x);
! 1155: /* a ameliorer */
! 1156: av=avma; setvalp(x,0); p1=gconvsp(x,0); setvalp(x,ex);
! 1157: p2=gsubst(p1,v,y); tetpil=avma; z=tayl(p2,vx,lx-2);
! 1158: if (ex)
! 1159: {
! 1160: p1=gpuigs(polx[vx],ex); tetpil=avma; z=gmul(z,p1);
! 1161: }
! 1162: return gerepile(av,tetpil,z);
! 1163: }
! 1164: switch(ty) /* here vx == v */
! 1165: {
! 1166: case t_SER:
! 1167: ey=valp(y); vy=varn(y);
! 1168: if (ey<1) return zeroser(vy,ey*(ex+lx-2));
! 1169: l=(lx-2)*ey+2;
! 1170: if (ex) { if (l>ly) l=ly; }
! 1171: else
! 1172: {
! 1173: if (gcmp0(y)) l=ey+2;
! 1174: else { if (l>ly) l=ey+ly; }
! 1175: }
! 1176: if (vy!=vx)
! 1177: {
! 1178: av=avma; z = zeroser(vy,0);
! 1179: for (i=lx-1; i>=2; i--)
! 1180: z = gadd((GEN)x[i], gmul(y,z));
! 1181: if (ex) z = gmul(z, gpuigs(y,ex));
! 1182: return gerepileupto(av,z);
! 1183: }
! 1184:
! 1185: av=avma; limite=stack_lim(av,1);
! 1186: t=cgetg(ly,t_SER);
! 1187: z = scalarser((GEN)x[2],varn(y),l-2);
! 1188: for (i=2; i<ly; i++) t[i]=y[i];
! 1189:
! 1190: for (i=3,jb=ey; jb<=l-2; i++,jb+=ey)
! 1191: {
! 1192: for (j=jb+2; j<l; j++)
! 1193: {
! 1194: p1=gmul((GEN)x[i],(GEN)t[j-jb]);
! 1195: z[j]=ladd((GEN)z[j],p1);
! 1196: }
! 1197: for (j=l-1-jb-ey; j>1; j--)
! 1198: {
! 1199: p1=gzero;
! 1200: for (k=2; k<j; k++)
! 1201: p1=gadd(p1,gmul((GEN)t[j-k+2],(GEN)y[k]));
! 1202: p2=gmul((GEN)t[2],(GEN)y[j]);
! 1203: t[j]=ladd(p1,p2);
! 1204: }
! 1205: if (low_stack(limite, stack_lim(av,1)))
! 1206: {
! 1207: GEN *gptr[2];
! 1208: if(DEBUGMEM>1) err(warnmem,"gsubst");
! 1209: gptr[0]=&z; gptr[1]=&t; gerepilemany(av,gptr,2);
! 1210: }
! 1211: }
! 1212: if (!ex) { tetpil=avma; return gerepile(av,tetpil,gcopy(z)); }
! 1213:
! 1214: if (l<ly) { setlg(y,l); p1=gpuigs(y,ex); setlg(y,ly); }
! 1215: else p1=gpuigs(y,ex);
! 1216: tetpil=avma; return gerepile(av,tetpil,gmul(z,p1));
! 1217:
! 1218: case t_POL: case t_RFRAC: case t_RFRACN:
! 1219: if (isexactzero(y)) return scalarser((GEN)x[2],v,lx-2);
! 1220: vy=gvar(y); e=gval(y,vy);
! 1221: if (e<=0)
! 1222: err(talker,"non positive valuation in a series substitution");
! 1223: av=avma; p1=gconvsp(x,0); p2=gsubst(p1,v,y); tetpil=avma;
! 1224: return gerepile(av,tetpil,tayl(p2,vy,e*(lx-2+ex)));
! 1225:
! 1226: default:
! 1227: err(talker,"non polynomial or series type substituted in a series");
! 1228: }
! 1229: break;
! 1230:
! 1231: case t_RFRAC: case t_RFRACN: av=avma;
! 1232: p1=gsubst((GEN)x[1],v,y);
! 1233: p2=gsubst((GEN)x[2],v,y); tetpil=avma;
! 1234: return gerepile(av,tetpil,gdiv(p1,p2));
! 1235:
! 1236: case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
! 1237: for (i=1; i<lx; i++) z[i]=lsubst((GEN)x[i],v,y);
! 1238: }
! 1239: return z;
! 1240: }
! 1241:
! 1242: /*******************************************************************/
! 1243: /* */
! 1244: /* SERIE RECIPROQUE D'UNE SERIE */
! 1245: /* */
! 1246: /*******************************************************************/
! 1247:
! 1248: GEN
! 1249: recip(GEN x)
! 1250: {
! 1251: long tetpil, av=avma, v=varn(x);
! 1252: GEN p1,p2,a,y,u;
! 1253:
! 1254: if (typ(x)!=t_SER) err(talker,"not a series in serreverse");
! 1255: if (valp(x)!=1) err(talker,"valuation not equal to 1 in serreverse");
! 1256:
! 1257: a=(GEN)x[2];
! 1258: if (gcmp1(a))
! 1259: {
! 1260: long i,j,k, lx=lg(x), lim=stack_lim(av,2);
! 1261:
! 1262: u=cgetg(lx,t_SER); y=cgetg(lx,t_SER);
! 1263: u[1]=y[1]=evalsigne(1) | evalvalp(1) | evalvarn(v);
! 1264: u[2]=un; u[3]=lmulsg(-2,(GEN)x[3]);
! 1265: y[2]=un; y[3]=lneg((GEN)x[3]);
! 1266: for (i=3; i<lx-1; )
! 1267: {
! 1268: for (j=3; j<i+1; j++)
! 1269: {
! 1270: p1=(GEN)u[j];
! 1271: for (k=j-1; k>2; k--)
! 1272: p1=gsub(p1,gmul((GEN)u[k],(GEN)x[j-k+2]));
! 1273: u[j]=lsub(p1,(GEN)x[j]);
! 1274: }
! 1275: p1=gmulsg(i,(GEN)x[i+1]);
! 1276: for (k=2; k<i; k++)
! 1277: {
! 1278: p2=gmul((GEN)x[k+1],(GEN)u[i-k+2]);
! 1279: p1=gadd(p1,gmulsg(k,p2));
! 1280: }
! 1281: i++; u[i]=lneg(p1); y[i]=ldivgs((GEN)u[i],i-1);
! 1282: if (low_stack(lim, stack_lim(av,2)))
! 1283: {
! 1284: GEN *gptr[2];
! 1285: if(DEBUGMEM>1) err(warnmem,"recip");
! 1286: for(k=i+1; k<lx; k++) u[k]=y[k]=zero; /* dummy */
! 1287: gptr[0]=&u; gptr[1]=&y; gerepilemany(av,gptr,2);
! 1288: }
! 1289: }
! 1290: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
! 1291: }
! 1292: y=gdiv(x,a); y[2]=un; y=recip(y);
! 1293: a=gdiv(polx[v],a); tetpil=avma;
! 1294: return gerepile(av,tetpil,gsubst(y,v,a));
! 1295: }
! 1296:
! 1297: /*******************************************************************/
! 1298: /* */
! 1299: /* DERIVATION ET INTEGRATION */
! 1300: /* */
! 1301: /*******************************************************************/
! 1302: GEN
! 1303: derivpol(GEN x)
! 1304: {
! 1305: long i,lx = lgef(x)-1;
! 1306: GEN y;
! 1307:
! 1308: if (lx<3) return gzero;
! 1309: y = cgetg(lx,t_POL);
! 1310: for (i=2; i<lx ; i++) y[i] = lmulsg(i-1,(GEN)x[i+1]);
! 1311: y[1] = x[1]; return normalizepol_i(y,i);
! 1312: }
! 1313:
! 1314: GEN
! 1315: derivser(GEN x)
! 1316: {
! 1317: long i,j,vx = varn(x), e = valp(x), lx = lg(x);
! 1318: GEN y;
! 1319: if (gcmp0(x)) return zeroser(vx,e-1);
! 1320: if (e)
! 1321: {
! 1322: y=cgetg(lx,t_SER); y[1] = evalsigne(1) | evalvalp(e-1) | evalvarn(vx);
! 1323: for (i=2; i<lx; i++) y[i]=lmulsg(i+e-2,(GEN)x[i]);
! 1324: return y;
! 1325: }
! 1326: i=3; while (i<lx && gcmp0((GEN)x[i])) i++;
! 1327: if (i==lx) return zeroser(vx,lx-3);
! 1328: lx--; if (lx<3) lx=3;
! 1329: lx = lx-i+3; y=cgetg(lx,t_SER);
! 1330: y[1]=evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
! 1331: for (j=2; j<lx; j++) y[j]=lmulsg(j+i-4,(GEN)x[i+j-2]);
! 1332: return y;
! 1333: }
! 1334:
! 1335: GEN
! 1336: deriv(GEN x, long v)
! 1337: {
! 1338: long lx,vx,tx,e,i,j,l,av,tetpil;
! 1339: GEN y,p1,p2;
! 1340:
! 1341: tx=typ(x); if (is_const_t(tx)) return gzero;
! 1342: if (v < 0) v = gvar(x);
! 1343: switch(tx)
! 1344: {
! 1345: case t_POLMOD:
! 1346: if (v<=varn(x[1])) return gzero;
! 1347: y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
! 1348: y[2]=lderiv((GEN)x[2],v); return y;
! 1349:
! 1350: case t_POL:
! 1351: vx=varn(x); if (vx>v) return gzero;
! 1352: if (vx<v)
! 1353: {
! 1354: lx = lgef(x);
! 1355: y = cgetg(lx,t_POL);
! 1356: for (i=2; i<lx; i++) y[i] = lderiv((GEN)x[i],v);
! 1357: y[1] = evalvarn(vx);
! 1358: return normalizepol_i(y,i);
! 1359: }
! 1360: return derivpol(x);
! 1361:
! 1362: case t_SER:
! 1363: vx=varn(x); if (vx>v) return gzero;
! 1364: if (vx<v)
! 1365: {
! 1366: if (!signe(x)) return gcopy(x);
! 1367: lx=lg(x); e=valp(x);
! 1368: l=avma;
! 1369: for (i=2; i<lx; i++)
! 1370: {
! 1371: if (!gcmp0(deriv((GEN)x[i],v))) break;
! 1372: avma=l;
! 1373: }
! 1374: if (i==lx) return ggrando(polx[vx],e+lx-2);
! 1375: y=cgetg(lx-i+2,t_SER);
! 1376: y[1] = evalsigne(1) | evalvalp(e+i-2) | evalvarn(vx);
! 1377: for (j=2; i<lx; j++,i++) y[j]=lderiv((GEN)x[i],v);
! 1378: return y;
! 1379: }
! 1380: return derivser(x);
! 1381:
! 1382: case t_RFRAC: case t_RFRACN: av=avma; y=cgetg(3,tx);
! 1383: y[2]=lsqr((GEN)x[2]); l=avma;
! 1384: p1=gmul((GEN)x[2],deriv((GEN)x[1],v));
! 1385: p2=gmul(gneg_i((GEN)x[1]),deriv((GEN)x[2],v));
! 1386: tetpil=avma; p1=gadd(p1,p2);
! 1387: if (tx==t_RFRACN) { y[1]=lpile(l,tetpil,p1); return y; }
! 1388: y[1]=(long)p1; return gerepileupto(av,gred_rfrac(y));
! 1389:
! 1390: case t_VEC: case t_COL: case t_MAT: lx=lg(x); y=cgetg(lx,tx);
! 1391: for (i=1; i<lx; i++) y[i]=lderiv((GEN)x[i],v);
! 1392: return y;
! 1393: }
! 1394: err(typeer,"deriv");
! 1395: return NULL; /* not reached */
! 1396: }
! 1397:
! 1398: /*******************************************************************/
! 1399: /* */
! 1400: /* INTEGRATION FORMELLE */
! 1401: /* */
! 1402: /*******************************************************************/
! 1403:
! 1404: static GEN
! 1405: triv_integ(GEN x, long v, long tx, long lx)
! 1406: {
! 1407: GEN y=cgetg(lx,tx);
! 1408: long i;
! 1409:
! 1410: y[1]=x[1];
! 1411: for (i=2; i<lx; i++) y[i]=linteg((GEN)x[i],v);
! 1412: return y;
! 1413: }
! 1414:
! 1415: GEN
! 1416: integ(GEN x, long v)
! 1417: {
! 1418: long lx,tx,e,i,j,vx,n,av=avma,tetpil;
! 1419: GEN y,p1;
! 1420:
! 1421: tx = typ(x);
! 1422: if (v < 0) v = gvar(x);
! 1423: if (is_scalar_t(tx))
! 1424: {
! 1425: if (tx == t_POLMOD && v>varn(x[1]))
! 1426: {
! 1427: y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
! 1428: y[2]=linteg((GEN)x[2],v); return y;
! 1429: }
! 1430: if (gcmp0(x)) return gzero;
! 1431:
! 1432: y=cgetg(4,t_POL); y[1] = evalsigne(1) | evallgef(4) | evalvarn(v);
! 1433: y[2]=zero; y[3]=lcopy(x); return y;
! 1434: }
! 1435:
! 1436: switch(tx)
! 1437: {
! 1438: case t_POL:
! 1439: vx=varn(x); lx=lgef(x);
! 1440: if (lx==2) return zeropol(min(v,vx));
! 1441: if (vx>v)
! 1442: {
! 1443: y=cgetg(4,t_POL);
! 1444: y[1] = signe(x)? evallgef(4) | evalvarn(v) | evalsigne(1)
! 1445: : evallgef(4) | evalvarn(v);
! 1446: y[2]=zero; y[3]=lcopy(x); return y;
! 1447: }
! 1448: if (vx<v) return triv_integ(x,v,tx,lx);
! 1449: y=cgetg(lx+1,tx); y[2]=zero;
! 1450: for (i=3; i<=lx; i++)
! 1451: y[i]=ldivgs((GEN)x[i-1],i-2);
! 1452: y[1] = signe(x)? evallgef(lx+1) | evalvarn(v) | evalsigne(1)
! 1453: : evallgef(lx+1) | evalvarn(v);
! 1454: return y;
! 1455:
! 1456: case t_SER:
! 1457: lx=lg(x); e=valp(x); vx=varn(x);
! 1458: if (!signe(x))
! 1459: {
! 1460: if (vx == v) e++; else if (vx < v) v = vx;
! 1461: return zeroser(v,e);
! 1462: }
! 1463: if (vx>v)
! 1464: {
! 1465: y=cgetg(4,t_POL);
! 1466: y[1] = evallgef(4) | evalvarn(v) | evalsigne(1);
! 1467: y[2]=zero; y[3]=lcopy(x); return y;
! 1468: }
! 1469: if (vx<v) return triv_integ(x,v,tx,lx);
! 1470: y=cgetg(lx,tx);
! 1471: for (i=2; i<lx; i++)
! 1472: {
! 1473: j=i+e-1;
! 1474: if (!j)
! 1475: {
! 1476: if (!gcmp0((GEN)x[i])) err(inter2);
! 1477: y[i]=zero;
! 1478: }
! 1479: else y[i] = ldivgs((GEN)x[i],j);
! 1480: }
! 1481: y[1]=x[1]+1; return y;
! 1482:
! 1483: case t_RFRAC: case t_RFRACN:
! 1484: vx = gvar(x);
! 1485: if (vx>v)
! 1486: {
! 1487: y=cgetg(4,t_POL);
! 1488: y[1] = signe(x[1])? evallgef(4) | evalvarn(v) | evalsigne(1)
! 1489: : evallgef(4) | evalvarn(v);
! 1490: y[2]=zero; y[3]=lcopy(x); return y;
! 1491: }
! 1492: if (vx<v)
! 1493: {
! 1494: p1=cgetg(v+2,t_VEC);
! 1495: for (i=0; i<vx; i++) p1[i+1] = lpolx[i];
! 1496: for (i=vx+1; i<v; i++) p1[i+1] = lpolx[i];
! 1497: p1[v+1] = lpolx[vx]; p1[vx+1] = lpolx[v];
! 1498: y=integ(changevar(x,p1),vx); tetpil=avma;
! 1499: return gerepile(av,tetpil,changevar(y,p1));
! 1500: }
! 1501:
! 1502: n=lgef(x[1])+lgef(x[2])-4;
! 1503: y=gdiv(gtrunc(gmul((GEN)x[2], integ(tayl(x,v,n),v))), (GEN)x[2]);
! 1504: if (!gegal(deriv(y,v),x)) err(inter2);
! 1505: if (typ(y)==t_RFRAC && lgef(y[1]) == lgef(y[2]))
! 1506: {
! 1507: GEN p2;
! 1508: tx=typ(y[1]); p1=is_scalar_t(tx)? (GEN)y[1]: leading_term(y[1]);
! 1509: tx=typ(y[2]); p2=is_scalar_t(tx)? (GEN)y[2]: leading_term(y[2]);
! 1510: y=gsub(y, gdiv(p1,p2));
! 1511: }
! 1512: return gerepileupto(av,y);
! 1513:
! 1514: case t_VEC: case t_COL: case t_MAT:
! 1515: lx=lg(x); y=cgetg(lx,tx);
! 1516: for (i=1; i<lg(x); i++) y[i]=linteg((GEN)x[i],v);
! 1517: return y;
! 1518: }
! 1519: err(typeer,"integ");
! 1520: return NULL; /* not reached */
! 1521: }
! 1522:
! 1523: /*******************************************************************/
! 1524: /* */
! 1525: /* PARTIES ENTIERES */
! 1526: /* */
! 1527: /*******************************************************************/
! 1528:
! 1529: GEN
! 1530: gfloor(GEN x)
! 1531: {
! 1532: GEN y;
! 1533: long i,lx, tx = typ(x);
! 1534:
! 1535: switch(tx)
! 1536: {
! 1537: case t_INT: case t_POL:
! 1538: return gcopy(x);
! 1539:
! 1540: case t_REAL:
! 1541: return mpent(x);
! 1542:
! 1543: case t_FRAC: case t_FRACN:
! 1544: return truedvmdii((GEN)x[1],(GEN)x[2],NULL);
! 1545:
! 1546: case t_RFRAC: case t_RFRACN:
! 1547: return gdeuc((GEN)x[1],(GEN)x[2]);
! 1548:
! 1549: case t_VEC: case t_COL: case t_MAT:
! 1550: lx=lg(x); y=cgetg(lx,tx);
! 1551: for (i=1; i<lx; i++) y[i]=lfloor((GEN)x[i]);
! 1552: return y;
! 1553: }
! 1554: err(typeer,"gfloor");
! 1555: return NULL; /* not reached */
! 1556: }
! 1557:
! 1558: GEN
! 1559: gfrac(GEN x)
! 1560: {
! 1561: long av = avma,tetpil;
! 1562: GEN p1 = gneg_i(gfloor(x));
! 1563:
! 1564: tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
! 1565: }
! 1566:
! 1567: GEN
! 1568: gceil(GEN x)
! 1569: {
! 1570: GEN y,p1;
! 1571: long i,lx,tx=typ(x),av,tetpil;
! 1572:
! 1573: switch(tx)
! 1574: {
! 1575: case t_INT: case t_POL:
! 1576: return gcopy(x);
! 1577:
! 1578: case t_REAL:
! 1579: av=avma; y=mpent(x);
! 1580: if (!gegal(x,y))
! 1581: {
! 1582: tetpil=avma; return gerepile(av,tetpil,addsi(1,y));
! 1583: }
! 1584: return y;
! 1585:
! 1586: case t_FRAC: case t_FRACN:
! 1587: av=avma; y=dvmdii((GEN)x[1],(GEN)x[2],&p1);
! 1588: if (p1 != gzero && gsigne(x)>0)
! 1589: {
! 1590: cgiv(p1); tetpil=avma;
! 1591: return gerepile(av,tetpil,addsi(1,y));
! 1592: }
! 1593: return y;
! 1594:
! 1595: case t_RFRAC: case t_RFRACN:
! 1596: return gdeuc((GEN)x[1],(GEN)x[2]);
! 1597:
! 1598: case t_VEC: case t_COL: case t_MAT:
! 1599: lx=lg(x); y=cgetg(lx,tx);
! 1600: for (i=1; i<lx; i++) y[i]=lceil((GEN)x[i]);
! 1601: return y;
! 1602: }
! 1603: err(typeer,"gceil");
! 1604: return NULL; /* not reached */
! 1605: }
! 1606:
! 1607: GEN
! 1608: round0(GEN x, GEN *pte)
! 1609: {
! 1610: if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
! 1611: return ground(x);
! 1612: }
! 1613:
! 1614: GEN
! 1615: ground(GEN x)
! 1616: {
! 1617: GEN y,p1;
! 1618: long i,lx,tx=typ(x),av,tetpil;
! 1619:
! 1620: switch(tx)
! 1621: {
! 1622: case t_INT: case t_INTMOD: case t_QUAD:
! 1623: return gcopy(x);
! 1624:
! 1625: case t_REAL:
! 1626: {
! 1627: long ex, s = signe(x);
! 1628: if (s==0 || (ex=expo(x)) < -1) return gzero;
! 1629: if (ex < 0) return s>0? gun: negi(gun);
! 1630: av=avma; p1 = realun(3 + (ex>>TWOPOTBITS_IN_LONG));
! 1631: setexpo(p1,-1); /* p1 = 0.5 */
! 1632: p1 = addrr(x,p1); tetpil = avma;
! 1633: return gerepile(av,tetpil,mpent(p1));
! 1634: }
! 1635: case t_FRAC: case t_FRACN:
! 1636: av=avma; p1 = addii(shifti((GEN)x[2], -1), (GEN)x[1]);
! 1637: return gerepileuptoint(av, truedvmdii(p1, (GEN)x[2], NULL));
! 1638:
! 1639: case t_POLMOD: y=cgetg(3,t_POLMOD);
! 1640: copyifstack(x[1],y[1]);
! 1641: y[2]=lround((GEN)x[2]); return y;
! 1642:
! 1643: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
! 1644: case t_VEC: case t_COL: case t_MAT:
! 1645: lx = (tx==t_POL)? lgef(x): lg(x);
! 1646: y=cgetg(lx,tx);
! 1647: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
! 1648: for ( ; i<lx; i++) y[i]=lround((GEN)x[i]);
! 1649: if (tx==t_POL) return normalizepol_i(y, lx);
! 1650: if (tx==t_SER) return normalize(y);
! 1651: return y;
! 1652: }
! 1653: err(typeer,"ground");
! 1654: return NULL; /* not reached */
! 1655: }
! 1656:
! 1657: /* e = number of error bits on integral part */
! 1658: GEN
! 1659: grndtoi(GEN x, long *e)
! 1660: {
! 1661: GEN y,p1;
! 1662: long i,tx=typ(x), lx,av,ex,e1;
! 1663:
! 1664: *e = -HIGHEXPOBIT;
! 1665: switch(tx)
! 1666: {
! 1667: case t_INT: case t_INTMOD: case t_QUAD:
! 1668: case t_FRAC: case t_FRACN:
! 1669: return ground(x);
! 1670:
! 1671: case t_REAL:
! 1672: av=avma; p1=gadd(ghalf,x); ex=expo(p1);
! 1673: if (ex<0)
! 1674: {
! 1675: if (signe(p1)>=0) { *e=expo(x); avma=av; return gzero; }
! 1676: *e=expo(addsr(1,x)); avma=av; return negi(gun);
! 1677: }
! 1678: lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
! 1679: settyp(p1,t_INT); setlgefint(p1,lx);
! 1680: y=shifti(p1,e1); if (signe(x)<0) y=addsi(-1,y);
! 1681: y = gerepileupto(av,y);
! 1682:
! 1683: if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
! 1684: *e=e1; return y;
! 1685:
! 1686: case t_POLMOD: y=cgetg(3,t_POLMOD);
! 1687: copyifstack(x[1],y[1]);
! 1688: y[2]=lrndtoi((GEN)x[2],e); return y;
! 1689:
! 1690: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
! 1691: case t_VEC: case t_COL: case t_MAT:
! 1692: lx=(tx==t_POL)? lgef(x): lg(x);
! 1693: y=cgetg(lx,tx);
! 1694: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
! 1695: for ( ; i<lx; i++)
! 1696: {
! 1697: y[i]=lrndtoi((GEN)x[i],&e1);
! 1698: if (e1>*e) *e=e1;
! 1699: }
! 1700: return y;
! 1701: }
! 1702: err(typeer,"grndtoi");
! 1703: return NULL; /* not reached */
! 1704: }
! 1705:
! 1706: /* e = number of error bits on integral part */
! 1707: GEN
! 1708: gcvtoi(GEN x, long *e)
! 1709: {
! 1710: long tx=typ(x), lx,i,ex,av,e1;
! 1711: GEN y;
! 1712:
! 1713: *e = -HIGHEXPOBIT;
! 1714: if (tx == t_REAL)
! 1715: {
! 1716: long x0,x1;
! 1717: ex=expo(x); if (ex<0) { *e=ex; return gzero; }
! 1718: lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
! 1719: x0=x[0]; x1=x[1]; settyp(x,t_INT); setlgefint(x,lx);
! 1720: y=shifti(x,e1); x[0]=x0; x[1]=x1;
! 1721: if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
! 1722: *e=e1; return y;
! 1723: }
! 1724: if (is_matvec_t(tx))
! 1725: {
! 1726: lx=lg(x); y=cgetg(lx,tx);
! 1727: for (i=1; i<lx; i++)
! 1728: {
! 1729: y[i]=lcvtoi((GEN)x[i],&e1);
! 1730: if (e1>*e) *e=e1;
! 1731: }
! 1732: return y;
! 1733: }
! 1734: return gtrunc(x);
! 1735: }
! 1736:
! 1737: GEN
! 1738: gtrunc(GEN x)
! 1739: {
! 1740: long tx=typ(x),av,tetpil,i,v;
! 1741: GEN y;
! 1742:
! 1743: switch(tx)
! 1744: {
! 1745: case t_INT: case t_POL:
! 1746: return gcopy(x);
! 1747:
! 1748: case t_REAL:
! 1749: return mptrunc(x);
! 1750:
! 1751: case t_FRAC: case t_FRACN:
! 1752: return divii((GEN)x[1],(GEN)x[2]);
! 1753:
! 1754: case t_PADIC:
! 1755: if (!signe(x[4])) return gzero;
! 1756: v = valp(x);
! 1757: if (!v) return gcopy((GEN)x[4]);
! 1758: if (v>0)
! 1759: { /* here p^v is an integer */
! 1760: av=avma; y=gpuigs((GEN)x[2],v); tetpil=avma;
! 1761: return gerepile(av,tetpil, mulii(y,(GEN)x[4]));
! 1762: }
! 1763: y=cgetg(3,t_FRAC);
! 1764: y[1]=licopy((GEN)x[4]);
! 1765: y[2]=lpuigs((GEN)x[2],-v);
! 1766: return y;
! 1767:
! 1768: case t_RFRAC: case t_RFRACN:
! 1769: return gdeuc((GEN)x[1],(GEN)x[2]);
! 1770:
! 1771: case t_SER:
! 1772: return gconvsp(x,1);
! 1773:
! 1774: case t_VEC: case t_COL: case t_MAT:
! 1775: {
! 1776: long lx=lg(x);
! 1777:
! 1778: y=cgetg(lx,tx);
! 1779: for (i=1; i<lx; i++) y[i]=ltrunc((GEN)x[i]);
! 1780: return y;
! 1781: }
! 1782: }
! 1783: err(typeer,"gtrunc");
! 1784: return NULL; /* not reached */
! 1785: }
! 1786:
! 1787: GEN
! 1788: trunc0(GEN x, GEN *pte)
! 1789: {
! 1790: if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
! 1791: return gtrunc(x);
! 1792: }
! 1793:
! 1794: /*******************************************************************/
! 1795: /* */
! 1796: /* CONVERSIONS GEN --> POLYNOMIALS & SERIES */
! 1797: /* */
! 1798: /*******************************************************************/
! 1799:
! 1800: GEN
! 1801: zeropol(long v)
! 1802: {
! 1803: GEN x = cgetg(2,t_POL);
! 1804: x[1] = evallgef(2) | evalvarn(v); return x;
! 1805: }
! 1806:
! 1807: GEN
! 1808: scalarpol(GEN x, long v)
! 1809: {
! 1810: GEN y=cgetg(3,t_POL);
! 1811: y[1] = gcmp0(x)? evallgef(3) | evalvarn(v)
! 1812: : evallgef(3) | evalvarn(v) | evalsigne(1);
! 1813: y[2]=lcopy(x); return y;
! 1814: }
! 1815:
! 1816: static GEN
! 1817: gtopoly0(GEN x, long v, int reverse)
! 1818: {
! 1819: long tx=typ(x),lx,i,j;
! 1820: GEN y;
! 1821:
! 1822: if (v<0) v = 0;
! 1823: if (isexactzero(x)) return zeropol(v);
! 1824: if (is_scalar_t(tx)) return scalarpol(x,v);
! 1825:
! 1826: switch(tx)
! 1827: {
! 1828: case t_POL:
! 1829: y=gcopy(x); break;
! 1830: case t_SER:
! 1831: y=gconvsp(x,1); break;
! 1832: case t_RFRAC: case t_RFRACN:
! 1833: y=gdeuc((GEN)x[1],(GEN)x[2]); break;
! 1834: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
! 1835: lx=lg(x);
! 1836: if (reverse)
! 1837: {
! 1838: while (lx-- && isexactzero((GEN)x[lx]));
! 1839: i=lx+2; y=cgetg(i,t_POL);
! 1840: y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
! 1841: for (j=2; j<i; j++) y[j]=lcopy((GEN)x[j-1]);
! 1842: }
! 1843: else
! 1844: {
! 1845: i=1; j=lx; while (lx-- && isexactzero((GEN)x[i++]));
! 1846: i=lx+2; y=cgetg(i,t_POL);
! 1847: y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
! 1848: lx = j-1;
! 1849: for (j=2; j<i; j++) y[j]=lcopy((GEN)x[lx--]);
! 1850: }
! 1851: break;
! 1852: default: err(typeer,"gtopoly");
! 1853: }
! 1854: setvarn(y,v); return y;
! 1855: }
! 1856:
! 1857: GEN
! 1858: gtopolyrev(GEN x, long v) { return gtopoly0(x,v,1); }
! 1859:
! 1860: GEN
! 1861: gtopoly(GEN x, long v) { return gtopoly0(x,v,0); }
! 1862:
! 1863: GEN
! 1864: zeroser(long v, long val)
! 1865: {
! 1866: GEN x = cgetg(2,t_SER);
! 1867: x[1]=evalvalp(val) | evalvarn(v); return x;
! 1868: }
! 1869:
! 1870: GEN
! 1871: scalarser(GEN x, long v, long prec)
! 1872: {
! 1873: GEN y=cgetg(prec+2,t_SER);
! 1874: long i;
! 1875:
! 1876: y[1]=evalsigne(1) | evalvalp(0) | evalvarn(v);
! 1877: y[2]=lcopy(x); for (i=3; i<=prec+1; i++) y[i]=zero;
! 1878: return y;
! 1879: }
! 1880:
! 1881: GEN
! 1882: gtoser(GEN x, long v)
! 1883: {
! 1884: long tx=typ(x),lx,i,j,l,av,tetpil;
! 1885: GEN y,p1,p2;
! 1886:
! 1887: if (v<0) v = 0;
! 1888: if (tx==t_SER) { y=gcopy(x); setvarn(y,v); return y; }
! 1889: if (isexactzero(x)) return zeroser(v,precdl);
! 1890: if (is_scalar_t(tx)) return scalarser(x,v,precdl);
! 1891:
! 1892: switch(tx)
! 1893: {
! 1894: case t_POL:
! 1895: lx=lgef(x); i=2; while (i<lx && gcmp0((GEN)x[i])) i++;
! 1896: l=lx-i; if (precdl>l) l=precdl;
! 1897: y=cgetg(l+2,t_SER);
! 1898: y[1] = evalsigne(1) | evalvalp(i-2) | evalvarn(v);
! 1899: for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
! 1900: for ( ; j<=l+1; j++) y[j]=zero;
! 1901: break;
! 1902:
! 1903: case t_RFRAC: case t_RFRACN:
! 1904: av=avma; p1=gtoser((GEN)x[1],v); p2=gtoser((GEN)x[2],v);
! 1905: tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
! 1906:
! 1907: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
! 1908: lx=lg(x); i=1; while (i<lx && isexactzero((GEN)x[i])) i++;
! 1909: y = cgetg(lx-i+2,t_SER);
! 1910: y[1] = evalsigne(1) | evalvalp(i-1) | evalvarn(v);
! 1911: for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
! 1912: break;
! 1913:
! 1914: default: err(typeer,"gtoser");
! 1915: }
! 1916: return y;
! 1917: }
! 1918:
! 1919: GEN
! 1920: gtovec(GEN x)
! 1921: {
! 1922: long tx,lx,i;
! 1923: GEN y;
! 1924:
! 1925: if (!x) return cgetg(1,t_VEC);
! 1926: tx = typ(x);
! 1927: if (is_scalar_t(tx) || is_rfrac_t(tx))
! 1928: {
! 1929: y=cgetg(2,t_VEC); y[1]=lcopy(x);
! 1930: return y;
! 1931: }
! 1932: if (is_graphicvec_t(tx))
! 1933: {
! 1934: lx=lg(x); y=cgetg(lx,t_VEC);
! 1935: for (i=1; i<lx; i++) y[i]=lcopy((GEN)x[i]);
! 1936: return y;
! 1937: }
! 1938: if (tx==t_POL)
! 1939: {
! 1940: lx=lgef(x); y=cgetg(lx-1,t_VEC);
! 1941: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[lx-i]);
! 1942: return y;
! 1943: }
! 1944: if (tx==t_LIST)
! 1945: {
! 1946: lx=lgef(x); y=cgetg(lx-1,t_VEC); x++;
! 1947: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
! 1948: return y;
! 1949: }
! 1950: if (!signe(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
! 1951: lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
! 1952: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
! 1953: return y;
! 1954: }
! 1955:
! 1956: GEN
! 1957: compo(GEN x, long n)
! 1958: {
! 1959: long l,tx=typ(x);
! 1960:
! 1961: if (tx==t_POL && n+1 >= lgef(x)) return gzero;
! 1962: if (tx==t_SER && !signe(x)) return gzero;
! 1963: if (!is_recursive_t(tx))
! 1964: err(talker, "this object doesn't have components (is a leaf)");
! 1965: l=lontyp[tx]+n-1;
! 1966: if (n<1 || l>=lg(x))
! 1967: err(talker,"nonexistent component");
! 1968: return gcopy((GEN)x[l]);
! 1969: }
! 1970:
! 1971: /* with respect to the main variable if v<0, with respect to the variable v
! 1972: otherwise. v ignored if x is not a polynomial/series. */
! 1973:
! 1974: GEN
! 1975: polcoeff0(GEN x, long n, long v)
! 1976: {
! 1977: long tx=typ(x),lx,ex,w,av,tetpil;
! 1978: GEN xinit;
! 1979:
! 1980: if (is_scalar_t(tx)) return n? gzero: gcopy(x);
! 1981:
! 1982: switch(tx)
! 1983: {
! 1984: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
! 1985: if (n<1 || n>=lg(x)) break;
! 1986: return gcopy((GEN)x[n]);
! 1987:
! 1988: case t_POL:
! 1989: if (n<0) return gzero;
! 1990: w=varn(x);
! 1991: if (v < 0 || v == w)
! 1992: return (n>=lgef(x)-2)? gzero: gcopy((GEN)x[n+2]);
! 1993: if (v < w) return n? gzero: gcopy(x);
! 1994: av=avma; xinit=x;
! 1995: x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
! 1996: if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
! 1997: if (typ(x) == t_POL)
! 1998: {
! 1999: if (n>=lgef(x)-2) { avma=av; return gzero; }
! 2000: x = (GEN)x[n+2];
! 2001: }
! 2002: else
! 2003: x = polcoeff0(x, n, 0);
! 2004: tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
! 2005:
! 2006: case t_SER:
! 2007: w=varn(x);
! 2008: if (v < 0 || v == w)
! 2009: {
! 2010: if (!signe(x)) return gzero;
! 2011: lx=lg(x); ex=valp(x); if (n<ex) return gzero;
! 2012: if (n>=ex+lx-2) break;
! 2013: return gcopy((GEN)x[n-ex+2]);
! 2014: }
! 2015: if (v < w) return n? gzero: gcopy(x);
! 2016: av=avma; xinit=x;
! 2017: x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
! 2018: if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
! 2019: if (gcmp0(x)) { avma=av; return gzero; }
! 2020: if (typ(x) == t_SER)
! 2021: {
! 2022: lx=lg(x); ex=valp(x); if (n<ex) return gzero;
! 2023: if (n>=ex+lx-2) break;
! 2024: x = (GEN)x[n-ex+2];
! 2025: }
! 2026: else
! 2027: x = polcoeff0(x, n, 0);
! 2028: tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
! 2029:
! 2030: case t_RFRAC: case t_RFRACN:
! 2031: w = precdl; av = avma;
! 2032: if (v<0) v = gvar(x);
! 2033: ex = ggval((GEN)x[2], polx[v]);
! 2034: precdl = n + ex + 1; x = gtoser(x, v); precdl = w;
! 2035: return gerepileupto(av, polcoeff0(x, n, v));
! 2036: }
! 2037: err(talker,"nonexistent component in truecoeff");
! 2038: return NULL; /* not reached */
! 2039: }
! 2040:
! 2041: GEN
! 2042: truecoeff(GEN x, long n)
! 2043: {
! 2044: return polcoeff0(x,n,-1);
! 2045: }
! 2046:
! 2047: GEN
! 2048: denom(GEN x)
! 2049: {
! 2050: long lx,i,av,tetpil;
! 2051: GEN s,t;
! 2052:
! 2053: switch(typ(x))
! 2054: {
! 2055: case t_INT: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
! 2056: return gun;
! 2057:
! 2058: case t_FRAC: case t_FRACN:
! 2059: return absi((GEN)x[2]);
! 2060:
! 2061: case t_COMPLEX:
! 2062: av=avma; t=denom((GEN)x[1]); s=denom((GEN)x[2]); tetpil=avma;
! 2063: return gerepile(av,tetpil,glcm(s,t));
! 2064:
! 2065: case t_QUAD:
! 2066: av=avma; t=denom((GEN)x[2]); s=denom((GEN)x[3]); tetpil=avma;
! 2067: return gerepile(av,tetpil,glcm(s,t));
! 2068:
! 2069: case t_POLMOD:
! 2070: return denom((GEN)x[2]);
! 2071:
! 2072: case t_RFRAC: case t_RFRACN:
! 2073: return gcopy((GEN)x[2]);
! 2074:
! 2075: case t_POL:
! 2076: return polun[varn(x)];
! 2077:
! 2078: case t_VEC: case t_COL: case t_MAT:
! 2079: lx=lg(x); if (lx==1) return gun;
! 2080: av = tetpil = avma; s = denom((GEN)x[1]);
! 2081: for (i=2; i<lx; i++)
! 2082: {
! 2083: t = denom((GEN)x[i]);
! 2084: /* t != gun est volontaire */
! 2085: if (t != gun) { tetpil=avma; s=glcm(s,t); }
! 2086: }
! 2087: return gerepile(av,tetpil,s);
! 2088: }
! 2089: err(typeer,"denom");
! 2090: return NULL; /* not reached */
! 2091: }
! 2092:
! 2093: GEN
! 2094: numer(GEN x)
! 2095: {
! 2096: long av,tetpil;
! 2097: GEN s;
! 2098:
! 2099: switch(typ(x))
! 2100: {
! 2101: case t_INT: case t_REAL: case t_INTMOD:
! 2102: case t_PADIC: case t_POL: case t_SER:
! 2103: return gcopy(x);
! 2104:
! 2105: case t_FRAC: case t_FRACN:
! 2106: return (signe(x[2])>0)? gcopy((GEN)x[1]): gneg((GEN)x[1]);
! 2107:
! 2108: case t_POLMOD:
! 2109: av=avma; s=numer((GEN)x[2]); tetpil=avma;
! 2110: return gerepile(av,tetpil,gmodulcp(s,(GEN)x[1]));
! 2111:
! 2112: case t_RFRAC: case t_RFRACN:
! 2113: return gcopy((GEN)x[1]);
! 2114:
! 2115: case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
! 2116: av=avma; s=denom(x); tetpil=avma;
! 2117: return gerepile(av,tetpil,gmul(s,x));
! 2118: }
! 2119: err(typeer,"numer");
! 2120: return NULL; /* not reached */
! 2121: }
! 2122:
! 2123: /* Lift only intmods if v does not occur in x, lift with respect to main
! 2124: * variable of x if v < 0, with respect to variable v otherwise.
! 2125: */
! 2126: GEN
! 2127: lift0(GEN x, long v)
! 2128: {
! 2129: long lx,tx=typ(x),i;
! 2130: GEN y;
! 2131:
! 2132: switch(tx)
! 2133: {
! 2134: case t_INT: case t_REAL:
! 2135: return gcopy(x);
! 2136:
! 2137: case t_INTMOD:
! 2138: return gcopy((GEN)x[2]);
! 2139:
! 2140: case t_POLMOD:
! 2141: if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
! 2142: y = cgetg(3,tx);
! 2143: y[1] = (long)lift0((GEN)x[1],v);
! 2144: y[2] = (long)lift0((GEN)x[2],v); return y;
! 2145:
! 2146: case t_SER:
! 2147: if (!signe(x)) return gcopy(x);
! 2148: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 2149: for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
! 2150: return y;
! 2151:
! 2152: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
! 2153: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 2154: lx=lg(x); y=cgetg(lx,tx);
! 2155: for (i=1; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
! 2156: return y;
! 2157:
! 2158: case t_POL:
! 2159: y=cgetg(lx=lgef(x),tx); y[1]=x[1];
! 2160: for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
! 2161: return y;
! 2162:
! 2163: case t_QUAD:
! 2164: y=cgetg(4,tx); copyifstack(x[1],y[1]);
! 2165: for (i=2; i<4; i++) y[i]=(long)lift0((GEN)x[i],v);
! 2166: return y;
! 2167: }
! 2168: err(typeer,"lift");
! 2169: return NULL; /* not reached */
! 2170: }
! 2171:
! 2172: GEN
! 2173: lift(GEN x)
! 2174: {
! 2175: return lift0(x,-1);
! 2176: }
! 2177:
! 2178: /* same as lift, without copy. May DESTROY x. For internal use only.
! 2179: Conventions on v as for lift. */
! 2180: GEN
! 2181: lift_intern0(GEN x, long v)
! 2182: {
! 2183: long i,lx,tx=typ(x);
! 2184:
! 2185: switch(tx)
! 2186: {
! 2187: case t_INT: case t_REAL:
! 2188: return x;
! 2189:
! 2190: case t_INTMOD:
! 2191: return (GEN)x[2];
! 2192:
! 2193: case t_POLMOD:
! 2194: if (v < 0 || v == varn((GEN)x[1])) return (GEN)x[2];
! 2195: x[1]=(long)lift_intern0((GEN)x[1],v);
! 2196: x[2]=(long)lift_intern0((GEN)x[2],v);
! 2197: return x;
! 2198:
! 2199: case t_SER: if (!signe(x)) return x; /* fall through */
! 2200: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_QUAD: case t_POL:
! 2201: case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 2202: lx = (tx==t_POL)? lgef(x): lg(x);
! 2203: for (i = lx-1; i>=lontyp[tx]; i--)
! 2204: x[i] = (long) lift_intern0((GEN)x[i],v);
! 2205: return x;
! 2206: }
! 2207: err(typeer,"lift_intern");
! 2208: return NULL; /* not reached */
! 2209: }
! 2210:
! 2211: /* memes conventions pour v que lift */
! 2212: GEN
! 2213: centerlift0(GEN x, long v)
! 2214: {
! 2215: long lx,tx=typ(x),i,av;
! 2216: GEN y;
! 2217:
! 2218: switch(tx)
! 2219: {
! 2220: case t_INT:
! 2221: return icopy(x);
! 2222:
! 2223: case t_INTMOD:
! 2224: av=avma; i=cmpii(shifti((GEN)x[2],1),(GEN)x[1]); avma=av;
! 2225: return (i>0)? subii((GEN)x[2],(GEN)x[1]): icopy((GEN)x[2]);
! 2226:
! 2227: case t_POLMOD:
! 2228: if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
! 2229: y=cgetg(3,tx);
! 2230: y[1]=(long)centerlift0((GEN)x[1],v); y[2]=(long)centerlift0((GEN)x[2],v);
! 2231: return y;
! 2232:
! 2233: case t_SER:
! 2234: if (!signe(x)) return gcopy(x);
! 2235: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
! 2236: for (i=2; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
! 2237: return y;
! 2238:
! 2239: case t_POL:
! 2240: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
! 2241: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
! 2242: lx = (tx==t_POL)? lgef(x): lg(x);
! 2243: y=cgetg(lx,tx); y[1]=x[1];
! 2244: for (i=lontyp[tx]; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
! 2245: return y;
! 2246:
! 2247: case t_QUAD:
! 2248: y=cgetg(4,tx); copyifstack(x[1],y[1]);
! 2249: for (i=2; i<4; i++) y[i]=(long)centerlift0((GEN)x[i],v);
! 2250: return y;
! 2251: }
! 2252: err(typeer,"centerlift");
! 2253: return NULL; /* not reached */
! 2254: }
! 2255:
! 2256: GEN
! 2257: centerlift(GEN x)
! 2258: {
! 2259: return centerlift0(x,-1);
! 2260: }
! 2261:
! 2262: /*******************************************************************/
! 2263: /* */
! 2264: /* PARTIES REELLE ET IMAGINAIRES */
! 2265: /* */
! 2266: /*******************************************************************/
! 2267:
! 2268: static GEN
! 2269: op_ReIm(GEN f(GEN), GEN x)
! 2270: {
! 2271: long lx,i,j,av,tetpil, tx = typ(x);
! 2272: GEN p1,p2,r2,i2,z;
! 2273:
! 2274: switch(tx)
! 2275: {
! 2276: case t_POL:
! 2277: lx=lgef(x); av=avma;
! 2278: for (i=lx-1; i>=2; i--)
! 2279: if (!gcmp0(f((GEN)x[i]))) break;
! 2280: avma=av; if (i==1) return zeropol(varn(x));
! 2281:
! 2282: z=cgetg(i+1,tx); z[1]=evalsigne(1)|evallgef(1+i)|evalvarn(varn(x));
! 2283: for (j=2; j<=i; j++) z[j] = (long)f((GEN)x[j]);
! 2284: return z;
! 2285:
! 2286: case t_SER:
! 2287: if (gcmp0(x)) { z=cgetg(2,tx); z[1]=x[1]; return z; }
! 2288: lx=lg(x); av=avma;
! 2289: for (i=2; i<lx; i++)
! 2290: if (!gcmp0(f((GEN)x[i]))) break;
! 2291: avma=av; if (i==lx) return zeroser(varn(x),lx-2+valp(x));
! 2292:
! 2293: z=cgetg(lx-i+2,tx); z[1]=x[1]; setvalp(z, valp(x)+i-2);
! 2294: for (j=2; i<lx; j++,i++) z[j] = (long) f((GEN)x[i]);
! 2295: return z;
! 2296:
! 2297: case t_RFRAC: case t_RFRACN:
! 2298: {
! 2299: av=avma; r2=greal((GEN)x[2]); i2=gimag((GEN)x[2]);
! 2300: if (f==greal)
! 2301: p1 = gadd(gmul(greal((GEN)x[1]),r2),
! 2302: gmul(gimag((GEN)x[1]),i2));
! 2303: else
! 2304: p1 = gsub(gmul(gimag((GEN)x[1]),r2),
! 2305: gmul(greal((GEN)x[1]),i2));
! 2306: p2=gadd(gsqr(r2), gsqr(i2));
! 2307: tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
! 2308: }
! 2309:
! 2310: case t_VEC: case t_COL: case t_MAT:
! 2311: lx=lg(x); z=cgetg(lx,tx);
! 2312: for (i=1; i<lx; i++) z[i] = (long) f((GEN)x[i]);
! 2313: return z;
! 2314: }
! 2315: err(typeer,"greal/gimag");
! 2316: return NULL; /* not reached */
! 2317: }
! 2318:
! 2319: GEN
! 2320: greal(GEN x)
! 2321: {
! 2322: switch(typ(x))
! 2323: {
! 2324: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 2325: return gcopy(x);
! 2326:
! 2327: case t_COMPLEX:
! 2328: return gcopy((GEN)x[1]);
! 2329:
! 2330: case t_QUAD:
! 2331: return gcopy((GEN)x[2]);
! 2332: }
! 2333: return op_ReIm(greal,x);
! 2334: }
! 2335:
! 2336: GEN
! 2337: gimag(GEN x)
! 2338: {
! 2339: switch(typ(x))
! 2340: {
! 2341: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 2342: return gzero;
! 2343:
! 2344: case t_COMPLEX:
! 2345: return gcopy((GEN)x[2]);
! 2346:
! 2347: case t_QUAD:
! 2348: return gcopy((GEN)x[3]);
! 2349: }
! 2350: return op_ReIm(gimag,x);
! 2351: }
! 2352:
! 2353: /*******************************************************************/
! 2354: /* */
! 2355: /* LOGICAL OPERATIONS */
! 2356: /* */
! 2357: /*******************************************************************/
! 2358:
! 2359: GEN
! 2360: glt(GEN x, GEN y) { return gcmp(x,y)<0? gun: gzero; }
! 2361:
! 2362: GEN
! 2363: gle(GEN x, GEN y) { return gcmp(x,y)<=0? gun: gzero; }
! 2364:
! 2365: GEN
! 2366: gge(GEN x, GEN y) { return gcmp(x,y)>=0? gun: gzero; }
! 2367:
! 2368: GEN
! 2369: ggt(GEN x, GEN y) { return gcmp(x,y)>0? gun: gzero; }
! 2370:
! 2371: GEN
! 2372: geq(GEN x, GEN y) { return gegal(x,y)? gun: gzero; }
! 2373:
! 2374: GEN
! 2375: gne(GEN x, GEN y) { return gegal(x,y)? gzero: gun; }
! 2376:
! 2377: GEN
! 2378: gand(GEN x, GEN y) { return gcmp0(x)? gzero: (gcmp0(y)? gzero: gun); }
! 2379:
! 2380: GEN
! 2381: gor(GEN x, GEN y) { return gcmp0(x)? (gcmp0(y)? gzero: gun): gun; }
! 2382:
! 2383: GEN
! 2384: gnot(GEN x) { return gcmp0(x)? gun: gzero; }
! 2385:
! 2386: /*******************************************************************/
! 2387: /* */
! 2388: /* FORMAL SIMPLIFICATIONS */
! 2389: /* */
! 2390: /*******************************************************************/
! 2391:
! 2392: GEN
! 2393: geval(GEN x)
! 2394: {
! 2395: long av,tetpil,lx,i, tx = typ(x);
! 2396: GEN y,z;
! 2397:
! 2398: if (is_const_t(tx)) return gcopy(x);
! 2399: if (is_graphicvec_t(tx) || tx == t_RFRACN)
! 2400: {
! 2401: lx=lg(x); y=cgetg(lx, tx);
! 2402: for (i=1; i<lx; i++) y[i] = (long)geval((GEN)x[i]);
! 2403: return y;
! 2404: }
! 2405:
! 2406: switch(tx)
! 2407: {
! 2408: case t_STR:
! 2409: return flisexpr(GSTR(x));
! 2410:
! 2411: case t_POLMOD: y=cgetg(3,tx);
! 2412: y[1]=(long)geval((GEN)x[1]);
! 2413: av=avma; z=geval((GEN)x[2]); tetpil=avma;
! 2414: y[2]=lpile(av,tetpil,gmod(z,(GEN)y[1]));
! 2415: return y;
! 2416:
! 2417: case t_POL:
! 2418: lx=lgef(x); if (lx==2) return gzero;
! 2419: {
! 2420: entree *ep = varentries[varn(x)];
! 2421: if (!ep) return gcopy(x);
! 2422: z = (GEN)ep->value;
! 2423: }
! 2424: y=gzero; av=avma;
! 2425: for (i=lx-1; i>1; i--)
! 2426: {
! 2427: tetpil=avma;
! 2428: y = gadd(geval((GEN)x[i]), gmul(z,y));
! 2429: }
! 2430: return gerepile(av,tetpil,y);
! 2431:
! 2432: case t_SER:
! 2433: err(impl, "evaluation of a power series");
! 2434:
! 2435: case t_RFRAC:
! 2436: return gdiv(geval((GEN)x[1]),geval((GEN)x[2]));
! 2437: }
! 2438: err(typeer,"geval");
! 2439: return NULL; /* not reached */
! 2440: }
! 2441:
! 2442: GEN
! 2443: simplify(GEN x)
! 2444: {
! 2445: long tx=typ(x),av,tetpil,i,lx;
! 2446: GEN p1,y;
! 2447:
! 2448: switch(tx)
! 2449: {
! 2450: case t_INT: case t_REAL: case t_FRAC:
! 2451: return is_universal_constant(x)? x: gcopy(x);
! 2452:
! 2453: case t_INTMOD: case t_PADIC: case t_QFR: case t_QFI:
! 2454: case t_LIST: case t_STR: case t_VECSMALL:
! 2455: return gcopy(x);
! 2456:
! 2457: case t_FRACN:
! 2458: return gred(x);
! 2459:
! 2460: case t_COMPLEX:
! 2461: if (is_universal_constant(x)) return x;
! 2462: if (isexactzero((GEN)x[2])) return simplify((GEN)x[1]);
! 2463: y=cgetg(3,tx);
! 2464: y[1]=(long)simplify((GEN)x[1]);
! 2465: y[2]=(long)simplify((GEN)x[2]); return y;
! 2466:
! 2467: case t_QUAD:
! 2468: if (isexactzero((GEN)x[3])) return simplify((GEN)x[2]);
! 2469: y=cgetg(4,tx);
! 2470: y[1]=lcopy((GEN)x[1]);
! 2471: y[2]=(long)simplify((GEN)x[2]);
! 2472: y[3]=(long)simplify((GEN)x[3]); return y;
! 2473:
! 2474: case t_POLMOD: y=cgetg(3,tx);
! 2475: y[1]=(long)simplify((GEN)x[1]);
! 2476: av=avma; p1=gmod((GEN)x[2],(GEN)y[1]); tetpil=avma;
! 2477: y[2]=lpile(av,tetpil,simplify(p1)); return y;
! 2478:
! 2479: case t_POL:
! 2480: lx=lgef(x); if (lx==2) return gzero;
! 2481: if (lx==3) return simplify((GEN)x[2]);
! 2482: y=cgetg(lx,tx); y[1]=x[1];
! 2483: for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
! 2484: return y;
! 2485:
! 2486: case t_SER:
! 2487: if (!signe(x)) return gcopy(x);
! 2488: lx=lg(x);
! 2489: y=cgetg(lx,tx); y[1]=x[1];
! 2490: for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
! 2491: return y;
! 2492:
! 2493: case t_RFRAC: y=cgetg(3,tx);
! 2494: y[1]=(long)simplify((GEN)x[1]);
! 2495: y[2]=(long)simplify((GEN)x[2]); return y;
! 2496:
! 2497: case t_RFRACN:
! 2498: av=avma; y=gred_rfrac(x); tetpil=avma;
! 2499: return gerepile(av,tetpil,simplify(y));
! 2500:
! 2501: case t_VEC: case t_COL: case t_MAT:
! 2502: lx=lg(x); y=cgetg(lx,tx);
! 2503: for (i=1; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
! 2504: return y;
! 2505: }
! 2506: err(typeer,"simplify");
! 2507: return NULL; /* not reached */
! 2508: }
! 2509:
! 2510: /*******************************************************************/
! 2511: /* */
! 2512: /* EVALUATION OF SOME SIMPLE FORMAL OBJECTS */
! 2513: /* */
! 2514: /*******************************************************************/
! 2515: static GEN
! 2516: qfeval0_i(GEN q, GEN x, long n)
! 2517: {
! 2518: long i,j,av=avma;
! 2519: GEN res=gzero;
! 2520:
! 2521: for (i=2;i<n;i++)
! 2522: for (j=1;j<i;j++)
! 2523: res = gadd(res, gmul(gcoeff(q,i,j), mulii((GEN)x[i],(GEN)x[j])) );
! 2524: res=gshift(res,1);
! 2525: for (i=1;i<n;i++)
! 2526: res = gadd(res, gmul(gcoeff(q,i,i), sqri((GEN)x[i])) );
! 2527: return gerepileupto(av,res);
! 2528: }
! 2529:
! 2530: #if 0
! 2531: static GEN
! 2532: qfeval0(GEN q, GEN x, long n)
! 2533: {
! 2534: long i,j,av=avma;
! 2535: GEN res=gzero;
! 2536:
! 2537: for (i=2;i<n;i++)
! 2538: for (j=1;j<i;j++)
! 2539: res = gadd(res, gmul(gcoeff(q,i,j), gmul((GEN)x[i],(GEN)x[j])) );
! 2540: res=gshift(res,1);
! 2541: for (i=1;i<n;i++)
! 2542: res = gadd(res, gmul(gcoeff(q,i,i), gsqr((GEN)x[i])) );
! 2543: return gerepileupto(av,res);
! 2544: }
! 2545: #else
! 2546: static GEN
! 2547: qfeval0(GEN q, GEN x, long n)
! 2548: {
! 2549: long i,j,av=avma;
! 2550: GEN res = gmul(gcoeff(q,1,1), gsqr((GEN)x[1]));
! 2551:
! 2552: for (i=2; i<n; i++)
! 2553: {
! 2554: GEN l = (GEN)q[i];
! 2555: GEN sx = gmul((GEN)l[1], (GEN)x[1]);
! 2556: for (j=2; j<i; j++)
! 2557: sx = gadd(sx, gmul((GEN)l[j],(GEN)x[j]));
! 2558: sx = gadd(gshift(sx,1), gmul((GEN)l[i],(GEN)x[i]));
! 2559:
! 2560: res = gadd(res, gmul((GEN)x[i], sx));
! 2561: }
! 2562: return gerepileupto(av,res);
! 2563: }
! 2564: #endif
! 2565:
! 2566: /* We assume q is a real symetric matrix */
! 2567: GEN
! 2568: qfeval(GEN q, GEN x)
! 2569: {
! 2570: long n=lg(q);
! 2571:
! 2572: if (n==1)
! 2573: {
! 2574: if (typ(q) != t_MAT || lg(x) != 1)
! 2575: err(talker,"invalid data in qfeval");
! 2576: return gzero;
! 2577: }
! 2578: if (typ(q) != t_MAT || lg(q[1]) != n)
! 2579: err(talker,"invalid quadratic form in qfeval");
! 2580: if (typ(x) != t_COL || lg(x) != n)
! 2581: err(talker,"invalid vector in qfeval");
! 2582:
! 2583: return qfeval0(q,x,n);
! 2584: }
! 2585:
! 2586: /* the Horner-type evaluation (mul x 2/3) would force us to use gmul and not
! 2587: * mulii (more than 4 x slower for small entries). Not worth it.
! 2588: */
! 2589: static GEN
! 2590: qfbeval0_i(GEN q, GEN x, GEN y, long n)
! 2591: {
! 2592: long i,j,av=avma;
! 2593: GEN res = gmul(gcoeff(q,1,1), mulii((GEN)x[1],(GEN)y[1]));
! 2594:
! 2595: for (i=2;i<n;i++)
! 2596: {
! 2597: for (j=1;j<i;j++)
! 2598: {
! 2599: GEN p1 = addii(mulii((GEN)x[i],(GEN)y[j]), mulii((GEN)x[j],(GEN)y[i]));
! 2600: res = gadd(res, gmul(gcoeff(q,i,j),p1));
! 2601: }
! 2602: res = gadd(res, gmul(gcoeff(q,i,i), mulii((GEN)x[i],(GEN)y[i])));
! 2603: }
! 2604: return gerepileupto(av,res);
! 2605: }
! 2606:
! 2607: #if 0
! 2608: static GEN
! 2609: qfbeval0(GEN q, GEN x, GEN y, long n)
! 2610: {
! 2611: long i,j,av=avma;
! 2612: GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1],(GEN)y[1]));
! 2613:
! 2614: for (i=2;i<n;i++)
! 2615: {
! 2616: for (j=1;j<i;j++)
! 2617: {
! 2618: GEN p1 = gadd(gmul((GEN)x[i],(GEN)y[j]), gmul((GEN)x[j],(GEN)y[i]));
! 2619: res = gadd(res, gmul(gcoeff(q,i,j),p1));
! 2620: }
! 2621: res = gadd(res, gmul(gcoeff(q,i,i), gmul((GEN)x[i],(GEN)y[i])));
! 2622: }
! 2623: return gerepileupto(av,res);
! 2624: }
! 2625: #else
! 2626: static GEN
! 2627: qfbeval0(GEN q, GEN x, GEN y, long n)
! 2628: {
! 2629: long i,j,av=avma;
! 2630: GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1], (GEN)y[1]));
! 2631:
! 2632: for (i=2; i<n; i++)
! 2633: {
! 2634: GEN l = (GEN)q[i];
! 2635: GEN sx = gmul((GEN)l[1], (GEN)y[1]);
! 2636: GEN sy = gmul((GEN)l[1], (GEN)x[1]);
! 2637: for (j=2; j<i; j++)
! 2638: {
! 2639: sx = gadd(sx, gmul((GEN)l[j],(GEN)y[j]));
! 2640: sy = gadd(sy, gmul((GEN)l[j],(GEN)x[j]));
! 2641: }
! 2642: sx = gadd(sx, gmul((GEN)l[i],(GEN)y[i]));
! 2643:
! 2644: sx = gmul((GEN)x[i], sx);
! 2645: sy = gmul((GEN)y[i], sy);
! 2646: res = gadd(res, gadd(sx,sy));
! 2647: }
! 2648: return gerepileupto(av,res);
! 2649: }
! 2650: #endif
! 2651:
! 2652: /* We assume q is a real symetric matrix */
! 2653: GEN
! 2654: qfbeval(GEN q, GEN x, GEN y)
! 2655: {
! 2656: long n=lg(q);
! 2657:
! 2658: if (n==1)
! 2659: {
! 2660: if (typ(q) != t_MAT || lg(x) != 1 || lg(y) != 1)
! 2661: err(talker,"invalid data in qfbeval");
! 2662: return gzero;
! 2663: }
! 2664: if (typ(q) != t_MAT || lg(q[1]) != n)
! 2665: err(talker,"invalid quadratic form in qfbeval");
! 2666: if (typ(x) != t_COL || lg(x) != n || typ(y) != t_COL || lg(y) != n)
! 2667: err(talker,"invalid vector in qfbeval");
! 2668:
! 2669: return qfbeval0(q,x,y,n);
! 2670: }
! 2671:
! 2672: /* yield X = M'.q.M, assuming q is symetric.
! 2673: * X_ij are X_ji identical, not copies
! 2674: * if flag is set, M has integer entries
! 2675: */
! 2676: GEN
! 2677: qf_base_change(GEN q, GEN M, int flag)
! 2678: {
! 2679: long i,j, k = lg(M), n=lg(q);
! 2680: GEN res = cgetg(k,t_MAT);
! 2681: GEN (*qf)(GEN,GEN,long) = flag? qfeval0_i: qfeval0;
! 2682: GEN (*qfb)(GEN,GEN,GEN,long) = flag? qfbeval0_i: qfbeval0;
! 2683:
! 2684: if (n==1)
! 2685: {
! 2686: if (typ(q) != t_MAT || k != 1)
! 2687: err(talker,"invalid data in qf_base_change");
! 2688: return res;
! 2689: }
! 2690: if (typ(M) != t_MAT || k == 1 || lg(M[1]) != n)
! 2691: err(talker,"invalid base change matrix in qf_base_change");
! 2692:
! 2693: for (i=1;i<k;i++)
! 2694: {
! 2695: res[i] = lgetg(k,t_COL);
! 2696: coeff(res,i,i) = (long) qf(q,(GEN)M[i],n);
! 2697: }
! 2698: for (i=2;i<k;i++)
! 2699: for (j=1;j<i;j++)
! 2700: coeff(res,i,j)=coeff(res,j,i) = (long) qfb(q,(GEN)M[i],(GEN)M[j],n);
! 2701: return res;
! 2702: }
! 2703:
! 2704: /* compute M'.M */
! 2705: GEN
! 2706: gram_matrix(GEN M)
! 2707: {
! 2708: long n=lg(M),i,j,k, av;
! 2709: GEN res = cgetg(n,t_MAT),p1;
! 2710:
! 2711: if (n==1)
! 2712: {
! 2713: if (typ(M) != t_MAT)
! 2714: err(talker,"invalid data in gram_matrix");
! 2715: return res;
! 2716: }
! 2717: if (typ(M) != t_MAT || lg(M[1]) != n)
! 2718: err(talker,"not a square matrix in gram_matrix");
! 2719:
! 2720: for (i=1;i<n;i++) res[i]=lgetg(n,t_COL);
! 2721: av=avma;
! 2722: for (i=1;i<n;i++)
! 2723: {
! 2724: p1 = gzero;
! 2725: for (k=1;k<n;k++)
! 2726: p1 = gadd(p1, gsqr(gcoeff(M,k,i)));
! 2727: coeff(res,i,i) = (long) gerepileupto(av,p1);
! 2728: av=avma;
! 2729: }
! 2730: for (i=2;i<n;i++)
! 2731: for (j=1;j<i;j++)
! 2732: {
! 2733: p1=gzero;
! 2734: for (k=1;k<n;k++)
! 2735: p1 = gadd(p1, gmul(gcoeff(M,k,i),gcoeff(M,k,j)));
! 2736: coeff(res,i,j)=coeff(res,j,i) = lpileupto(av,p1);
! 2737: av=avma;
! 2738: }
! 2739: return res;
! 2740: }
! 2741:
! 2742: /* return Re(x * y), x and y scalars */
! 2743: static GEN
! 2744: mul_real(GEN x, GEN y)
! 2745: {
! 2746: if (typ(x) == t_COMPLEX)
! 2747: {
! 2748: if (typ(y) == t_COMPLEX)
! 2749: {
! 2750: long av=avma, tetpil;
! 2751: GEN p1 = gmul((GEN)x[1], (GEN) y[1]);
! 2752: GEN p2 = gneg(gmul((GEN)x[2], (GEN) y[2]));
! 2753: tetpil=avma; return gerepile(av,tetpil,gadd(p1,p2));
! 2754: }
! 2755: x = (GEN)x[1];
! 2756: }
! 2757: else if (typ(y) == t_COMPLEX) y = (GEN)y[1];
! 2758: return gmul(x,y);
! 2759: }
! 2760:
! 2761: /* Compute Re(x * y), x and y matrices of compatible dimensions
! 2762: * assume lx, ly > 1, and scalar entries */
! 2763: GEN
! 2764: mulmat_real(GEN x, GEN y)
! 2765: {
! 2766: long av,tetpil,i,j,k, lx = lg(x), ly = lg(y), l = lg(x[1]);
! 2767: GEN p1,p2, z = cgetg(ly,t_MAT);
! 2768:
! 2769: for (j=1; j<ly; j++)
! 2770: {
! 2771: z[j] = lgetg(l,t_COL);
! 2772: for (i=1; i<l; i++)
! 2773: {
! 2774: p1=gzero; av=avma;
! 2775: for (k=1; k<lx; k++)
! 2776: {
! 2777: p2=mul_real(gcoeff(x,i,k),gcoeff(y,k,j));
! 2778: tetpil=avma; p1=gadd(p1,p2);
! 2779: }
! 2780: coeff(z,i,j)=lpile(av,tetpil,p1);
! 2781: }
! 2782: }
! 2783: return z;
! 2784: }
! 2785:
! 2786: static GEN
! 2787: hqfeval0(GEN q, GEN x, long n)
! 2788: {
! 2789: long i,j, av=avma;
! 2790: GEN res=gzero;
! 2791:
! 2792: for (i=2;i<n;i++)
! 2793: for (j=1;j<i;j++)
! 2794: {
! 2795: GEN p1 = gmul((GEN)x[i],gconj((GEN)x[j]));
! 2796: res = gadd(res, mul_real(gcoeff(q,i,j),p1));
! 2797: }
! 2798: res=gshift(res,1);
! 2799: for (i=1;i<n;i++)
! 2800: res = gadd(res, gmul(gcoeff(q,i,i), gnorm((GEN)x[i])) );
! 2801: return gerepileupto(av,res);
! 2802: }
! 2803:
! 2804: /* We assume q is a hermitian complex matrix */
! 2805: GEN
! 2806: hqfeval(GEN q, GEN x)
! 2807: {
! 2808: long n=lg(q);
! 2809:
! 2810: if (n==1)
! 2811: {
! 2812: if (typ(q) != t_MAT || lg(x) != 1)
! 2813: err(talker,"invalid data in hqfeval");
! 2814: return gzero;
! 2815: }
! 2816: if (typ(q) != t_MAT || lg(q[1]) != n)
! 2817: err(talker,"invalid quadratic form in hqfeval");
! 2818: if (typ(x) != t_COL || lg(x) != n)
! 2819: err(talker,"invalid vector in hqfeval");
! 2820:
! 2821: return hqfeval0(q,x,n);
! 2822: }
! 2823:
! 2824: GEN
! 2825: poleval(GEN x, GEN y)
! 2826: {
! 2827: long av,tetpil,i,j,imin,tx=typ(x);
! 2828: GEN p1,p2,p3,r,s;
! 2829:
! 2830: if (is_scalar_t(tx)) return gcopy(x);
! 2831: switch(tx)
! 2832: {
! 2833: case t_POL:
! 2834: i=lgef(x)-1; imin=2; break;
! 2835:
! 2836: case t_RFRAC: case t_RFRACN: av=avma;
! 2837: p1=poleval((GEN)x[1],y);
! 2838: p2=poleval((GEN)x[2],y); tetpil=avma;
! 2839: return gerepile(av,tetpil,gdiv(p1,p2));
! 2840:
! 2841: case t_VEC: case t_COL:
! 2842: i=lg(x)-1; imin=1; break;
! 2843: default: err(typeer,"poleval");
! 2844: }
! 2845: if (i<=imin)
! 2846: return (i==imin)? gcopy((GEN)x[imin]): gzero;
! 2847:
! 2848: av=avma; p1=(GEN)x[i]; i--;
! 2849: if (typ(y)!=t_COMPLEX)
! 2850: {
! 2851: #if 0 /* standard Horner's rule */
! 2852: for ( ; i>=imin; i--)
! 2853: p1 = gadd(gmul(p1,y),(GEN)x[i]);
! 2854: #endif
! 2855: /* specific attention to sparse polynomials */
! 2856: for ( ; i>=imin; i=j-1)
! 2857: {
! 2858: for (j=i; gcmp0((GEN)x[j]); j--)
! 2859: if (j==imin)
! 2860: {
! 2861: if (i!=j) y = gpuigs(y,i-j+1);
! 2862: tetpil=avma; return gerepile(av,tetpil,gmul(p1,y));
! 2863: }
! 2864: r = (i==j)? y: gpuigs(y,i-j+1);
! 2865: p1 = gadd(gmul(p1,r), (GEN)x[j]);
! 2866: }
! 2867: return gerepileupto(av,p1);
! 2868: }
! 2869:
! 2870: p2=(GEN)x[i]; i--; r=gtrace(y); s=gneg_i(gnorm(y));
! 2871: for ( ; i>=imin; i--)
! 2872: {
! 2873: p3=gadd(p2,gmul(r,p1));
! 2874: p2=gadd((GEN)x[i],gmul(s,p1)); p1=p3;
! 2875: }
! 2876: p1=gmul(y,p1); tetpil=avma;
! 2877: return gerepile(av,tetpil,gadd(p1,p2));
! 2878: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen3.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