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