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