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