Annotation of OpenXM_contrib/pari/src/basemath/gen1.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** GENERIC OPERATIONS **/
! 4: /** (first part) **/
! 5: /** **/
! 6: /********************************************************************/
! 7: /* $Id: gen1.c,v 1.2 1999/09/21 10:57:55 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
! 11: GEN quickmul(GEN a, GEN b, long na, long nb);
! 12:
! 13: #define cpifstack(x) isonstack(x)?gcopy(x):x
! 14: /* y is a polmod, f is gadd or gmul */
! 15: static GEN
! 16: op_polmod(GEN f(GEN,GEN), GEN x, GEN y, long tx)
! 17: {
! 18: GEN mod,k,l, z=cgetg(3,t_POLMOD);
! 19: long av,tetpil;
! 20:
! 21: l=(GEN)y[1];
! 22: if (tx==t_POLMOD)
! 23: {
! 24: k=(GEN)x[1];
! 25: if (gegal(k,l))
! 26: { mod=cpifstack(k); x=(GEN)x[2]; y=(GEN)y[2]; }
! 27: else
! 28: {
! 29: long vx=varn(k), vy=varn(l);
! 30: if (vx==vy) { mod=srgcd(k,l); x=(GEN)x[2]; y=(GEN)y[2]; }
! 31: else
! 32: if (vx<vy) { mod=cpifstack(k); x=(GEN)x[2]; }
! 33: else { mod=cpifstack(l); y=(GEN)y[2]; }
! 34: }
! 35: }
! 36: else
! 37: {
! 38: mod=cpifstack(l); y=(GEN)y[2];
! 39: if (is_scalar_t(tx))
! 40: {
! 41: z[2] = (long)f(x,y);
! 42: z[1] = (long)mod; return z;
! 43: }
! 44: }
! 45: av=avma; x = f(x,y); tetpil=avma;
! 46: z[2] = lpile(av,tetpil,gmod(x,mod));
! 47: z[1] = (long)mod; return z;
! 48: }
! 49: #undef cpifstack
! 50:
! 51: /********************************************************************/
! 52: /** **/
! 53: /** SUBTRACTION **/
! 54: /** **/
! 55: /********************************************************************/
! 56:
! 57: GEN
! 58: gsub(GEN x, GEN y)
! 59: {
! 60: long tetpil, av = avma;
! 61: y=gneg_i(y); tetpil=avma;
! 62: return gerepile(av,tetpil,gadd(x,y));
! 63: }
! 64:
! 65: /********************************************************************/
! 66: /** **/
! 67: /** ADDITION **/
! 68: /** **/
! 69: /********************************************************************/
! 70:
! 71: static GEN
! 72: addpadic(GEN x, GEN y)
! 73: {
! 74: long c,e,r,d,r1,r2,av,tetpil;
! 75: GEN z,p1,p2, p = (GEN)x[2];
! 76:
! 77: z=cgetg(5,t_PADIC); icopyifstack(p, z[2]); av=avma;
! 78: e=valp(x); r=valp(y); d = r-e;
! 79: if (d<0) { p1=x; x=y; y=p1; e=r; d = -d; }
! 80: r1=precp(x); r2=precp(y);
! 81: if (d)
! 82: {
! 83: r = d+r2;
! 84: p1 = (d==1)? p: gclone(gpuigs(p,d));
! 85: avma=av;
! 86: if (r<r1) z[3]=lmulii(p1,(GEN)y[3]);
! 87: else
! 88: {
! 89: r=r1; z[3]=licopy((GEN)x[3]);
! 90: }
! 91: av=avma; p2=mulii(p1,(GEN)y[4]);
! 92: if (d!=1) gunclone(p1);
! 93: p1=addii(p2,(GEN)x[4]); tetpil=avma;
! 94: z[4]=lpile(av,tetpil, modii(p1,(GEN)z[3]));
! 95: z[1]=evalprecp(r) | evalvalp(e); return z;
! 96: }
! 97: if (r2<r1) { r=r2; p1=x; x=y; y=p1; } else r=r1;
! 98: p1 = addii((GEN)x[4],(GEN)y[4]);
! 99: if (!signe(p1) || (c = pvaluation(p1,p,&p2)) >=r)
! 100: {
! 101: avma=av; z[4]=zero; z[3]=un;
! 102: z[1]=evalvalp(e+r); return z;
! 103: }
! 104: if (c)
! 105: {
! 106: p2=gclone(p2); avma=av;
! 107: if (c==1)
! 108: z[3] = ldivii((GEN)x[3], p);
! 109: else
! 110: {
! 111: p1 = gpuigs(p,c); tetpil=avma;
! 112: z[3] = lpile(av,tetpil, divii((GEN)x[3], p1));
! 113: }
! 114: z[4]=lmodii(p2,(GEN)z[3]); gunclone(p2);
! 115: z[1]=evalprecp(r-c) | evalvalp(e+c); return z;
! 116: }
! 117: tetpil=avma;
! 118: z[4]=lpile(av,tetpil,modii(p1,(GEN)x[3]));
! 119: z[3]=licopy((GEN)x[3]);
! 120: z[1]=evalprecp(r) | evalvalp(e); return z;
! 121: }
! 122:
! 123: /* return x + y, where x is t_INT or t_FRAC(N), y t_PADIC */
! 124: static GEN
! 125: gaddpex(GEN x, GEN y)
! 126: {
! 127: long tx,e1,e2,e3,av,tetpil;
! 128: GEN z,p,p1,p2;
! 129:
! 130: if (gcmp0(x)) return gcopy(y);
! 131:
! 132: av=avma; p=(GEN)y[2]; tx=typ(x);
! 133: z=cgetg(5,t_PADIC); z[2]=(long)p;
! 134: e3 = (tx == t_INT)? pvaluation(x,p,&p1)
! 135: : pvaluation((GEN)x[1],p,&p1) -
! 136: pvaluation((GEN)x[2],p,&p2);
! 137: e1 = valp(y)-e3; e2 = signe(y[4])? e1+precp(y): e1;
! 138: if (e2<=0)
! 139: {
! 140: z[1] = evalprecp(0) | evalvalp(e3);
! 141: z[3] = un;
! 142: z[4] = zero;
! 143: }
! 144: else
! 145: {
! 146: if (tx != t_INT && !is_pm1(p2)) p1 = gdiv(p1,p2);
! 147: z[1] = evalprecp(e2) | evalvalp(e3);
! 148: z[3] = e1? lmul((GEN)y[3], gpuigs(p,e1)): y[3];
! 149: z[4] = lmod(p1,(GEN)z[3]);
! 150: }
! 151: tetpil=avma; return gerepile(av,tetpil,addpadic(z,y));
! 152: }
! 153:
! 154: static long
! 155: kro_quad(GEN x, GEN y)
! 156: {
! 157: long k, av=avma;
! 158:
! 159: x = subii(sqri((GEN)x[3]), shifti((GEN)x[2],2));
! 160: k = kronecker(x,y); avma=av; return k;
! 161: }
! 162:
! 163: static GEN
! 164: addfrac(GEN x, GEN y)
! 165: {
! 166: GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
! 167: GEN y1 = (GEN)y[1], y2 = (GEN)y[2], p1,p2,n,d,delta,z;
! 168:
! 169: z = cgetg(3,t_FRAC);
! 170: (void)new_chunk((lgefint(x1) + lgefint(x2) +
! 171: lgefint(y1) + lgefint(y2)) << 1);
! 172: delta = mppgcd(x2,y2);
! 173: if (is_pm1(delta))
! 174: {
! 175: p1 = mulii(x1,y2);
! 176: p2 = mulii(y1,x2); avma = (long)z;
! 177: z[1] = laddii(p1,p2);
! 178: z[2] = lmulii(x2,y2); return z;
! 179: }
! 180: x2 = divii(x2,delta);
! 181: y2 = divii(y2,delta);
! 182: n = addii(mulii(x1,y2), mulii(y1,x2));
! 183: if (!signe(n)) { avma = (long)(z+3); return gzero; }
! 184: d = mulii(x2, y2);
! 185: p1 = dvmdii(n, delta, &p2);
! 186: if (p2 == gzero)
! 187: {
! 188: if (is_pm1(d)) { avma = (long)(z+3); return icopy(p1); }
! 189: avma = (long)z;
! 190: z[1] = licopy(p1);
! 191: z[2] = licopy(d); return z;
! 192: }
! 193: p1 = mppgcd(delta, p2);
! 194: if (is_pm1(p1))
! 195: {
! 196: avma = (long)z;
! 197: z[1] = licopy(n);
! 198: }
! 199: else
! 200: {
! 201: delta = divii(delta, p1);
! 202: avma = (long)z;
! 203: z[1] = ldivii(n,p1);
! 204: }
! 205: z[2] = lmulii(d,delta); return z;
! 206: }
! 207:
! 208: static GEN
! 209: addrfrac(GEN x, GEN y)
! 210: {
! 211: GEN z = cgetg(3,t_RFRAC);
! 212: GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
! 213: GEN y1 = (GEN)y[1], y2 = (GEN)y[2], p1,p2,n,d,delta;
! 214: long tetpil;
! 215:
! 216: delta = ggcd(x2,y2);
! 217: if (gcmp1(delta))
! 218: {
! 219: p1 = gmul(x1,y2);
! 220: p2 = gmul(y1,x2);
! 221: tetpil = avma; /* numerator is non-zero */
! 222: z[1] = lpile((long)z,tetpil, gadd(p1,p2));
! 223: z[2] = lmul(x2, y2); return z;
! 224: }
! 225: x2 = gdeuc(x2,delta);
! 226: y2 = gdeuc(y2,delta);
! 227: n = gadd(gmul(x1,y2), gmul(y1,x2));
! 228: if (!signe(n)) return gerepileupto((long)(z+3), n);
! 229: tetpil = avma; d = gmul(x2, y2);
! 230: p1 = poldivres(n, delta, &p2);
! 231: if (!signe(p2))
! 232: {
! 233: if (gcmp1(d)) return gerepileupto((long)(z+3), p1);
! 234: z[1]=(long)p1; z[2]=(long)d;
! 235: gerepilemanyvec((long)z,tetpil,z+1,2); return z;
! 236: }
! 237: p1 = ggcd(delta, p2);
! 238: if (gcmp1(p1))
! 239: {
! 240: tetpil = avma;
! 241: z[1] = lcopy(n);
! 242: }
! 243: else
! 244: {
! 245: delta = gdeuc(delta, p1);
! 246: tetpil = avma;
! 247: z[1] = ldeuc(n,p1);
! 248: }
! 249: z[2] = lmul(d,delta);
! 250: gerepilemanyvec((long)z,tetpil,z+1,2); return z;
! 251: }
! 252:
! 253: static GEN
! 254: addscalrfrac(GEN x, GEN y)
! 255: {
! 256: GEN p1,num, z = cgetg(3,t_RFRAC);
! 257: long tetpil, av;
! 258:
! 259: p1 = gmul(x,(GEN)y[2]); tetpil = avma;
! 260: num = gadd(p1,(GEN)y[1]);
! 261: av = avma;
! 262: p1 = content((GEN)y[2]);
! 263: if (!gcmp1(p1))
! 264: {
! 265: p1 = ggcd(p1, content(num));
! 266: if (!gcmp1(p1))
! 267: {
! 268: tetpil = avma;
! 269: z[1] = ldiv(num, p1);
! 270: z[2] = ldiv((GEN)y[2], p1);
! 271: gerepilemanyvec((long)z,tetpil,z+1,2); return z;
! 272: }
! 273: }
! 274: avma = av;
! 275: z[1]=lpile((long)z,tetpil, num);
! 276: z[2]=lcopy((GEN)y[2]); return z;
! 277: }
! 278:
! 279: /* assume gvar(x) = varn(mod) */
! 280: static GEN
! 281: to_polmod(GEN x, GEN mod)
! 282: {
! 283: long tx = typ(x);
! 284: GEN z = cgetg(3, t_POLMOD);
! 285:
! 286: if (tx == t_RFRACN) { x = gred_rfrac(x); tx = t_RFRAC; }
! 287: if (tx == t_RFRAC) x = gmul((GEN)x[1], ginvmod((GEN)x[2],mod));
! 288: z[1] = (long)mod;
! 289: z[2] = (long)x;
! 290: return z;
! 291: }
! 292:
! 293: GEN
! 294: gadd(GEN x, GEN y)
! 295: {
! 296: long vx,vy,lx,ly,tx,ty,i,j,k,l,av,tetpil;
! 297: GEN z,p1,p2;
! 298:
! 299: tx=typ(x); ty=typ(y);
! 300: if (is_const_t(tx) && is_const_t(ty))
! 301: {
! 302: if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
! 303: }
! 304: else
! 305: {
! 306: vx=gvar(x); vy=gvar(y);
! 307: if (vx<vy || (vx==vy && tx>ty))
! 308: {
! 309: p1=x; x=y; y=p1;
! 310: i=tx; tx=ty; ty=i;
! 311: i=vx; vx=vy; vy=i;
! 312: }
! 313: if (ty==t_POLMOD) return op_polmod(gadd,x,y,tx);
! 314: }
! 315:
! 316: /* here vx >= vy */
! 317: if (is_scalar_t(ty))
! 318: {
! 319: switch(tx)
! 320: {
! 321: case t_INT:
! 322: switch(ty)
! 323: {
! 324: case t_INT: return addii(x,y);
! 325: case t_REAL: return addir(x,y);
! 326:
! 327: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
! 328: (void)new_chunk(lgefint(p2)+1);
! 329: p1 = addii(modii(x,p2),(GEN)y[2]); avma = (long)z;
! 330: z[2] = (cmpii(p1,p2) >=0)? lsubii(p1,p2): licopy(p1);
! 331: icopyifstack(p2,z[1]); return z;
! 332:
! 333: case t_FRAC: case t_FRACN: z=cgetg(3,ty);
! 334: (void)new_chunk(lgefint(x) + lgefint(y[1]) + lgefint(y[2]) + 1);
! 335: p1 = mulii((GEN)y[2],x); avma = (long)z;
! 336: z[1] = laddii((GEN)y[1], p1);
! 337: z[2] = licopy((GEN)y[2]); return z;
! 338:
! 339: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 340: z[1]=ladd(x,(GEN)y[1]);
! 341: z[2]=lcopy((GEN)y[2]); return z;
! 342:
! 343: case t_PADIC:
! 344: return gaddpex(x,y);
! 345:
! 346: case t_QUAD: z=cgetg(4,t_QUAD);
! 347: copyifstack(y[1], z[1]);
! 348: z[2]=ladd(x,(GEN)y[2]);
! 349: z[3]=lcopy((GEN)y[3]); return z;
! 350: }
! 351:
! 352: case t_REAL:
! 353: switch(ty)
! 354: {
! 355: case t_REAL: return addrr(x,y);
! 356:
! 357: case t_FRAC: case t_FRACN:
! 358: if (!signe(y[1])) return rcopy(x);
! 359: if (!signe(x))
! 360: {
! 361: lx = expi((GEN)y[1]) - expi((GEN)y[2]) - expo(x);
! 362: if (lx < 0) return rcopy(x);
! 363: lx >>= TWOPOTBITS_IN_LONG;
! 364: z=cgetr(lx+3); diviiz((GEN)y[1],(GEN)y[2],z);
! 365: return z;
! 366: }
! 367: av=avma; z=addir((GEN)y[1],mulir((GEN)y[2],x)); tetpil=avma;
! 368: return gerepile(av,tetpil,divri(z,(GEN)y[2]));
! 369:
! 370: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 371: z[1]=ladd(x,(GEN)y[1]);
! 372: z[2]=lcopy((GEN)y[2]); return z;
! 373:
! 374: case t_QUAD:
! 375: if (gcmp0(y)) return rcopy(x);
! 376:
! 377: av=avma; i=gexpo(y)-expo(x);
! 378: if (i<=0) i=0; else i >>= TWOPOTBITS_IN_LONG;
! 379: p1=co8(y,lg(x)+i); tetpil=avma;
! 380: return gerepile(av,tetpil,gadd(p1,x));
! 381:
! 382: case t_INTMOD: case t_PADIC: err(gadderf,tx,ty);
! 383: }
! 384:
! 385: case t_INTMOD:
! 386: switch(ty)
! 387: {
! 388: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
! 389: if (p1==p2 || egalii(p1,p2))
! 390: {
! 391: icopyifstack(p2,z[1]);
! 392: if (!is_bigint(p2))
! 393: {
! 394: z[2] = lstoi(addssmod(itos((GEN)x[2]),itos((GEN)y[2]), p2[2]));
! 395: return z;
! 396: }
! 397: }
! 398: else
! 399: { p2 = mppgcd(p1,p2); z[1] = (long)p2; }
! 400: av=avma; (void)new_chunk((lgefint(p1)<<1) + lgefint(x[1]));
! 401: p1=addii((GEN)x[2],(GEN)y[2]); avma=av;
! 402: z[2]=lmodii(p1,p2); return z;
! 403:
! 404: case t_FRAC: case t_FRACN: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
! 405: (void)new_chunk(lgefint(p2)<<2);
! 406: p1 = mulii((GEN)y[1], mpinvmod((GEN)y[2],p2));
! 407: p1 = addii(modii(p1,p2), (GEN)x[2]); avma=(long)z;
! 408: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 409:
! 410: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 411: z[1]=ladd(x,(GEN)y[1]);
! 412: z[2]=lcopy((GEN)y[2]); return z;
! 413:
! 414: case t_PADIC:
! 415: l=avma; p1=cgetg(3,t_INTMOD);
! 416: p1[1]=x[1]; p1[2]=lgeti(lgefint(x[1]));
! 417: gaffect(y,p1); tetpil=avma;
! 418: return gerepile(l,tetpil,gadd(p1,x));
! 419:
! 420: case t_QUAD: z=cgetg(4,t_QUAD);
! 421: copyifstack(y[1], z[1]);
! 422: z[2]=ladd(x,(GEN)y[2]);
! 423: z[3]=lcopy((GEN)y[3]); return z;
! 424: }
! 425:
! 426: case t_FRAC: case t_FRACN:
! 427: switch (ty)
! 428: {
! 429: case t_FRAC: return addfrac(x,y);
! 430: case t_FRACN: z=cgetg(3,t_FRACN); l=avma;
! 431: p1=mulii((GEN)x[1],(GEN)y[2]);
! 432: p2=mulii((GEN)x[2],(GEN)y[1]);
! 433: tetpil=avma; z[1]=lpile(l,tetpil,addii(p1,p2));
! 434: z[2]=lmulii((GEN)x[2],(GEN)y[2]);
! 435: return z;
! 436:
! 437: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 438: z[1]=ladd((GEN)y[1],x);
! 439: z[2]=lcopy((GEN)y[2]); return z;
! 440:
! 441: case t_PADIC:
! 442: return gaddpex(x,y);
! 443:
! 444: case t_QUAD: z=cgetg(4,t_QUAD);
! 445: z[1]=lcopy((GEN)y[1]);
! 446: z[2]=ladd((GEN)y[2],x);
! 447: z[3]=lcopy((GEN)y[3]); return z;
! 448: }
! 449:
! 450: case t_COMPLEX:
! 451: switch(ty)
! 452: {
! 453: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 454: z[1]=ladd((GEN)x[1],(GEN)y[1]);
! 455: z[2]=ladd((GEN)x[2],(GEN)y[2]); return z;
! 456:
! 457: case t_PADIC:
! 458: if (krosg(-1,(GEN)y[2])== -1)
! 459: {
! 460: z=cgetg(3,t_COMPLEX);
! 461: z[1]=ladd((GEN)x[1],y);
! 462: z[2]=lcopy((GEN)x[2]); return z;
! 463: }
! 464: av=avma; l = signe(y[4])? precp(y): 1;
! 465: p1=cvtop(x,(GEN)y[2], l + valp(y)); tetpil=avma;
! 466: return gerepile(av,tetpil,gadd(p1,y));
! 467:
! 468: case t_QUAD:
! 469: lx=precision(x); if (!lx) err(gadderi,tx,ty);
! 470: if (gcmp0(y)) return gcopy(x);
! 471:
! 472: av=avma; i=gexpo(y)-gexpo(x);
! 473: if (i<=0) i=0; else i >>= TWOPOTBITS_IN_LONG;
! 474: p1=co8(y,lx+i); tetpil=avma;
! 475: return gerepile(av,tetpil,gadd(p1,x));
! 476: }
! 477:
! 478: case t_PADIC:
! 479: switch(ty)
! 480: {
! 481: case t_PADIC:
! 482: if (!egalii((GEN)x[2],(GEN)y[2])) err(gadderi,tx,ty);
! 483: return addpadic(x,y);
! 484:
! 485: case t_QUAD:
! 486: if (kro_quad((GEN)y[1],(GEN)x[2]) == -1)
! 487: {
! 488: z=cgetg(4,t_QUAD);
! 489: copyifstack(y[1], z[1]);
! 490: z[2]=ladd((GEN)y[2],x);
! 491: z[3]=lcopy((GEN)y[3]); return z;
! 492: }
! 493: av=avma; l = signe(x[4])? precp(x): 1;
! 494: p1=cvtop(y,(GEN)x[2],valp(x)+l); tetpil=avma;
! 495: return gerepile(av,tetpil,gadd(p1,x));
! 496: }
! 497:
! 498: case t_QUAD: z=cgetg(4,t_QUAD); k=x[1]; l=y[1];
! 499: if (!gegal((GEN)k,(GEN)l)) err(gadderi,tx,ty);
! 500: copyifstack(l, z[1]);
! 501: z[2]=ladd((GEN)x[2],(GEN)y[2]);
! 502: z[3]=ladd((GEN)x[3],(GEN)y[3]); return z;
! 503: }
! 504: err(bugparier,"gadd");
! 505: }
! 506:
! 507: /* here !isscalar(y) */
! 508: if ( (vx>vy && (!is_matvec_t(tx) || !is_matvec_t(ty)))
! 509: || (vx==vy && is_scalar_t(tx)) )
! 510: {
! 511: if (tx == t_POLMOD && vx == vy && ty != t_SER)
! 512: {
! 513: av = avma;
! 514: return gerepileupto(av, op_polmod(gadd, x, to_polmod(y,(GEN)x[1]), tx));
! 515: }
! 516:
! 517: switch(ty)
! 518: {
! 519: case t_POL: ly=lgef(y);
! 520: if (ly==2) return isexactzero(x)? zeropol(vy): scalarpol(x,vy);
! 521:
! 522: z = cgetg(ly,t_POL); z[1] = y[1];
! 523: z[2] = ladd(x,(GEN)y[2]);
! 524: for (i=3; i<ly; i++) z[i]=lcopy((GEN)y[i]);
! 525: return normalizepol_i(z, ly);
! 526:
! 527: case t_SER: l=valp(y); ly=lg(y);
! 528: if (l<3-ly) return gcopy(y);
! 529: if (l<0)
! 530: {
! 531: z=cgetg(ly,t_SER); z[1]=y[1];
! 532: for (i=2; i<=1-l; i++) z[i]=lcopy((GEN)y[i]);
! 533: for (i=3-l; i<ly; i++) z[i]=lcopy((GEN)y[i]);
! 534: z[2-l]=ladd(x,(GEN)y[2-l]); return z;
! 535: }
! 536: if (l>0)
! 537: {
! 538: if (gcmp0(x)) return gcopy(y);
! 539: if (gcmp0(y)) ly=2;
! 540:
! 541: ly += l; z=cgetg(ly,t_SER);
! 542: z[1]=evalsigne(1) | evalvalp(0) | evalvarn(vy);
! 543: for (i=3; i<=l+1; i++) z[i]=zero;
! 544: for ( ; i<ly; i++) z[i]=lcopy((GEN)y[i-l]);
! 545: z[2]=lcopy(x); return z;
! 546: }
! 547: av=avma; z=cgetg(ly,t_SER);
! 548: p1=signe(y)? gadd(x,(GEN)y[2]): x;
! 549: if (!isexactzero(p1))
! 550: {
! 551: z[1] = evalvalp(0) | evalvarn(vy);
! 552: if (signe(y))
! 553: {
! 554: z[1] |= evalsigne(1); z[2]=(long)p1;
! 555: for (i=3; i<ly; i++) z[i]=lcopy((GEN)y[i]);
! 556: }
! 557: return z;
! 558: }
! 559: avma=av; /* first coeff is 0 */
! 560: i=3; while (i<ly && gcmp0((GEN)y[i])) i++;
! 561: if (i==ly) return zeroser(vy,i-2);
! 562:
! 563: z=cgetg(ly-i+2,t_SER); z[1]=evalvalp(i-2)|evalvarn(vy)|evalsigne(1);
! 564: for (j=2; j<=ly-i+1; j++) z[j]=lcopy((GEN)y[j+i-2]);
! 565: return z;
! 566:
! 567: case t_RFRAC: return addscalrfrac(x,y);
! 568: case t_RFRACN: z=cgetg(3,t_RFRACN);
! 569: av=avma; p1=gmul(x,(GEN)y[2]); tetpil=avma;
! 570: z[1]=lpile(av,tetpil, gadd(p1,(GEN)y[1]));
! 571: z[2]=lcopy((GEN)y[2]); return z;
! 572:
! 573: case t_VEC: case t_COL: case t_MAT:
! 574: if (isexactzero(x)) return gcopy(y);
! 575: if (ty == t_MAT) return gaddmat(x,y);
! 576: /* fall through */
! 577: case t_QFR: case t_QFI: err(gadderf,tx,ty);
! 578: }
! 579: err(typeer,"addition");
! 580: }
! 581:
! 582: /* here !isscalar(x) && isscalar(y) && (vx=vy || ismatvec(x and y)) */
! 583: if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
! 584: switch(tx)
! 585: {
! 586: case t_POL:
! 587: switch (ty)
! 588: {
! 589: case t_POL:
! 590: lx = lgef(x); ly = lgef(y); if (lx < ly) swapspec(x,y, lx,ly);
! 591: z = cgetg(lx,t_POL); z[1] = x[1];
! 592: for (i=2; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 593: for ( ; i<lx; i++) z[i]=lcopy((GEN)x[i]);
! 594: (void)normalizepol_i(z, lx);
! 595: if (lgef(z) == 2) { avma = (long)(z + lx); z = zeropol(vx); }
! 596: return z;
! 597:
! 598: case t_SER:
! 599: if (gcmp0(x)) return gcopy(y);
! 600: ly = signe(y)? lg(y): 3;
! 601: i = ly+valp(y)-gval(x,vx);
! 602: if (i<3) return gcopy(y);
! 603:
! 604: p1=greffe(x,i,0); y=gadd(p1,y);
! 605: free(p1); return y;
! 606:
! 607: case t_RFRAC: return addscalrfrac(x,y);
! 608: case t_RFRACN: z=cgetg(3,t_RFRACN);
! 609: av=avma; p1=gmul(x,(GEN)y[2]); tetpil=avma;
! 610: z[1]=lpile(av,tetpil, gadd(p1,(GEN)y[1]));
! 611: z[2]=lcopy((GEN)y[2]); return z;
! 612:
! 613: case t_QFR: case t_QFI:
! 614: case t_VEC: case t_COL: case t_MAT: err(gadderf,tx,ty);
! 615: default: err(typeer,"addition");
! 616: }
! 617:
! 618: case t_SER:
! 619: switch(ty)
! 620: {
! 621: case t_SER:
! 622: l=valp(y)-valp(x);
! 623: if (l<0) { l= -l; p1=x; x=y; y=p1; }
! 624: if (gcmp0(x)) return gcopy(x);
! 625: lx = lg(x);
! 626: ly = signe(y)? lg(y): 2;
! 627: ly += l; if (lx<ly) ly=lx;
! 628: z=cgetg(ly,t_SER);
! 629: if (l)
! 630: {
! 631: if (l>=ly-2)
! 632: for (i=2; i<ly; i++) z[i]=lcopy((GEN)x[i]);
! 633: else
! 634: {
! 635: for (i=2; i<=l+1; i++) z[i]=lcopy((GEN)x[i]);
! 636: for ( ; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i-l]);
! 637: }
! 638: z[1]=x[1]; return z;
! 639: }
! 640: if (ly>2)
! 641: {
! 642: l=avma;
! 643: for (i=2; i<ly; i++)
! 644: {
! 645: p1 = gadd((GEN)x[i],(GEN)y[i]);
! 646: if (!isexactzero(p1))
! 647: {
! 648: l = i-2; stackdummy(z,l); z += l;
! 649: z[0]=evaltyp(t_SER) | evallg(ly-l);
! 650: z[1]=evalvalp(valp(x)+i-2) | evalsigne(1) | evalvarn(vx);
! 651: for (j=i+1; j<ly; j++)
! 652: z[j-l]=ladd((GEN)x[j],(GEN)y[j]);
! 653: z[2]=(long)p1; return z;
! 654: }
! 655: avma=l;
! 656: }
! 657: }
! 658: return zeroser(vx,ly-2+valp(y));
! 659:
! 660: case t_RFRAC: case t_RFRACN:
! 661: if (gcmp0(y)) return gcopy(x);
! 662:
! 663: l = valp(x)-gval(y,vy); l += gcmp0(x)? 3: lg(x);
! 664: if (l<3) return gcopy(x);
! 665:
! 666: av=avma; ty=typ(y[2]);
! 667: p1 = is_scalar_t(ty)? (GEN)y[2]: greffe((GEN)y[2],l,1);
! 668: p1 = gdiv((GEN)y[1], p1); tetpil=avma;
! 669: return gerepile(av,tetpil,gadd(p1,x));
! 670:
! 671: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
! 672: err(gadderf,tx,ty);
! 673:
! 674: default: err(typeer,"addition");
! 675: }
! 676:
! 677: case t_RFRAC:
! 678: if (!is_rfrac_t(ty)) err(gadderi,tx,ty);
! 679: return addrfrac(x,y);
! 680: case t_RFRACN:
! 681: if (!is_rfrac_t(ty)) err(gadderi,tx,ty);
! 682: z=cgetg(3,t_RFRACN); av=avma;
! 683: p1=gmul((GEN)x[1],(GEN)y[2]);
! 684: p2=gmul((GEN)x[2],(GEN)y[1]); tetpil=avma;
! 685: z[1]=lpile(av,tetpil, gadd(p1,p2));
! 686: z[2]=lmul((GEN)x[2],(GEN)y[2]); return z;
! 687:
! 688: case t_VEC: case t_COL: case t_MAT:
! 689: lx = lg(x); ly = lg(y);
! 690: if (lx!=ly || tx!=ty) err(gadderi,tx,ty);
! 691: z=cgetg(ly,ty);
! 692: for (i=1; i<ly; i++)
! 693: z[i]=ladd((GEN)x[i],(GEN)y[i]);
! 694: return z;
! 695:
! 696: case t_QFR: case t_QFI: err(gadderf,tx,ty);
! 697: }
! 698: err(typeer,"addition");
! 699: return NULL; /* not reached */
! 700: }
! 701:
! 702: /********************************************************************/
! 703: /** **/
! 704: /** MULTIPLICATION **/
! 705: /** **/
! 706: /********************************************************************/
! 707: #define fix_frac(z) if (signe(z[2])<0)\
! 708: {\
! 709: setsigne(z[1],-signe(z[1]));\
! 710: setsigne(z[2],1);\
! 711: }
! 712:
! 713: /* assume z[1] was created last */
! 714: #define fix_frac_if_int(z) if (is_pm1(z[2]))\
! 715: z = gerepileupto((long)(z+3), (GEN)z[1]);
! 716:
! 717: /* assume z[1] was created last */
! 718: #define fix_frac_if_int_GC(z,tetpil) { if (is_pm1(z[2]))\
! 719: z = gerepileupto((long)(z+3), (GEN)z[1]);\
! 720: else\
! 721: gerepilemanyvec((long)z, tetpil, z+1, 2); }
! 722:
! 723: GEN
! 724: fix_rfrac_if_pol(GEN x, GEN y)
! 725: {
! 726: if (gcmp1(y)) return x;
! 727: if (typ(y) != t_POL)
! 728: {
! 729: if (typ(x) != t_POL || gvar2(y) > varn(x))
! 730: return gdiv(x,y);
! 731: }
! 732: else if (varn(y) > varn(x)) return gdiv(x,y);
! 733: return NULL;
! 734: }
! 735:
! 736: static long
! 737: mingvar(GEN x, GEN y)
! 738: {
! 739: long i = gvar(x);
! 740: long j = gvar(y);
! 741: return min(i,j);
! 742: }
! 743:
! 744: GEN
! 745: mulscalrfrac(GEN x, GEN y)
! 746: {
! 747: GEN p1,z,y1,y2,cx,cy1,cy2;
! 748: long tetpil,tx;
! 749:
! 750: if (gcmp0(x)) return gcopy(x);
! 751:
! 752: y1=(GEN)y[1]; if (gcmp0(y1)) return gcopy(y1);
! 753: y2=(GEN)y[2]; tx = typ(x);
! 754: z = cgetg(3, t_RFRAC);
! 755: if (is_const_t(tx) || varn(x) > mingvar(y1,y2)) { cx = x; x = gun; }
! 756: else
! 757: {
! 758: p1 = ggcd(x,y2); if (isnonscalar(p1)) { x=gdeuc(x,p1); y2=gdeuc(y2,p1); }
! 759: cx = content(x); if (!gcmp1(cx)) x = gdiv(x,cx);
! 760: }
! 761: cy1 = content(y1); if (!gcmp1(cy1)) y1 = gdiv(y1,cy1);
! 762: cy2 = content(y2); if (!gcmp1(cy2)) y2 = gdiv(y2,cy2);
! 763: if (x != gun) y1 = gmul(y1,x);
! 764: x = gdiv(gmul(cx,cy1), cy2);
! 765: cy1 = numer(x);
! 766: cy2 = denom(x); tetpil = avma;
! 767: z[2] = lmul(y2, cy2);
! 768: z[1] = lmul(y1, cy1);
! 769: p1 = fix_rfrac_if_pol((GEN)z[1],(GEN)z[2]);
! 770: if (p1) return gerepileupto((long)(z+3), p1);
! 771: gerepilemanyvec((long)z,tetpil,z+1,2); return z;
! 772: }
! 773:
! 774: static GEN
! 775: mulrfrac(GEN x, GEN y)
! 776: {
! 777: GEN z = cgetg(3,t_RFRAC), p1;
! 778: GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
! 779: GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
! 780: long tetpil;
! 781:
! 782: p1 = ggcd(x1, y2); if (!gcmp1(p1)) { x1 = gdiv(x1,p1); y2 = gdiv(y2,p1); }
! 783: p1 = ggcd(x2, y1); if (!gcmp1(p1)) { x2 = gdiv(x2,p1); y1 = gdiv(y1,p1); }
! 784: tetpil = avma;
! 785: z[2] = lmul(x2,y2);
! 786: z[1] = lmul(x1,y1);
! 787: p1 = fix_rfrac_if_pol((GEN)z[1],(GEN)z[2]);
! 788: if (p1) return gerepileupto((long)(z+3), p1);
! 789: gerepilemanyvec((long)z,tetpil,z+1,2); return z;
! 790: }
! 791:
! 792: GEN
! 793: to_Kronecker(GEN P, GEN Q)
! 794: {
! 795: /* P(X) = sum Mod(.,Q(Y)) * X^i, lift then set X := Y^(2n-1) */
! 796: long i,j,k,l, lx = lgef(P), N = ((lgef(Q)-3)<<1) + 1, vQ = varn(Q);
! 797: GEN p1, y = cgetg((N-2)*(lx-2) + 2, t_POL);
! 798: for (k=i=2; i<lx; i++)
! 799: {
! 800: p1 = (GEN)P[i];
! 801: if ((l=typ(p1)) == t_POLMOD) { p1 = (GEN)p1[2]; l = typ(p1); }
! 802: if (is_scalar_t(l) || varn(p1)<vQ) { y[k++] = (long)p1; j = 3; }
! 803: else
! 804: {
! 805: l = lgef(p1);
! 806: for (j=2; j < l; j++) y[k++] = p1[j];
! 807: }
! 808: if (i == lx-1) break;
! 809: for ( ; j < N; j++) y[k++] = zero;
! 810: }
! 811: y[1] = evalsigne(1)|evalvarn(vQ)|evallgef(k);
! 812: return y;
! 813: }
! 814:
! 815: int
! 816: ff_poltype(GEN *x, GEN *p, GEN *pol)
! 817: {
! 818: GEN Q, P = *x, pr,p1,p2,y;
! 819: long i, lx;
! 820:
! 821: if (!signe(P)) return 0;
! 822: lx = lgef(P); Q = *pol;
! 823: for (i=2; i<lx; i++)
! 824: {
! 825: p1 = (GEN)P[i]; if (typ(p1) != t_POLMOD) { Q = NULL; break; }
! 826: p2 = (GEN)p1[1];
! 827: if (Q==NULL) Q = p2;
! 828: else if (p2 != Q)
! 829: {
! 830: if (!gegal(p2, Q))
! 831: {
! 832: if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
! 833: return 0;
! 834: }
! 835: if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
! 836: }
! 837: }
! 838: if (Q) {
! 839: *x = P = to_Kronecker(P, Q);
! 840: *pol = Q; lx = lgef(P);
! 841: }
! 842: pr = *p; y = cgetg(lx, t_POL);
! 843: for (i=lx-1; i>1; i--)
! 844: {
! 845: p1 = (GEN)P[i];
! 846: switch(typ(p1))
! 847: {
! 848: case t_INTMOD: break;
! 849: case t_INT:
! 850: if (*p) p1 = modii(p1, *p);
! 851: y[i] = (long)p1; continue;
! 852: default:
! 853: return (Q && !pr)? 1: 0;
! 854: }
! 855: p2 = (GEN)p1[1];
! 856: if (pr==NULL) pr = p2;
! 857: else if (p2 != pr)
! 858: {
! 859: if (!egalii(p2, pr))
! 860: {
! 861: if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
! 862: return 0;
! 863: }
! 864: if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
! 865: }
! 866: y[i] = p1[2];
! 867: }
! 868: y[1] = evalsigne(1)|evalvarn(varn(P))|evallgef(lx);
! 869: *x = y; *p = pr; return (Q || pr);
! 870: }
! 871:
! 872: GEN
! 873: gmul(GEN x, GEN y)
! 874: {
! 875: long tx,ty,lx,ly,vx,vy,i,j,k,l,av,tetpil;
! 876: GEN z,p1,p2,p3,p4;
! 877:
! 878: if (x == y) return gsqr(x);
! 879: tx=typ(x); ty=typ(y);
! 880: if (is_const_t(tx) && is_const_t(ty))
! 881: {
! 882: if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
! 883: }
! 884: else
! 885: {
! 886: vx=gvar(x); vy=gvar(y);
! 887: if (!is_matvec_t(ty))
! 888: if (is_matvec_t(tx) || vx<vy || (vx==vy && tx>ty))
! 889: {
! 890: p1=x; x=y; y=p1;
! 891: i=tx; tx=ty; ty=i;
! 892: i=vx; vx=vy; vy=i;
! 893: }
! 894: if (ty==t_POLMOD) return op_polmod(gmul,x,y,tx);
! 895: }
! 896:
! 897: if (is_scalar_t(ty))
! 898: {
! 899: switch(tx)
! 900: {
! 901: case t_INT:
! 902: switch(ty)
! 903: {
! 904: case t_INT: return mulii(x,y);
! 905: case t_REAL: return mulir(x,y);
! 906:
! 907: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
! 908: (void)new_chunk(lgefint(p2)<<2);
! 909: p1=mulii(modii(x,p2),(GEN)y[2]); avma=(long)z;
! 910: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 911:
! 912: case t_FRAC:
! 913: if (!signe(x)) return gzero;
! 914: z=cgetg(3,t_FRAC);
! 915: p1 = mppgcd(x,(GEN)y[2]);
! 916: if (is_pm1(p1))
! 917: {
! 918: avma = (long)z;
! 919: z[2] = licopy((GEN)y[2]);
! 920: z[1] = lmulii((GEN)y[1], x);
! 921: }
! 922: else
! 923: {
! 924: x = divii(x,p1); tetpil = avma;
! 925: z[2] = ldivii((GEN)y[2], p1);
! 926: z[1] = lmulii((GEN)y[1], x);
! 927: fix_frac_if_int_GC(z,tetpil);
! 928: }
! 929: return z;
! 930:
! 931: case t_FRACN: z=cgetg(3,t_FRACN);
! 932: z[1]=lmulii(x,(GEN)y[1]);
! 933: z[2]=licopy((GEN)y[2]); return z;
! 934:
! 935: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 936: z[1]=lmul(x,(GEN)y[1]);
! 937: z[2]=lmul(x,(GEN)y[2]); return z;
! 938:
! 939: case t_PADIC:
! 940: if (!signe(x)) return gzero;
! 941: l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
! 942: return gerepile(l,tetpil,gmul(p1,y));
! 943:
! 944: case t_QUAD: z=cgetg(4,t_QUAD);
! 945: copyifstack(y[1], z[1]);
! 946: z[2]=lmul(x,(GEN)y[2]);
! 947: z[3]=lmul(x,(GEN)y[3]); return z;
! 948: }
! 949:
! 950: case t_REAL:
! 951: switch(ty)
! 952: {
! 953: case t_REAL: return mulrr(x,y);
! 954:
! 955: case t_FRAC: case t_FRACN:
! 956: l=avma; p1=cgetr(lg(x)); tetpil=avma; gaffect(y,p1);
! 957: p2=mulrr(p1,x); return gerepile(l,tetpil,p2);
! 958:
! 959: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 960: z[1]=lmul(x,(GEN)y[1]);
! 961: z[2]=lmul(x,(GEN)y[2]); return z;
! 962:
! 963: case t_QUAD:
! 964: l=avma; p1=co8(y,lg(x)); tetpil=avma;
! 965: return gerepile(l,tetpil,gmul(p1,x));
! 966:
! 967: case t_INTMOD: err(gmulerf,tx,ty);
! 968: }
! 969:
! 970: case t_INTMOD:
! 971: switch(ty)
! 972: {
! 973: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
! 974: if (p1==p2 || egalii(p1,p2))
! 975: {
! 976: icopyifstack(p2,z[1]);
! 977: if (!is_bigint(p2))
! 978: {
! 979: z[2] = lstoi(mulssmod(itos((GEN)x[2]),itos((GEN)y[2]), p2[2]));
! 980: return z;
! 981: }
! 982: }
! 983: else
! 984: { p2 = mppgcd(p1,p2); z[1] = (long)p2; }
! 985: av=avma;
! 986: (void)new_chunk(lgefint(x[1]) + (lgefint(p1)<<1));
! 987: p1=mulii((GEN)x[2],(GEN)y[2]); avma=av;
! 988: z[2]=lmodii(p1,p2); return z;
! 989:
! 990: case t_FRAC: case t_FRACN: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
! 991: (void)new_chunk(lgefint(p2)<<2);
! 992: p1 = mulii((GEN)y[1], mpinvmod((GEN)y[2],p2));
! 993: p1 = mulii(modii(p1,p2),(GEN)x[2]); avma=(long)z;
! 994: z[2] = lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 995:
! 996: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 997: z[1]=lmul(x,(GEN)y[1]);
! 998: z[2]=lmul(x,(GEN)y[2]); return z;
! 999:
! 1000: case t_PADIC:
! 1001: l=avma; p1=cgetg(3,t_INTMOD);
! 1002: p1[1]=x[1]; p1[2]=lgeti(lg(x[1]));
! 1003: gaffect(y,p1); tetpil=avma;
! 1004: return gerepile(l,tetpil,gmul(x,p1));
! 1005:
! 1006: case t_QUAD: z=cgetg(4,t_QUAD);
! 1007: copyifstack(y[1], z[1]);
! 1008: z[2]=lmul(x,(GEN)y[2]);
! 1009: z[3]=lmul(x,(GEN)y[3]); return z;
! 1010: }
! 1011:
! 1012: case t_FRAC: case t_FRACN:
! 1013: switch(ty)
! 1014: {
! 1015: case t_FRAC:
! 1016: {
! 1017: GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
! 1018: GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
! 1019: z=cgetg(3,t_FRAC);
! 1020: p1 = mppgcd(x1, y2);
! 1021: if (!is_pm1(p1)) { x1 = divii(x1,p1); y2 = divii(y2,p1); }
! 1022: p1 = mppgcd(x2, y1);
! 1023: if (!is_pm1(p1)) { x2 = divii(x2,p1); y1 = divii(y1,p1); }
! 1024: tetpil = avma;
! 1025: z[2] = lmulii(x2,y2);
! 1026: z[1] = lmulii(x1,y1);
! 1027: fix_frac_if_int_GC(z,tetpil); return z;
! 1028: }
! 1029: case t_FRACN: z=cgetg(3,t_FRACN);
! 1030: z[1]=lmulii((GEN)x[1],(GEN)y[1]);
! 1031: z[2]=lmulii((GEN)x[2],(GEN)y[2]); return z;
! 1032:
! 1033: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 1034: z[1]=lmul((GEN)y[1],x);
! 1035: z[2]=lmul((GEN)y[2],x); return z;
! 1036:
! 1037: case t_PADIC:
! 1038: if (!signe(x[1])) return gzero;
! 1039: l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
! 1040: return gerepile(l,tetpil,gmul(p1,y));
! 1041:
! 1042: case t_QUAD: z=cgetg(4,t_QUAD);
! 1043: copyifstack(y[1], z[1]);
! 1044: z[2]=lmul((GEN)y[2],x);
! 1045: z[3]=lmul((GEN)y[3],x); return z;
! 1046: }
! 1047:
! 1048: case t_COMPLEX:
! 1049: switch(ty)
! 1050: {
! 1051: case t_COMPLEX: z=cgetg(3,t_COMPLEX); l=avma;
! 1052: p1=gmul((GEN)x[1],(GEN)y[1]);
! 1053: p2=gmul((GEN)x[2],(GEN)y[2]);
! 1054: x=gadd((GEN)x[1],(GEN)x[2]);
! 1055: y=gadd((GEN)y[1],(GEN)y[2]);
! 1056: y=gmul(x,y); x=gadd(p1,p2);
! 1057: tetpil=avma; z[1]=lsub(p1,p2); z[2]=lsub(y,x);
! 1058: gerepilemanyvec(l,tetpil,z+1,2); return z;
! 1059:
! 1060: case t_PADIC:
! 1061: if (krosg(-1,(GEN)y[2]))
! 1062: {
! 1063: z=cgetg(3,t_COMPLEX);
! 1064: z[1]=lmul((GEN)x[1],y);
! 1065: z[2]=lmul((GEN)x[2],y); return z;
! 1066: }
! 1067: av=avma;
! 1068: if (signe(y[4])) l=precp(y);
! 1069: else
! 1070: {
! 1071: l=valp(y)+1; if (l<=0) l=1;
! 1072: }
! 1073: p1=cvtop(x,(GEN)y[2],l); tetpil=avma;
! 1074: return gerepile(av,tetpil,gmul(p1,y));
! 1075:
! 1076: case t_QUAD:
! 1077: lx=precision(x); if (!lx) err(gmuleri,tx,ty);
! 1078: l=avma; p1=co8(y,lx); tetpil=avma;
! 1079: return gerepile(l,tetpil,gmul(p1,x));
! 1080: }
! 1081:
! 1082: case t_PADIC:
! 1083: switch(ty)
! 1084: {
! 1085: case t_PADIC:
! 1086: if (!egalii((GEN)x[2],(GEN)y[2])) err(gmuleri,tx,ty);
! 1087: l = valp(x)+valp(y);
! 1088: if (!signe(x[4])) { z=gcopy(x); setvalp(z,l); return z; }
! 1089: if (!signe(y[4])) { z=gcopy(y); setvalp(z,l); return z; }
! 1090:
! 1091: p1 = (precp(x) > precp(y))? y: x;
! 1092: z=cgetp(p1); setvalp(z,l); av=avma;
! 1093: modiiz(mulii((GEN)x[4],(GEN)y[4]),(GEN)p1[3],(GEN)z[4]);
! 1094: avma=av; return z;
! 1095:
! 1096: case t_QUAD:
! 1097: if (kro_quad((GEN)y[1],(GEN)x[2])== -1)
! 1098: {
! 1099: z=cgetg(4,t_QUAD);
! 1100: copyifstack(y[1], z[1]);
! 1101: z[2]=lmul((GEN)y[2],x);
! 1102: z[3]=lmul((GEN)y[3],x); return z;
! 1103: }
! 1104: l = signe(x[4])? precp(x): valp(x)+1;
! 1105: av=avma; p1=cvtop(y,(GEN)x[2],l); tetpil=avma;
! 1106: return gerepile(av,tetpil,gmul(p1,x));
! 1107: }
! 1108:
! 1109: case t_QUAD: z=cgetg(4,t_QUAD);
! 1110: p1=(GEN)x[1]; p2=(GEN)y[1];
! 1111: if (!gegal(p1,p2)) err(gmuleri,tx,ty);
! 1112:
! 1113: copyifstack(p2, z[1]); l=avma;
! 1114: p2=gmul((GEN)x[2],(GEN)y[2]);
! 1115: p3=gmul((GEN)x[3],(GEN)y[3]);
! 1116: p4=gmul(gneg_i((GEN)p1[2]),p3);
! 1117:
! 1118: if (gcmp0((GEN)p1[3]))
! 1119: {
! 1120: tetpil=avma;
! 1121: z[2]=lpile(l,tetpil,gadd(p4,p2)); l=avma;
! 1122: p2=gmul((GEN)x[2],(GEN)y[3]);
! 1123: p3=gmul((GEN)x[3],(GEN)y[2]); tetpil=avma;
! 1124: z[3]=lpile(l,tetpil,gadd(p2,p3)); return z;
! 1125: }
! 1126:
! 1127: p1 = gadd(gmul((GEN)x[2],(GEN)y[3]), gmul((GEN)x[3],(GEN)y[2]));
! 1128: tetpil=avma;
! 1129: z[2]=ladd(p2,p4);
! 1130: z[3]=ladd(p1,p3);
! 1131: gerepilemanyvec(l,tetpil,z+2,2); return z;
! 1132: }
! 1133: err(bugparier,"multiplication");
! 1134: }
! 1135:
! 1136: /* here !isscalar(y) */
! 1137: if (is_matvec_t(ty))
! 1138: {
! 1139: ly=lg(y);
! 1140: if (!is_matvec_t(tx))
! 1141: {
! 1142: z=cgetg(ly,ty);
! 1143: for (i=1; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
! 1144: return z;
! 1145: }
! 1146: lx=lg(x);
! 1147:
! 1148: switch(tx)
! 1149: {
! 1150: case t_VEC:
! 1151: switch(ty)
! 1152: {
! 1153: case t_COL:
! 1154: if (lx!=ly) err(gmuleri,tx,ty);
! 1155: z=gzero; l=avma;
! 1156: for (i=1; i<lx; i++)
! 1157: {
! 1158: p1=gmul((GEN)x[i],(GEN)y[i]);
! 1159: tetpil=avma; z=gadd(z,p1);
! 1160: }
! 1161: return gerepile(l,tetpil,z);
! 1162:
! 1163: case t_MAT:
! 1164: if (ly==1) return cgetg(1,t_VEC);
! 1165: l=lg(y[1]); if (lx!=l) err(gmuleri,tx,ty);
! 1166:
! 1167: z=cgetg(ly,tx);
! 1168: for (i=1; i<ly; i++)
! 1169: {
! 1170: p1=gzero; av=avma;
! 1171: for (j=1; j<lx; j++)
! 1172: {
! 1173: p2=gmul((GEN)x[j],gcoeff(y,j,i));
! 1174: tetpil=avma; p1=gadd(p1,p2);
! 1175: }
! 1176: z[i]=lpile(av,tetpil,p1);
! 1177: }
! 1178: return z;
! 1179:
! 1180: case t_VEC: err(gmulerf,tx,ty);
! 1181: default: err(typeer,"multiplication");
! 1182: }
! 1183:
! 1184: case t_COL:
! 1185: switch(ty)
! 1186: {
! 1187: case t_VEC:
! 1188: z=cgetg(ly,t_MAT);
! 1189: for (i=1; i<ly; i++)
! 1190: {
! 1191: p1 = gmul((GEN)y[i],x);
! 1192: if (typ(p1) != t_COL) err(gmuleri,tx,ty);
! 1193: z[i]=(long)p1;
! 1194: }
! 1195: return z;
! 1196:
! 1197: case t_MAT:
! 1198: if (ly!=1 && lg(y[1])!=2) err(gmuleri,tx,ty);
! 1199:
! 1200: z=cgetg(ly,t_MAT);
! 1201: for (i=1; i<ly; i++) z[i]=lmul(gcoeff(y,1,i),x);
! 1202: return z;
! 1203:
! 1204: case t_COL: err(gmulerf,tx,ty);
! 1205: default: err(typeer,"multiplication");
! 1206: }
! 1207:
! 1208: case t_MAT:
! 1209: switch(ty)
! 1210: {
! 1211: case t_VEC:
! 1212: if (lx!=2) err(gmuleri,tx,ty);
! 1213: z=cgetg(ly,t_MAT);
! 1214: for (i=1; i<ly; i++) z[i]=lmul((GEN)y[i],(GEN)x[1]);
! 1215: return z;
! 1216:
! 1217: case t_COL:
! 1218: if (lx!=ly) err(gmuleri,tx,ty);
! 1219: if (lx==1) return gcopy(y);
! 1220:
! 1221: lx=lg(x[1]); z=cgetg(lx,t_COL);
! 1222: for (i=1; i<lx; i++)
! 1223: {
! 1224: p1=gzero; l=avma;
! 1225: for (j=1; j<ly; j++)
! 1226: {
! 1227: p2=gmul(gcoeff(x,i,j),(GEN)y[j]);
! 1228: tetpil=avma; p1=gadd(p1,p2);
! 1229: }
! 1230: z[i]=lpile(l,tetpil,p1);
! 1231: }
! 1232: return z;
! 1233:
! 1234: case t_MAT:
! 1235: if (ly==1) return cgetg(ly,t_MAT);
! 1236: if (lx != lg(y[1])) err(gmuleri,tx,ty);
! 1237: z=cgetg(ly,t_MAT);
! 1238: if (lx==1)
! 1239: {
! 1240: for (i=1; i<ly; i++) z[i]=lgetg(1,t_COL);
! 1241: return z;
! 1242: }
! 1243: l=lg(x[1]);
! 1244: for (j=1; j<ly; j++)
! 1245: {
! 1246: z[j] = lgetg(l,t_COL);
! 1247: for (i=1; i<l; i++)
! 1248: {
! 1249: p1=gzero; av=avma;
! 1250: for (k=1; k<lx; k++)
! 1251: {
! 1252: p2=gmul(gcoeff(x,i,k),gcoeff(y,k,j));
! 1253: tetpil=avma; p1=gadd(p1,p2);
! 1254: }
! 1255: coeff(z,i,j)=lpile(av,tetpil,p1);
! 1256: }
! 1257: }
! 1258: return z;
! 1259: }
! 1260: }
! 1261: err(bugparier,"multiplication");
! 1262: }
! 1263: /* now !ismatvec(x and y) */
! 1264:
! 1265: if (vx>vy || (vx==vy && is_scalar_t(tx)))
! 1266: {
! 1267: if (isexactzero(x)) return zeropol(vy);
! 1268: if (tx == t_INT && is_pm1(x))
! 1269: return (signe(x)>0) ? gcopy(y): gneg(y);
! 1270: if (tx == t_POLMOD && vx == vy && ty != t_SER)
! 1271: {
! 1272: av = avma;
! 1273: return gerepileupto(av, op_polmod(gmul, x, to_polmod(y,(GEN)x[1]), tx));
! 1274: }
! 1275: switch(ty)
! 1276: {
! 1277: case t_POL:
! 1278: if (isexactzero(y)) return zeropol(vy);
! 1279: ly = lgef(y); z = cgetg(ly,t_POL); z[1]=y[1];
! 1280: for (i=2; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
! 1281: return normalizepol_i(z,ly);
! 1282:
! 1283: case t_SER:
! 1284: if (!signe(y)) return gcopy(y);
! 1285: ly=lg(y); z=cgetg(ly,t_SER);
! 1286: for (i=2; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
! 1287: z[1]=y[1]; return normalize(z);
! 1288:
! 1289: case t_RFRAC: return mulscalrfrac(x,y);
! 1290: case t_RFRACN: av=avma; z=cgetg(3,t_RFRACN);
! 1291: z[1]=lmul(x,(GEN)y[1]);
! 1292: z[2]=lcopy((GEN)y[2]); return z;
! 1293: }
! 1294: err(typeer,"multiplication");
! 1295: }
! 1296:
! 1297: if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
! 1298: switch(tx)
! 1299: {
! 1300: case t_POL:
! 1301: switch (ty)
! 1302: {
! 1303: case t_POL:
! 1304: {
! 1305: GEN a = x,b = y, p = NULL, pol = NULL;
! 1306: long av = avma;
! 1307: if (ff_poltype(&x,&p,&pol) && ff_poltype(&y,&p,&pol))
! 1308: {
! 1309: /* fprintferr("HUM"); */
! 1310: if (pol && varn(x) != varn(y))
! 1311: x = to_Kronecker(x,pol);
! 1312: z = quickmul(x+2, y+2, lgef(x)-2, lgef(y)-2);
! 1313: if (p) z = Fp_pol(z,p);
! 1314: if (pol) z = from_Kronecker(z,pol);
! 1315: z = gerepileupto(av, z);
! 1316: }
! 1317: else
! 1318: {
! 1319: avma = av;
! 1320: z = quickmul(a+2, b+2, lgef(a)-2, lgef(b)-2);
! 1321: }
! 1322: setvarn(z,vx); return z;
! 1323: }
! 1324: case t_SER:
! 1325: if (gcmp0(x)) return zeropol(vx);
! 1326: if (gcmp0(y)) return zeroser(vx, valp(y)+gval(x,vx));
! 1327: p1=greffe(x,lg(y),0); p2=gmul(p1,y);
! 1328: free(p1); return p2;
! 1329:
! 1330: case t_RFRAC: return mulscalrfrac(x,y);
! 1331: case t_RFRACN: av=avma; z=cgetg(3,t_RFRACN);
! 1332: z[1]=lmul(x,(GEN)y[1]);
! 1333: z[2]=lcopy((GEN)y[2]); return z;
! 1334:
! 1335: default: err(typeer,"multiplication");
! 1336: }
! 1337:
! 1338: case t_SER:
! 1339: switch (ty)
! 1340: {
! 1341: case t_SER:
! 1342: if (gcmp0(x) || gcmp0(y)) return zeroser(vx, valp(x)+valp(y));
! 1343: lx=lg(x); ly=lg(y);
! 1344: if (lx>ly) { k=ly; ly=lx; lx=k; p1=y; y=x; x=p1; }
! 1345: z = cgetg(lx,t_SER);
! 1346: z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
! 1347: x += 2; y += 2; z += 2; lx -= 3;
! 1348: p2 = (GEN)gpmalloc((lx+1)*sizeof(long));
! 1349: for (i=0; i<=lx; i++)
! 1350: {
! 1351: p2[i] = !isexactzero((GEN)y[i]);
! 1352: p1 = gzero; av = avma;
! 1353: for (j=0; j<=i; j++)
! 1354: if (p2[j])
! 1355: p1 = gadd(p1, gmul((GEN)y[j],(GEN)x[i-j]));
! 1356: z[i] = lpileupto(av,p1);
! 1357: }
! 1358: z -= 2; /* back to normalcy */
! 1359: free(p2); return normalize(z);
! 1360:
! 1361: case t_RFRAC: case t_RFRACN:
! 1362: if (gcmp0(y)) return zeropol(vx);
! 1363: if (gcmp0(x)) return zeroser(vx, valp(x)+gval(y,vx));
! 1364: l=avma; p1=gmul((GEN)y[1],x); tetpil=avma;
! 1365: return gerepile(l,tetpil,gdiv(p1,(GEN)y[2]));
! 1366:
! 1367: default: err(typeer,"multiplication");
! 1368: }
! 1369:
! 1370: /* (tx,ty) == t_RFRAC <==> ty == t_RFRAC */
! 1371: case t_RFRAC: return mulrfrac(x,y);
! 1372: case t_RFRACN:
! 1373: if (!is_rfrac_t(ty)) err(typeer,"multiplication");
! 1374: av=avma; z=cgetg(3,ty);
! 1375: z[1]=lmul((GEN)x[1],(GEN)y[1]);
! 1376: z[2]=lmul((GEN)x[2],(GEN)y[2]); return z;
! 1377: }
! 1378: if (tx==ty)
! 1379: switch(tx)
! 1380: {
! 1381: case t_QFI: return compimag(x,y);
! 1382: case t_QFR: return compreal(x,y);
! 1383: }
! 1384: err(typeer,"multiplication");
! 1385: return NULL; /* not reached */
! 1386: }
! 1387:
! 1388: /********************************************************************/
! 1389: /** **/
! 1390: /** DIVISION **/
! 1391: /** **/
! 1392: /********************************************************************/
! 1393:
! 1394: static
! 1395: GEN divrfracscal(GEN x, GEN y)
! 1396: {
! 1397: long Y[3]; Y[1]=un; Y[2]=(long)y;
! 1398: return mulrfrac(x,Y);
! 1399: }
! 1400:
! 1401: static
! 1402: GEN divscalrfrac(GEN x, GEN y)
! 1403: {
! 1404: long Y[3]; Y[1]=y[2]; Y[2]=y[1];
! 1405: return mulscalrfrac(x,Y);
! 1406: }
! 1407:
! 1408: static
! 1409: GEN divrfrac(GEN x, GEN y)
! 1410: {
! 1411: long Y[3]; Y[1]=y[2]; Y[2]=y[1];
! 1412: return mulrfrac(x,Y);
! 1413: }
! 1414:
! 1415: GEN
! 1416: gdiv(GEN x, GEN y)
! 1417: {
! 1418: long tx = typ(x), ty = typ(y), lx,ly,vx,vy,i,j,k,l,av,tetpil;
! 1419: GEN z,p1,p2,p3;
! 1420:
! 1421: if (tx==t_INT && is_const_t(ty))
! 1422: {
! 1423: switch (signe(x))
! 1424: {
! 1425: case 0:
! 1426: if (gcmp0(y)) err(gdiver2);
! 1427: if (ty != t_INTMOD) return gzero;
! 1428: z = cgetg(3,t_INTMOD); icopyifstack(y[1],z[1]); z[2]=zero;
! 1429: return z;
! 1430: case 1:
! 1431: if (is_pm1(x)) return ginv(y);
! 1432: break;
! 1433: case -1:
! 1434: if (is_pm1(x)) { av = avma; return gerepileupto(av, ginv(gneg(y))); }
! 1435: }
! 1436: switch(ty)
! 1437: {
! 1438: case t_INT:
! 1439: av=avma; z=dvmdii(x,y,&p1);
! 1440: if (p1==gzero) return z;
! 1441: (void)new_chunk((lgefint(x) + lgefint(y)) << 2);
! 1442: p1 = mppgcd(y,p1);
! 1443: avma=av; z=cgetg(3,t_FRAC);
! 1444: if (is_pm1(p1)) {
! 1445: z[1]=licopy(x);
! 1446: z[2]=licopy(y);
! 1447: }
! 1448: else {
! 1449: z[1]=ldivii(x,p1);
! 1450: z[2]=ldivii(y,p1);
! 1451: }
! 1452: fix_frac(z); return z;
! 1453:
! 1454: case t_REAL:
! 1455: return divir(x,y);
! 1456:
! 1457: case t_INTMOD:
! 1458: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
! 1459: (void)new_chunk(lgefint(p2)<<2);
! 1460: p1=mulii(modii(x,p2), mpinvmod((GEN)y[2],p2)); avma=(long)z;
! 1461: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 1462:
! 1463: case t_FRAC:
! 1464: z=cgetg(3,t_FRAC);
! 1465: p1 = mppgcd(x,(GEN)y[1]);
! 1466: if (is_pm1(p1))
! 1467: {
! 1468: avma = (long)z; tetpil = 0;
! 1469: z[2] = licopy((GEN)y[1]);
! 1470: }
! 1471: else
! 1472: {
! 1473: x = divii(x,p1); tetpil = avma;
! 1474: z[2] = ldivii((GEN)y[1], p1);
! 1475: }
! 1476: z[1] = lmulii((GEN)y[2], x);
! 1477: fix_frac(z);
! 1478: if (tetpil)
! 1479: { fix_frac_if_int_GC(z,tetpil); }
! 1480: else
! 1481: fix_frac_if_int(z);
! 1482: return z;
! 1483:
! 1484: case t_FRACN:
! 1485: z=cgetg(3,t_FRACN);
! 1486: z[1]=lmulii((GEN)y[2], x);
! 1487: z[2]=licopy((GEN)y[1]);
! 1488: fix_frac(z); return z;
! 1489:
! 1490: case t_PADIC:
! 1491: l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
! 1492: return gerepile(l,tetpil,gdiv(p1,y));
! 1493:
! 1494: case t_COMPLEX: case t_QUAD:
! 1495: l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
! 1496: return gerepile(l,tetpil,gdiv(p2,p1));
! 1497: }
! 1498: }
! 1499: if (gcmp0(y)) err(gdiver2);
! 1500:
! 1501: if (is_const_t(tx) && is_const_t(ty))
! 1502: {
! 1503: switch(tx)
! 1504: {
! 1505: case t_REAL:
! 1506: switch(ty)
! 1507: {
! 1508: case t_INT:
! 1509: return divri(x,y);
! 1510:
! 1511: case t_REAL:
! 1512: return divrr(x,y);
! 1513:
! 1514: case t_FRAC: case t_FRACN:
! 1515: l=avma; p1=cgetg(lg(x),t_REAL); gaffect(y,p1);
! 1516: return gerepile(l,(long)p1,divrr(x,p1));
! 1517:
! 1518: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 1519: l=avma; p1=gnorm(y);
! 1520: p2=gmul(x,(GEN)y[1]);
! 1521: p3=gmul(x,(GEN)y[2]);
! 1522: if (!gcmp0(p3)) p3 = gneg_i(p3);
! 1523: tetpil=avma;
! 1524: z[1]=ldiv(p2,p1);
! 1525: z[2]=ldiv(p3,p1);
! 1526: gerepilemanyvec(l,tetpil,z+1,2); return z;
! 1527:
! 1528: case t_QUAD:
! 1529: l=avma; p1=co8(y,lg(x)); tetpil=avma;
! 1530: return gerepile(l,tetpil,gdiv(x,p1));
! 1531:
! 1532: case t_INTMOD: case t_PADIC: err(gdiverf,tx,ty);
! 1533: }
! 1534:
! 1535: case t_INTMOD:
! 1536: switch(ty)
! 1537: {
! 1538: case t_INT: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
! 1539: (void)new_chunk(lgefint(p2)<<2);
! 1540: p1=mulii((GEN)x[2], mpinvmod(y,p2)); avma=(long)z;
! 1541: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 1542:
! 1543: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
! 1544: if (p1==p2 || egalii(p1,p2))
! 1545: { icopyifstack(p2,z[1]); }
! 1546: else
! 1547: { p2 = mppgcd(p1,p2); z[1] = (long)p2; }
! 1548: av=avma; (void)new_chunk(lgefint(x[1]) + (lgefint(p1) << 1));
! 1549: p1=mulii((GEN)x[2], mpinvmod((GEN)y[2],p2)); avma=av;
! 1550: z[2]=lmodii(p1,p2); return z;
! 1551:
! 1552: case t_FRAC: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
! 1553: (void)new_chunk(lgefint(p2)<<2);
! 1554: p1=mulii((GEN)y[2], mpinvmod((GEN)y[1],p2));
! 1555: p1=mulii(modii(p1,p2),(GEN)x[2]); avma=(long)z;
! 1556: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 1557:
! 1558: case t_FRACN:
! 1559: l=avma; p1=gred(y); tetpil=avma;
! 1560: return gerepile(l,tetpil,gdiv(x,p1));
! 1561:
! 1562: case t_COMPLEX: case t_QUAD:
! 1563: l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
! 1564: return gerepile(l,tetpil,gdiv(p2,p1));
! 1565:
! 1566: case t_PADIC:
! 1567: l=avma; p1=cgetg(3,t_INTMOD); p1[1]=x[1]; p1[2]=lgeti(lg(x[1]));
! 1568: gaffect(y,p1); tetpil=avma; return gerepile(l,tetpil,gdiv(x,p1));
! 1569:
! 1570: case t_REAL: err(gdiverf,tx,ty);
! 1571: }
! 1572:
! 1573: case t_FRAC: case t_FRACN:
! 1574: switch(ty)
! 1575: {
! 1576: case t_INT:
! 1577: z = cgetg(3, tx);
! 1578: if (tx == t_FRAC)
! 1579: {
! 1580: p1 = mppgcd(y,(GEN)x[1]);
! 1581: if (is_pm1(p1))
! 1582: {
! 1583: avma = (long)z; tetpil = 0;
! 1584: z[1] = licopy((GEN)x[1]);
! 1585: }
! 1586: else
! 1587: {
! 1588: y = divii(y,p1); tetpil = avma;
! 1589: z[1] = ldivii((GEN)x[1], p1);
! 1590: }
! 1591: }
! 1592: else
! 1593: {
! 1594: tetpil = 0;
! 1595: z[1] = licopy((GEN)x[1]);
! 1596: }
! 1597: z[2] = lmulii((GEN)x[2],y);
! 1598: fix_frac(z);
! 1599: if (tetpil) fix_frac_if_int_GC(z,tetpil);
! 1600: return z;
! 1601:
! 1602: case t_REAL:
! 1603: l=avma; p1=cgetg(lg(y),t_REAL); gaffect(x,p1);
! 1604: p2=divrr(p1,y); return gerepile(l,(long)p1,p2);
! 1605:
! 1606: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
! 1607: (void)new_chunk(lgefint(p2)<<2);
! 1608: p1=mulii((GEN)y[2],(GEN)x[2]);
! 1609: p1=mulii(mpinvmod(p1,p2), modii((GEN)x[1],p2)); avma=(long)z;
! 1610: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
! 1611:
! 1612: case t_FRAC: if (tx == t_FRACN) ty=t_FRACN;
! 1613: case t_FRACN:
! 1614: z = cgetg(3,ty);
! 1615: if (ty == t_FRAC)
! 1616: {
! 1617: GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
! 1618: GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
! 1619: p1 = mppgcd(x1, y1);
! 1620: if (!is_pm1(p1)) { x1 = divii(x1,p1); y1 = divii(y1,p1); }
! 1621: p1 = mppgcd(x2, y2);
! 1622: if (!is_pm1(p1)) { x2 = divii(x2,p1); y2 = divii(y2,p1); }
! 1623: tetpil = avma;
! 1624: z[2] = lmulii(x2,y1);
! 1625: z[1] = lmulii(x1,y2);
! 1626: fix_frac(z);
! 1627: fix_frac_if_int_GC(z,tetpil);
! 1628: }
! 1629: else
! 1630: {
! 1631: z[1]=lmulii((GEN)x[1],(GEN)y[2]);
! 1632: z[2]=lmulii((GEN)x[2],(GEN)y[1]);
! 1633: fix_frac(z);
! 1634: }
! 1635: return z;
! 1636:
! 1637: case t_COMPLEX: z=cgetg(3,t_COMPLEX);
! 1638: l=avma; p1=gnorm(y);
! 1639: p2=gmul(x,(GEN)y[1]);
! 1640: p3=gmul(x,(GEN)y[2]);
! 1641: if(!gcmp0(p3)) p3 = gneg_i(p3);
! 1642: tetpil=avma;
! 1643: z[1]=ldiv(p2,p1); z[2]=ldiv(p3,p1);
! 1644: gerepilemanyvec(l,tetpil,z+1,2); return z;
! 1645:
! 1646: case t_PADIC:
! 1647: if (!signe(x[1])) return gzero;
! 1648:
! 1649: l=avma; p1=cgetp(y); gaffect(x,p1);
! 1650: tetpil=avma; return gerepile(l,tetpil,gdiv(p1,y));
! 1651:
! 1652: case t_QUAD:
! 1653: l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
! 1654: return gerepile(l,tetpil,gdiv(p2,p1));
! 1655: }
! 1656:
! 1657: case t_COMPLEX:
! 1658: switch(ty)
! 1659: {
! 1660: case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FRACN:
! 1661: z=cgetg(3,t_COMPLEX);
! 1662: z[1]=ldiv((GEN)x[1],y);
! 1663: z[2]=ldiv((GEN)x[2],y); return z;
! 1664:
! 1665: case t_COMPLEX:
! 1666: l=avma; p1=gnorm(y); p2=gconj(y); p2=gmul(x,p2); tetpil=avma;
! 1667: return gerepile(l,tetpil, gdiv(p2,p1));
! 1668:
! 1669: case t_PADIC:
! 1670: if (krosg(-1,(GEN)y[2])== -1)
! 1671: {
! 1672: z=cgetg(3,t_COMPLEX);
! 1673: z[1]=ldiv((GEN)x[1],y);
! 1674: z[2]=ldiv((GEN)x[2],y); return z;
! 1675: }
! 1676: av=avma; p1=cvtop(x,(GEN)y[2],precp(y)); tetpil=avma;
! 1677: return gerepile(av,tetpil,gdiv(p1,y));
! 1678:
! 1679: case t_QUAD:
! 1680: lx=precision(x); if (!lx) err(gdiveri,tx,ty);
! 1681: l=avma; p1=co8(y,lx); tetpil=avma;
! 1682: return gerepile(l,tetpil,gdiv(x,p1));
! 1683: }
! 1684:
! 1685: case t_PADIC:
! 1686: switch(ty)
! 1687: {
! 1688: case t_INT: case t_FRAC: case t_FRACN:
! 1689: l=avma;
! 1690: if (signe(x[4])) { p1=cgetp(x); gaffect(y,p1); }
! 1691: else p1=cvtop(y,(GEN)x[2],(valp(x)>0)?valp(x):1);
! 1692: tetpil=avma; return gerepile(l,tetpil,gdiv(x,p1));
! 1693:
! 1694: case t_INTMOD:
! 1695: l=avma; p1=cgetg(3,t_INTMOD);
! 1696: p1[1]=y[1]; p1[2]=lgeti(lg(y[1]));
! 1697: gaffect(x,p1); tetpil=avma;
! 1698: return gerepile(l,tetpil,gdiv(p1,y));
! 1699:
! 1700: case t_PADIC:
! 1701: if (!egalii((GEN)x[2],(GEN)y[2])) err(gdiveri,tx,ty);
! 1702: if (!signe(x[4]))
! 1703: {
! 1704: z=gcopy(x); setvalp(z,valp(x)-valp(y));
! 1705: return z;
! 1706: }
! 1707:
! 1708: p1=(precp(x)>precp(y)) ? y : x;
! 1709: z=cgetp(p1); l=avma;
! 1710: setvalp(z,valp(x)-valp(y));
! 1711: p2=mpinvmod((GEN)y[4],(GEN)p1[3]);
! 1712: modiiz(mulii((GEN)x[4],p2),(GEN)p1[3],(GEN)z[4]);
! 1713: avma=l; return z;
! 1714:
! 1715: case t_COMPLEX: case t_QUAD:
! 1716: l=avma; p1=gmul(x,gconj(y)); p2=gnorm(y); tetpil=avma;
! 1717: return gerepile(l,tetpil,gdiv(p1,p2));
! 1718:
! 1719: case t_REAL:
! 1720: err(talker,"forbidden division p-adic/R");
! 1721: }
! 1722:
! 1723: case t_QUAD:
! 1724: switch (ty)
! 1725: {
! 1726: case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN:
! 1727: z=cgetg(4,t_QUAD);
! 1728: copyifstack(x[1], z[1]);
! 1729: for (i=2; i<4; i++) z[i]=ldiv((GEN)x[i],y);
! 1730: return z;
! 1731:
! 1732: case t_REAL:
! 1733: l=avma; p1=co8(x,lg(y)); tetpil=avma;
! 1734: return gerepile(l,tetpil,gdiv(p1,y));
! 1735:
! 1736: case t_PADIC:
! 1737: l=avma; p1=cvtop(x,(GEN)y[2],precp(y));
! 1738: tetpil=avma; return gerepile(l,tetpil,gdiv(p1,y));
! 1739:
! 1740: case t_COMPLEX:
! 1741: ly=precision(y); if (!ly) err(gdiveri,tx,ty);
! 1742: l=avma; p1=co8(x,ly); tetpil=avma;
! 1743: return gerepile(l,tetpil,gdiv(p1,y));
! 1744:
! 1745: case t_QUAD:
! 1746: k=x[1]; l=y[1];
! 1747: if (!gegal((GEN)k,(GEN)l)) err(gdiveri,tx,ty);
! 1748: l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
! 1749: return gerepile(l,tetpil,gdiv(p2,p1));
! 1750: }
! 1751: }
! 1752: err(bugparier,"division");
! 1753: }
! 1754:
! 1755: vx=gvar(x); vy=gvar(y);
! 1756: if (ty==t_POLMOD && (tx==t_POLMOD || vy<vx))
! 1757: {
! 1758: z=cgetg(3,t_POLMOD);
! 1759: if (tx==t_POLMOD)
! 1760: {
! 1761: k=x[1]; l=y[1];
! 1762: if (gegal((GEN)k,(GEN)l))
! 1763: {
! 1764: copyifstack(k, z[1]); av=avma;
! 1765: p1 = ginvmod((GEN)y[2],(GEN)z[1]);
! 1766: p2 = gmul((GEN)x[2],p1);
! 1767: }
! 1768: else
! 1769: {
! 1770: vx=varn(x[1]); vy=varn(y[1]);
! 1771: if (vx==vy)
! 1772: {
! 1773: z[1]=lgcd((GEN)k,(GEN)l); av=avma;
! 1774: p1=ginvmod((GEN)y[2],(GEN)z[1]);
! 1775: p2=gmul((GEN)x[2],p1);
! 1776: }
! 1777: else
! 1778: {
! 1779: if (vx<vy)
! 1780: { copyifstack(k,z[1]); av=avma; p2=gdiv((GEN)x[2],y); }
! 1781: else
! 1782: {
! 1783: copyifstack(l,z[1]); av=avma;
! 1784: p1 = ginvmod((GEN)y[2],(GEN)z[1]);
! 1785: p2 = gmul(x, p1);
! 1786: }
! 1787: }
! 1788: }
! 1789: p2 = gmod(p2,(GEN)z[1]);
! 1790: }
! 1791: else
! 1792: {
! 1793: copyifstack(y[1],z[1]); av=avma;
! 1794: p1 = ginvmod((GEN)y[2],(GEN)y[1]);
! 1795: p2 = gmul(x,p1);
! 1796: }
! 1797: z[2]=lpileupto(av, p2); return z;
! 1798: }
! 1799: if (tx == t_POLMOD && vx<vy)
! 1800: {
! 1801: z=cgetg(3,t_POLMOD);
! 1802: copyifstack(x[1],z[1]);
! 1803: z[2]=ldiv((GEN)x[2],y); return z;
! 1804: }
! 1805: if (vx == vy)
! 1806: {
! 1807: av = avma;
! 1808: if (tx == t_POLMOD)
! 1809: return gerepileupto(av, gdiv(x, to_polmod(y,(GEN)x[1])));
! 1810: if (ty == t_POLMOD)
! 1811: return gerepileupto(av, gdiv(to_polmod(x,(GEN)y[1]), y));
! 1812: }
! 1813: /* now x and y are not both is_scalar_t */
! 1814:
! 1815: lx = lg(x);
! 1816: if ((vx<vy && (!is_matvec_t(tx) || !is_matvec_t(ty)))
! 1817: || (vx==vy && is_scalar_t(ty)) || (is_matvec_t(tx) && !is_matvec_t(ty)))
! 1818: {
! 1819: if (tx == t_RFRAC) return divrfracscal(x,y);
! 1820: z = cgetg(lx,tx);
! 1821: if (tx == t_RFRACN)
! 1822: {
! 1823: z[2]=lmul((GEN)x[2],y);
! 1824: z[1]=lcopy((GEN)x[1]); return z;
! 1825: }
! 1826: switch(tx)
! 1827: {
! 1828: case t_POL: lx = lgef(x);
! 1829: case t_SER: z[1] = x[1];
! 1830: case t_VEC: case t_COL: case t_MAT:
! 1831: if (ty == t_POLMOD || ty == t_INTMOD)
! 1832: {
! 1833: if (!gcmp1(y)) y = ginv(y); /* garbage, left alone */
! 1834: for (i=lontyp[tx]; i<lx; i++) z[i]=lmul((GEN)x[i],y);
! 1835: return z;
! 1836: }
! 1837: else
! 1838: for (i=lontyp[tx]; i<lx; i++) z[i]=ldiv((GEN)x[i],y);
! 1839: return z;
! 1840: }
! 1841: err(typeer,"division");
! 1842: }
! 1843:
! 1844: ly=lg(y);
! 1845: if (vy<vx || (vy==vx && is_scalar_t(tx)))
! 1846: {
! 1847: switch(ty)
! 1848: {
! 1849: case t_POL:
! 1850: if (lgef(y)==3) return gdiv(x,(GEN)y[2]);
! 1851: if (isexactzero(x)) return zeropol(vy);
! 1852: av=avma; z=cgetg(3,t_RFRAC); z[1]=(long)x; z[2]=(long)y;
! 1853: return gerepileupto(av,gred_rfrac(z));
! 1854:
! 1855: case t_SER:
! 1856: if (gcmp0(x))
! 1857: {
! 1858: l=avma; p1=ginv(y); tetpil=avma; /* a ameliorer !!!! */
! 1859: return gerepile(l,tetpil,gmul(x,p1));
! 1860: }
! 1861: p1 = (GEN)gpmalloc(ly*sizeof(long));
! 1862: p1[0] = evaltyp(t_SER) | evallg(ly);
! 1863: p1[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
! 1864: p1[2] = (long)x; for (i=3; i<ly; i++) p1[i]=zero;
! 1865: y = gdiv(p1,y); free(p1); return y;
! 1866:
! 1867: case t_RFRAC: return divscalrfrac(x,y);
! 1868: case t_RFRACN: z=cgetg(ly,t_RFRACN);
! 1869: z[1]=lmul(x,(GEN)y[2]);
! 1870: z[2]=lcopy((GEN)y[1]); return z;
! 1871:
! 1872: case t_MAT:
! 1873: if (ly==1 || lg(y[1])!=ly) err(gdiveri,tx,ty);
! 1874: l=avma; p1=invmat(y); tetpil=avma;
! 1875: return gerepile(l,tetpil,gmul(x,p1));
! 1876:
! 1877: case t_VEC: case t_COL: err(gdiverf,tx,ty);
! 1878: }
! 1879: err(typeer,"division");
! 1880: }
! 1881:
! 1882: /* ici vx=vy et tx>=10 et ty>=10*/
! 1883: switch(tx)
! 1884: {
! 1885: case t_POL:
! 1886: switch(ty)
! 1887: {
! 1888: case t_POL:
! 1889: if (lgef(y)==3) return gdiv(x,(GEN)y[2]);
! 1890: if (isexactzero(x)) return zeropol(vy);
! 1891: av=avma; z=cgetg(3,t_RFRAC); z[1]=(long)x; z[2]=(long)y;
! 1892: return gerepileupto(av,gred_rfrac(z));
! 1893:
! 1894: case t_SER:
! 1895: if (gcmp0(x)) return zeropol(vx);
! 1896: p1=greffe(x,ly,0); p2=gdiv(p1,y);
! 1897: free(p1); return p2;
! 1898:
! 1899: case t_RFRAC: return divscalrfrac(x,y);
! 1900: case t_RFRACN: z=cgetg(ly,t_RFRACN);
! 1901: z[1]=lmul(x,(GEN)y[2]);
! 1902: z[2]=lcopy((GEN)y[1]); return z;
! 1903:
! 1904: case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
! 1905: default: err(typeer,"division");
! 1906: }
! 1907:
! 1908: case t_SER:
! 1909: switch(ty)
! 1910: {
! 1911: case t_POL:
! 1912: p1=greffe(y,lx,0); p2=gdiv(x,p1);
! 1913: free(p1); return p2;
! 1914:
! 1915: case t_SER:
! 1916: {
! 1917: GEN y_lead;
! 1918:
! 1919: l = valp(x) - valp(y);
! 1920: if (gcmp0(x)) return zeroser(vx,l);
! 1921: y_lead = (GEN)y[2];
! 1922: if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */
! 1923: {
! 1924: err(warner,"normalizing a series with 0 leading term");
! 1925: for (i=3,y++; i<ly; i++,y++)
! 1926: {
! 1927: y_lead = (GEN)y[2]; ly--; l--;
! 1928: if (!gcmp0(y_lead)) break;
! 1929: }
! 1930: if (i==ly) err(gdiver2);
! 1931: }
! 1932: if (ly < lx) lx = ly;
! 1933: p2 = (GEN)gpmalloc(lx*sizeof(long));
! 1934: for (i=3; i<lx; i++)
! 1935: {
! 1936: p1 = (GEN)y[i];
! 1937: if (isexactzero(p1)) p2[i] = 0;
! 1938: else
! 1939: {
! 1940: av = avma; p2[i] = lclone(gneg_i(p1));
! 1941: avma = av;
! 1942: }
! 1943: }
! 1944: z = cgetg(lx,t_SER);
! 1945: z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
! 1946: z[2] = ldiv((GEN)x[2], y_lead);
! 1947: for (i=3; i<lx; i++)
! 1948: {
! 1949: av=avma; p1 = (GEN)x[i];
! 1950: for (j=2; j<i; j++)
! 1951: {
! 1952: l = i-j+2;
! 1953: if (p2[l])
! 1954: p1 = gadd(p1, gmul((GEN)z[j], (GEN)p2[l]));
! 1955: }
! 1956: tetpil=avma; z[i]=lpile(av,tetpil, gdiv(p1,y_lead));
! 1957: }
! 1958: for (i=3; i<lx; i++)
! 1959: if (p2[i]) gunclone((GEN)p2[i]);
! 1960: free(p2); return z;
! 1961: }
! 1962:
! 1963: case t_RFRAC: case t_RFRACN:
! 1964: l=avma; p2=gmul(x,(GEN)y[2]); tetpil=avma;
! 1965: return gerepile(l,tetpil,gdiv(p2,(GEN)y[1]));
! 1966:
! 1967: case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
! 1968: default: err(typeer,"division");
! 1969: }
! 1970:
! 1971: case t_RFRAC: case t_RFRACN:
! 1972: switch(ty)
! 1973: {
! 1974: case t_POL:
! 1975: if (tx==t_RFRAC) return divrfracscal(x,y);
! 1976: z=cgetg(3,t_RFRACN);
! 1977: z[2]=lmul((GEN)x[2],y);
! 1978: z[1]=lcopy((GEN)x[1]); return z;
! 1979:
! 1980: case t_SER:
! 1981: l=avma; p2=gmul((GEN)x[2],y); tetpil=avma;
! 1982: return gerepile(l,tetpil, gdiv((GEN)x[1],p2));
! 1983:
! 1984: case t_RFRAC: case t_RFRACN:
! 1985: if (tx == t_RFRACN) ty=t_RFRACN;
! 1986: if (ty != t_RFRACN) return divrfrac(x,y);
! 1987: z=cgetg(3,t_RFRACN);
! 1988: z[1]=lmul((GEN)x[1],(GEN)y[2]);
! 1989: z[2]=lmul((GEN)x[2],(GEN)y[1]); return z;
! 1990:
! 1991: case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
! 1992: default: err(typeer,"division");
! 1993: }
! 1994:
! 1995: case t_VEC: case t_COL: case t_MAT:
! 1996: if (!is_matvec_t(ty))
! 1997: {
! 1998: z=cgetg(lx,tx);
! 1999: for (i=1; i<lx; i++) z[i]=ldiv((GEN)x[i],y);
! 2000: return z;
! 2001: }
! 2002: if (ty!=t_MAT || ly==1 || lg(y[1])!=ly) err(gdiveri,tx,ty);
! 2003: l=avma; p1=invmat(y); tetpil=avma;
! 2004: return gerepile(l,tetpil,gmul(x,p1));
! 2005: }
! 2006: if (tx==ty)
! 2007: {
! 2008: l=signe(y[2]); setsigne(y[2],-l);
! 2009: switch(tx)
! 2010: {
! 2011: case t_QFI: z = compimag(x,y);
! 2012: setsigne(y[2],l); return z;
! 2013: case t_QFR:
! 2014: k=signe(y[4]); setsigne(y[4],-k); z=compreal(x,y);
! 2015: setsigne(y[2],l); setsigne(y[4],k); return z;
! 2016: }
! 2017: }
! 2018: err(typeer,"division");
! 2019: return NULL; /* not reached */
! 2020: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen1.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