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