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