Annotation of OpenXM_contrib/pari/src/basemath/gen3.c, Revision 1.1.1.1
1.1 maekawa 1: /********************************************************************/
2: /** **/
3: /** GENERIC OPERATIONS **/
4: /** (third part) **/
5: /** **/
6: /********************************************************************/
7: /* $Id: gen3.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
8: #include "pari.h"
9:
10: /********************************************************************/
11: /** **/
12: /** PRINCIPAL VARIABLE NUMBER **/
13: /** **/
14: /********************************************************************/
15:
16: int
17: gvar(GEN x)
18: {
19: long tx=typ(x),i,v,w;
20:
21: if (is_polser_t(tx)) return varn(x);
22: if (tx==t_POLMOD) return varn(x[1]);
23: if (is_const_t(tx) || is_qf_t(tx) || tx == t_STR || tx == t_LIST)
24: return BIGINT;
25: v=BIGINT;
26: for (i=1; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
27: return v;
28: }
29:
30: int
31: gvar2(GEN x)
32: {
33: long tx=typ(x),i,v,w;
34:
35: if (is_const_t(tx) || is_qf_t(tx)) return BIGINT;
36: v = BIGINT;
37: switch(tx)
38: {
39: case t_POLMOD:
40: v=gvar2((GEN)x[1]);
41: w=gvar2((GEN)x[2]); if (w<v) v=w;
42: return v;
43:
44: case t_SER:
45: for (i=2; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
46: return v;
47:
48: case t_POL:
49: for (i=2; i<lgef(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
50: return v;
51: }
52: for (i=1; i<lg(x); i++) { w=gvar2((GEN)x[i]); if (w<v) v=w; }
53: return v;
54: }
55:
56: GEN
57: gpolvar(GEN x)
58: {
59: if (typ(x)==t_PADIC) x = (GEN)x[2];
60: else
61: {
62: long v=gvar(x);
63: if (v==BIGINT) err(typeer,"polvar");
64: x = polx[v];
65: }
66: return gcopy(x);
67: }
68:
69: /*******************************************************************/
70: /* */
71: /* PRECISION OF SCALAR OBJECTS */
72: /* */
73: /*******************************************************************/
74:
75: long
76: precision(GEN x)
77: {
78: long tx=typ(x),k,l;
79:
80: if (tx==t_REAL)
81: {
82: k=2-(expo(x)>>TWOPOTBITS_IN_LONG);
83: l=lg(x); if (l>k) k=l;
84: return k;
85: }
86: if (tx==t_COMPLEX)
87: {
88: k=precision((GEN)x[1]);
89: l=precision((GEN)x[2]); if (l && l<k) k=l;
90: return k;
91: }
92: return 0;
93: }
94:
95: long
96: gprecision(GEN x)
97: {
98: long tx=typ(x),lx=lg(x),i,k,l;
99:
100: if (is_scalar_t(tx)) return precision(x);
101: switch(tx)
102: {
103: case t_POL: lx=lgef(x);
104: case t_VEC: case t_COL: case t_MAT:
105: k=VERYBIGINT;
106: for (i=lontyp[tx]; i<lx; i++)
107: {
108: l = gprecision((GEN)x[i]);
109: if (l && l<k) k = l;
110: }
111: return (k==VERYBIGINT)? 0: k;
112:
113: case t_RFRAC: case t_RFRACN:
114: {
115: k=gprecision((GEN)x[1]);
116: l=gprecision((GEN)x[2]); if (l && l<k) k=l;
117: return k;
118: }
119: case t_QFR:
120: return gprecision((GEN)x[4]);
121: }
122: return 0;
123: }
124:
125: GEN
126: ggprecision(GEN x)
127: {
128: long a=gprecision(x);
129: return stoi(a ? (long) ((a-2)*pariK): VERYBIGINT);
130: }
131:
132: GEN
133: precision0(GEN x, long n)
134: {
135: if (n) return gprec(x,n);
136: return ggprecision(x);
137: }
138:
139: /* attention: precision p-adique absolue */
140: long
141: padicprec(GEN x, GEN p)
142: {
143: long i,s,t,lx=lg(x),tx=typ(x);
144:
145: switch(tx)
146: {
147: case t_INT: case t_FRAC: case t_FRACN:
148: return VERYBIGINT;
149:
150: case t_INTMOD:
151: return ggval((GEN)x[1],p);
152:
153: case t_PADIC:
154: if (!gegal((GEN)x[2],p))
155: err(talker,"not the same prime in padicprec");
156: return precp(x)+valp(x);
157:
158: case t_POL:
159: lx=lgef(x);
160:
161: case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_SER: case t_RFRAC:
162: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
163: for (s=VERYBIGINT, i=lontyp[tx]; i<lx; i++)
164: {
165: t = padicprec((GEN)x[i],p); if (t<s) s = t;
166: }
167: return s;
168: }
169: err(typeer,"padicprec");
170: return 0; /* not reached */
171: }
172:
173: /* Degre de x par rapport a la variable v si v>=0, par rapport a la variable
174: * principale si v<0. On suppose x fraction rationnelle ou polynome.
175: * Convention deg(0)=-1.
176: */
177: long
178: poldegree(GEN x, long v)
179: {
180: long tx=typ(x), av, w, d;
181:
182: if (is_scalar_t(tx)) return gcmp0(x)? -1: 0;
183: switch(tx)
184: {
185: case t_POL:
186: w = varn(x);
187: if (v < 0 || v == w) return lgef(x)-3;
188: if (v < w) return signe(x)? 0: -1;
189: av = avma; x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
190: if (gvar(x)) { d = gcmp0(x)? -1: 0; } else d = lgef(x)-3;
191: avma = av; return d;
192:
193: case t_RFRAC: case t_RFRACN:
194: if (gcmp0((GEN)x[1])) return -1;
195: return poldegree((GEN)x[1],v) - poldegree((GEN)x[2],v);
196: }
197: err(typeer,"degree");
198: return 0; /* not reached */
199: }
200:
201: long
202: degree(GEN x)
203: {
204: return poldegree(x,-1);
205: }
206:
207: /* si v<0, par rapport a la variable principale, sinon par rapport a v.
208: * On suppose que x est un polynome ou une serie.
209: */
210: GEN
211: pollead(GEN x, long v)
212: {
213: long l,tx = typ(x),av,tetpil,w;
214: GEN xinit;
215:
216: if (is_scalar_t(tx)) return gcopy(x);
217: w = varn(x);
218: switch(tx)
219: {
220: case t_POL:
221: if (v < 0 || v == w)
222: {
223: l=lgef(x);
224: return (l==2)? gzero: gcopy((GEN)x[l-1]);
225: }
226: if (v < w) return gcopy(x);
227: av = avma; xinit = x;
228: x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
229: if (gvar(x)) { avma = av; return gcopy(xinit); }
230: l = lgef(x); if (l == 2) { avma = av; return gzero; }
231: tetpil = avma; x = gsubst((GEN)x[l-1],MAXVARN,polx[w]);
232: return gerepile(av,tetpil,x);
233:
234: case t_SER:
235: if (v < 0 || v == w) return (!signe(x))? gzero: gcopy((GEN)x[2]);
236: if (v < w) return gcopy(x);
237: av = avma; xinit = x;
238: x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
239: if (gvar(x)) { avma = av; return gcopy(xinit);}
240: if (!signe(x)) { avma = av; return gzero;}
241: tetpil = avma; x = gsubst((GEN)x[2],MAXVARN,polx[w]);
242: return gerepile(av,tetpil,x);
243: }
244: err(typeer,"pollead");
245: return NULL; /* not reached */
246: }
247:
248: /* returns 1 if there's a real component in the structure, 0 otherwise */
249: int
250: isinexactreal(GEN x)
251: {
252: long tx=typ(x),lx,i;
253:
254: if (is_scalar_t(tx))
255: {
256: switch(tx)
257: {
258: case t_REAL:
259: return 1;
260:
261: case t_COMPLEX: case t_QUAD:
262: return (typ(x[1])==t_REAL || typ(x[2])==t_REAL);
263: }
264: return 0;
265: }
266: switch(tx)
267: {
268: case t_QFR: case t_QFI:
269: return 0;
270:
271: case t_RFRAC: case t_RFRACN:
272: return isinexactreal((GEN)x[1]) || isinexactreal((GEN)x[2]);
273: }
274: lx = (tx==t_POL)? lgef(x): lg(x);
275: for (i=lontyp[tx]; i<lx; i++)
276: if (isinexactreal((GEN)x[i])) return 1;
277: return 0;
278: }
279:
280: int
281: isexactzero(GEN g)
282: {
283: long i;
284: switch (typ(g))
285: {
286: case t_INT:
287: return !signe(g);
288: case t_REAL: case t_PADIC: case t_SER:
289: return 0;
290: case t_INTMOD: case t_POLMOD:
291: return isexactzero((GEN)g[2]);
292: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
293: return isexactzero((GEN)g[1]);
294: case t_COMPLEX:
295: return isexactzero((GEN)g[1]) && isexactzero((GEN)g[2]);
296: case t_QUAD:
297: return isexactzero((GEN)g[2]) && isexactzero((GEN)g[3]);
298:
299: case t_POL:
300: for (i=lgef(g)-1; i>1; i--)
301: if (!isexactzero((GEN)g[i])) return 0;
302: return 1;
303: case t_VEC: case t_COL: case t_MAT:
304: for (i=lg(g)-1; i; i--)
305: if (!isexactzero((GEN)g[i])) return 0;
306: return 1;
307: }
308: return 0;
309: }
310:
311: int
312: iscomplex(GEN x)
313: {
314: switch(typ(x))
315: {
316: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
317: return 0;
318: case t_COMPLEX:
319: return !gcmp0((GEN)x[2]);
320: case t_QUAD:
321: return signe(mael(x,1,2)) > 0;
322: }
323: err(typeer,"iscomplex");
324: return 0; /* not reached */
325: }
326:
327: int
328: ismonome(GEN x)
329: {
330: long i;
331: if (typ(x)!=t_POL || !signe(x)) return 0;
332: for (i=lgef(x)-2; i>1; i--)
333: if (!isexactzero((GEN)x[i])) return 0;
334: return 1;
335: }
336:
337: /********************************************************************/
338: /** **/
339: /** MULTIPLICATION SIMPLE **/
340: /** **/
341: /********************************************************************/
342: #define fix_frac(z) if (signe(z[2])<0)\
343: {\
344: setsigne(z[1],-signe(z[1]));\
345: setsigne(z[2],1);\
346: }
347:
348: /* assume z[1] was created last */
349: #define fix_frac_if_int(z) if (is_pm1(z[2]))\
350: z = gerepileupto((long)(z+3), (GEN)z[1]);
351:
352: GEN
353: gmulsg(long s, GEN y)
354: {
355: long ty=typ(y),ly=lg(y),i,av,tetpil;
356: GEN z,p1,p2;
357:
358: switch(ty)
359: {
360: case t_INT:
361: return mulsi(s,y);
362:
363: case t_REAL:
364: return mulsr(s,y);
365:
366: case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
367: (void)new_chunk(lgefint(p2)<<2);
368: p1=mulsi(s,(GEN)y[2]); avma=(long)z;
369: z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
370:
371: case t_FRAC:
372: if (!s) return gzero;
373: z = cgetg(3,t_FRAC);
374: i = cgcd(s, smodis((GEN)y[2], s));
375: if (i == 1)
376: {
377: z[2] = licopy((GEN)y[2]);
378: z[1] = lmulis((GEN)y[1], s);
379: }
380: else
381: {
382: z[2] = ldivis((GEN)y[2], i);
383: z[1] = lmulis((GEN)y[1], s/i);
384: fix_frac_if_int(z);
385: }
386: return z;
387:
388: case t_FRACN: z=cgetg(3,ty);
389: z[1]=lmulsi(s,(GEN)y[1]);
390: z[2]=licopy((GEN)y[2]);
391: return z;
392:
393: case t_COMPLEX: z=cgetg(ly,ty);
394: z[1]=lmulsg(s,(GEN)y[1]);
395: z[2]=lmulsg(s,(GEN)y[2]); return z;
396:
397: case t_PADIC:
398: if (!s) return gzero;
399: av=avma; p1=cgetp(y); gaffsg(s,p1); tetpil=avma;
400: return gerepile(av,tetpil,gmul(p1,y));
401:
402: case t_QUAD: z=cgetg(ly,ty);
403: copyifstack(y[1],z[1]);
404: z[2]=lmulsg(s,(GEN)y[2]);
405: z[3]=lmulsg(s,(GEN)y[3]); return z;
406:
407: case t_POLMOD: z=cgetg(ly,ty);
408: z[2]=lmulsg(s,(GEN)y[2]);
409: copyifstack(y[1],z[1]); return z;
410:
411: case t_POL:
412: if (!s || !signe(y)) return zeropol(varn(y));
413: ly=lgef(y); z=cgetg(ly,t_POL); z[1]=y[1];
414: for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
415: return normalizepol_i(z, ly);
416:
417: case t_SER:
418: if (!s) return zeropol(varn(y));
419: if (gcmp0(y)) return gcopy(y);
420: z=cgetg(ly,ty);
421: for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
422: z[1]=y[1]; return normalize(z);
423:
424: case t_RFRAC:
425: if (!s) return zeropol(gvar(y));
426: z = cgetg(3, t_RFRAC);
427: i = ggcd(stoi(s),(GEN)y[2])[2];
428: avma = (long)z;
429: if (i == 1)
430: {
431: z[1]=lmulgs((GEN)y[1], s);
432: z[2]= lcopy((GEN)y[2]);
433: }
434: else
435: {
436: z[1] = lmulgs((GEN)y[1], s/i);
437: z[2] = ldivgs((GEN)y[2], i);
438: }
439: return z;
440:
441: case t_RFRACN:
442: if (!s) return zeropol(gvar(y));
443: z=cgetg(ly,t_RFRACN);
444: z[1]=lmulsg(s,(GEN)y[1]);
445: z[2]=lcopy((GEN)y[2]); return z;
446:
447: case t_VEC: case t_COL: case t_MAT:
448: z=cgetg(ly,ty);
449: for (i=1; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
450: return z;
451: }
452: err(typeer,"gmulsg");
453: return NULL; /* not reached */
454: }
455:
456: /********************************************************************/
457: /** **/
458: /** DIVISION SIMPLE **/
459: /** **/
460: /********************************************************************/
461:
462: GEN
463: gdivgs(GEN x, long s)
464: {
465: static long court[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
466: long tx=typ(x),lx=lg(x),i,av;
467: GEN z,p1;
468:
469: if (!s) err(gdiver2);
470: switch(tx)
471: {
472: case t_INT:
473: av=avma; z=dvmdis(x,s,&p1);
474: if (p1==gzero) return z;
475:
476: i = cgcd(s, smodis(x,s));
477: avma=av; z=cgetg(3,t_FRAC);
478: if (i == 1)
479: x = icopy(x);
480: else
481: {
482: s /= i;
483: x = divis(x, i);
484: }
485: z[1]=(long)x; z[2]=lstoi(s);
486: fix_frac(z); return z;
487:
488: case t_REAL:
489: return divrs(x,s);
490:
491: case t_FRAC: z=cgetg(3,tx);
492: i = cgcd(s, smodis((GEN)x[1], s));
493: if (i == 1)
494: {
495: z[2] = lmulsi(s, (GEN)x[2]);
496: z[1] = licopy((GEN)x[1]);
497: }
498: else
499: {
500: z[2] = lmulsi(s/i, (GEN)x[2]);
501: z[1] = ldivis((GEN)x[1], i);
502: }
503: fix_frac(z);
504: fix_frac_if_int(z); return z;
505:
506: case t_FRACN: z = cgetg(3,tx);
507: z[1]=licopy((GEN)x[1]);
508: z[2]=lmulsi(s,(GEN)x[2]);
509: fix_frac(z); return z;
510:
511: case t_COMPLEX: z=cgetg(lx,tx);
512: z[1]=ldivgs((GEN)x[1],s);
513: z[2]=ldivgs((GEN)x[2],s);
514: return z;
515:
516: case t_QUAD: z=cgetg(lx,tx);
517: copyifstack(x[1],z[1]);
518: for (i=2; i<4; i++) z[i]=ldivgs((GEN)x[i],s);
519: return z;
520:
521: case t_POLMOD: z=cgetg(lx,tx);
522: copyifstack(x[1],z[1]);
523: z[2]=ldivgs((GEN)x[2],s);
524: return z;
525:
526: case t_POL: lx=lgef(x); z=cgetg(lx,tx);
527: for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
528: z[1]=x[1]; return z;
529:
530: case t_SER: z=cgetg(lx,tx);
531: for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
532: z[1]=x[1]; return z;
533:
534: case t_RFRAC:
535: z = cgetg(3, t_RFRAC);
536: i = ggcd(stoi(s),(GEN)x[1])[2];
537: avma = (long)z;
538: if (i == 1)
539: {
540: z[2]=lmulsg(s,(GEN)x[2]);
541: z[1]=lcopy((GEN)x[1]); return z;
542: }
543: z[1] = ldivgs((GEN)x[1], i);
544: z[2] = lmulgs((GEN)x[2], s/i); return z;
545:
546: case t_RFRACN: z=cgetg(3,t_RFRACN);
547: z[2]=lmulsg(s,(GEN)x[2]);
548: z[1]=lcopy((GEN)x[1]); return z;
549:
550: case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
551: for (i=1; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
552: return z;
553: }
554: affsi(s,court); return gdiv(x,court);
555: }
556:
557: /*******************************************************************/
558: /* */
559: /* MODULO GENERAL */
560: /* */
561: /*******************************************************************/
562:
563: GEN
564: gmod(GEN x, GEN y)
565: {
566: long av,tetpil,i, tx=typ(x);
567: GEN z,p1;
568:
569: if (is_matvec_t(tx))
570: {
571: i=lg(x); z=cgetg(i,tx);
572: for (i--; i; i--) z[i]=lmod((GEN)x[i],y);
573: return z;
574: }
575: switch(typ(y))
576: {
577: case t_INT:
578: switch(tx)
579: {
580: case t_INT:
581: return modii(x,y);
582:
583: case t_INTMOD: z=cgetg(3,tx);
584: z[1]=lmppgcd((GEN)x[1],y);
585: z[2]=lmodii((GEN)x[2],(GEN)z[1]); return z;
586:
587: case t_FRAC: case t_FRACN:
588: av=avma; if (tx==t_FRACN) x=gred(x);
589: p1=mulii((GEN)x[1],mpinvmod((GEN)x[2],y));
590: tetpil=avma; return gerepile(av,tetpil,modii(p1,y));
591:
592: case t_QUAD: z=cgetg(4,tx);
593: copyifstack(x[1],z[1]);
594: z[2]=lmod((GEN)x[2],y);
595: z[3]=lmod((GEN)x[3],y); return z;
596:
597: case t_PADIC:
598: {
599: long p1[] = {evaltyp(t_INTMOD)|m_evallg(3),0,0};
600: p1[1] = (long)y;
601: p1[2] = lgeti(lgefint(y));
602: gaffect(x,p1); return (GEN)p1[2];
603: }
604: case t_POLMOD: case t_POL:
605: return gzero;
606:
607: default: err(gmoder1);
608: }
609:
610: case t_REAL: case t_FRAC: case t_FRACN:
611: switch(tx)
612: {
613: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
614: av=avma; p1 = gfloor(gdiv(x,y)); p1 = gneg_i(gmul(p1,y));
615: tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
616:
617: case t_POLMOD: case t_POL:
618: return gzero;
619:
620: default: err(gmoder1);
621: }
622:
623: case t_POL:
624: if (is_scalar_t(tx))
625: {
626: if (tx!=t_POLMOD || varn(x[1]) > varn(y))
627: return (lgef(y) < 4)? gzero: gcopy(x);
628: if (varn(x[1])!=varn(y)) return gzero;
629: z=cgetg(3,t_POLMOD);
630: z[1]=lgcd((GEN)x[1],y);
631: z[2]=lres((GEN)x[2],(GEN)z[1]); return z;
632: }
633: switch(tx)
634: {
635: case t_POL:
636: return gres(x,y);
637:
638: case t_RFRAC: case t_RFRACN:
639: av=avma; if (tx==t_RFRACN) x=gred_rfrac(x);
640: p1=gmul((GEN)x[1],ginvmod((GEN)x[2],y)); tetpil=avma;
641: return gerepile(av,tetpil,gres(p1,y));
642:
643: default: err(talker,"type mod polynomial forbidden in gmod");
644: }
645: }
646: err(talker,"modulus type forbidden in gmod");
647: return NULL; /* not reached */
648: }
649:
650: GEN
651: gmodulsg(long x, GEN y)
652: {
653: GEN z;
654:
655: switch(typ(y))
656: {
657: case t_INT: z=cgetg(3,t_INTMOD);
658: z[1]=labsi(y); z[2]=lmodsi(x,y); return z;
659:
660: case t_POL: z=cgetg(3,t_POLMOD);
661: z[1]=lcopy(y); z[2]=lstoi(x); return z;
662: }
663: err(gmoder1); return NULL; /* not reached */
664: }
665:
666: GEN
667: gmodulss(long x, long y)
668: {
669: GEN z=cgetg(3,t_INTMOD);
670:
671: y=labs(y); z[1]=lstoi(y); z[2]=lstoi(x % y); return z;
672: }
673:
674: GEN
675: gmodulo(GEN x,GEN y)
676: {
677: long tx=typ(x),l,i;
678: GEN z;
679:
680: if (is_matvec_t(tx))
681: {
682: l=lg(x); z=cgetg(l,tx);
683: for (i=1; i<l; i++) z[i] = (long) gmodulo((GEN)x[i],y);
684: return z;
685: }
686: switch(typ(y))
687: {
688: case t_INT:
689: if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
690: z=cgetg(3,t_INTMOD);
691: if (!signe(y)) err(talker,"zero modulus in gmodulo");
692: y = gclone(y); setsigne(y,1);
693: z[1]=(long)y;
694: z[2]=lmod(x,y); return z;
695:
696: case t_POL: z=cgetg(3,t_POLMOD);
697: z[1] = lclone(y);
698: if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
699: if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
700: }
701: err(gmoder1); return NULL; /* not reached */
702: }
703:
704:
705: GEN
706: gmodulcp(GEN x,GEN y)
707: {
708: long tx=typ(x),l,i;
709: GEN z;
710:
711: if (is_matvec_t(tx))
712: {
713: l=lg(x); z=cgetg(l,tx);
714: for (i=1; i<l; i++) z[i] = lmodulcp((GEN)x[i],y);
715: return z;
716: }
717: switch(typ(y))
718: {
719: case t_INT:
720: if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
721: z=cgetg(3,t_INTMOD);
722: z[1]=labsi(y);
723: z[2]=lmod(x,y); return z;
724:
725: case t_POL: z=cgetg(3,t_POLMOD);
726: z[1]=lcopy(y);
727: if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
728: if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
729: }
730: err(gmoder1); return NULL; /* not reached */
731: }
732:
733: GEN
734: Mod0(GEN x,GEN y,long flag)
735: {
736: switch(flag)
737: {
738: case 0: return gmodulcp(x,y);
739: case 1: return gmodulo(x,y);
740: default: err(flagerr,"Mod");
741: }
742: return NULL; /* not reached */
743: }
744:
745: /*******************************************************************/
746: /* */
747: /* DIVISION ENTIERE GENERALE */
748: /* DIVISION ENTIERE AVEC RESTE GENERALE */
749: /* */
750: /*******************************************************************/
751:
752: GEN
753: gdivent(GEN x, GEN y)
754: {
755: long tx=typ(x),ty=typ(y);
756:
757: if (tx == t_INT)
758: {
759: if (ty == t_INT)
760: {
761: GEN z,p1;
762: long av,tetpil;
763:
764: if (signe(x)>=0) return divii(x,y);
765: av=avma; z=dvmdii(x,y,&p1);
766: if (p1==gzero) return z;
767: tetpil=avma; return gerepile(av,tetpil,addsi(-signe(y),z));
768: }
769: if (ty!=t_POL) err(typeer,"gdivent");
770: return gzero;
771: }
772: if (ty != t_POL) err(typeer,"gdivent");
773: if (tx == t_POL) return gdeuc(x,y);
774: if (! is_scalar_t(tx)) err(typeer,"gdivent");
775: return gzero;
776: }
777:
778: GEN
779: gdiventres(GEN x, GEN y)
780: {
781: long tx=typ(x),ty=typ(y);
782: GEN z = cgetg(3,t_COL);
783:
784: if (tx==t_INT)
785: {
786: if (ty==t_INT) z[1]=(long)truedvmdii(x,y,(GEN*)(z+2));
787: else
788: {
789: if (ty!=t_POL) err(typeer,"gdiventres");
790: z[1]=zero; z[2]=licopy(x);
791: }
792: return z;
793: }
794: if (ty != t_POL) err(typeer,"gdiventres");
795: if (tx == t_POL)
796: {
797: z[1]=ldivres(x,y,(GEN *)(z+2));
798: return z;
799: }
800: if (! is_scalar_t(tx)) err(typeer,"gdiventres");
801: z[1]=zero; z[2]=lcopy(x); return z;
802: }
803:
804: GEN
805: gdivmod(GEN x, GEN y, GEN *pr)
806: {
807: long ty,tx=typ(x);
808:
809: if (tx==t_INT)
810: {
811: ty=typ(y);
812: if (ty==t_INT) return dvmdii(x,y,pr);
813: if (ty==t_POL) { *pr=gcopy(x); return gzero; }
814: err(typeer,"gdivmod");
815: }
816: if (tx!=t_POL) err(typeer,"gdivmod");
817: return poldivres(x,y,pr);
818: }
819:
820: /* When x and y are integers, compute the quotient x/y, rounded to the
821: * nearest integer. If there is a tie, the quotient is rounded towards zero.
822: * If x and y are not both integers, same as gdivent.
823: */
824: GEN
825: gdivround(GEN x, GEN y)
826: {
827: long tx=typ(x),ty=typ(y);
828:
829: if (tx==t_INT)
830: {
831: if (ty==t_INT)
832: {
833: long av = avma, av1,fl;
834: GEN r, q=dvmdii(x,y,&r); /* q = x/y rounded towards 0, sqn(r)=sgn(x) */
835:
836: if (r==gzero) return q;
837: av1 = avma;
838: fl = absi_cmp(shifti(r,1),y);
839: avma = av1; cgiv(r);
840: if (fl >= 0) /* If 2*|r| >= |y| */
841: {
842: long sz = signe(x)*signe(y);
843: if (fl || sz > 0)
844: { av1=avma; q = gerepile(av,av1,addis(q,sz)); }
845: }
846: return q;
847: }
848: if (ty!=t_POL) err(typeer,"gdivround");
849: return gzero;
850: }
851: if (ty != t_POL) err(typeer,"gdivround");
852: if (tx == t_POL) return gdeuc(x,y);
853: if (! is_scalar_t(tx)) err(typeer,"gdivround");
854: return gzero;
855: }
856:
857: /*******************************************************************/
858: /* */
859: /* SHIFT D'UN GEN */
860: /* */
861: /*******************************************************************/
862:
863: /* Shift tronque si n<0 (multiplication tronquee par 2^n) */
864: GEN
865: gshift(GEN x, long n)
866: {
867: long i,l,lx, tx = typ(x);
868: GEN y;
869:
870: switch(tx)
871: {
872: case t_INT:
873: return shifti(x,n);
874: case t_REAL:
875: return shiftr(x,n);
876:
877: case t_VEC: case t_COL: case t_MAT:
878: lx=lg(x); y=cgetg(lx,tx); l=lontyp[tx];
879: for (i=1; i<l ; i++) y[i]=x[i];
880: for ( ; i<lx; i++) y[i]=lshift((GEN)x[i],n);
881: return y;
882: }
883: return gmul2n(x,n);
884: }
885:
886: GEN mulscalrfrac(GEN x, GEN y);
887:
888: /* Shift vrai (multiplication exacte par 2^n) */
889: GEN
890: gmul2n(GEN x, long n)
891: {
892: long tx,lx,i,k,l,av,tetpil;
893: GEN p2,p1,y;
894:
895: tx=typ(x);
896: switch(tx)
897: {
898: case t_INT:
899: if (n>=0) return shifti(x,n);
900: if (!signe(x)) return gzero;
901: l = vali(x); n = -n;
902: if (n<=l) return shifti(x,-n);
903: y=cgetg(3,t_FRAC);
904: y[1]=lshifti(x,-l);
905: y[2]=lshifti(gun,n-l); return y;
906:
907: case t_REAL:
908: return shiftr(x,n);
909:
910: case t_INTMOD:
911: if (n > 0)
912: {
913: y=cgetg(3,t_INTMOD); p2=(GEN)x[1];
914: av=avma; new_chunk(lgefint(p2) * (3 + (n>>TWOPOTBITS_IN_LONG)));
915: p1 = shifti((GEN)x[2],n); avma=av;
916: y[2]=lmodii(p1,p2); icopyifstack(p2,y[1]); return y;
917: }
918: l=avma; y=gmul2n(gun,n); tetpil=avma;
919: return gerepile(l,tetpil,gmul(y,x));
920:
921: case t_FRAC: case t_FRACN:
922: l = vali((GEN)x[1]);
923: k = vali((GEN)x[2]);
924: if (n+l-k>=0)
925: {
926: if (expi((GEN)x[2]) == k) /* x[2] power of 2 */
927: return shifti((GEN)x[1],n-k);
928: l = n-k; k = -k;
929: }
930: else
931: {
932: k = -l-n; l = -l;
933: }
934: y=cgetg(3,t_FRAC);
935: y[1]=lshifti((GEN)x[1],l);
936: y[2]=lshifti((GEN)x[2],k); return y;
937:
938: case t_QUAD: y=cgetg(4,t_QUAD);
939: copyifstack(x[1],y[1]);
940: for (i=2; i<4; i++) y[i]=lmul2n((GEN)x[i],n);
941: return y;
942:
943: case t_POLMOD: y=cgetg(3,t_POLMOD);
944: copyifstack(x[1],y[1]);
945: y[2]=lmul2n((GEN)x[2],n); return y;
946:
947: case t_POL: case t_COMPLEX: case t_SER:
948: case t_VEC: case t_COL: case t_MAT:
949: lx = (tx==t_POL)? lgef(x): lg(x);
950: y=cgetg(lx,tx); l=lontyp[tx];
951: for (i=1; i<l ; i++) y[i]=x[i];
952: for ( ; i<lx; i++) y[i]=lmul2n((GEN)x[i],n);
953: return y;
954:
955: case t_RFRAC: av=avma; p1 = gmul2n(gun,n); tetpil = avma;
956: return gerepile(av,tetpil, mulscalrfrac(p1,x));
957:
958: case t_RFRACN: y=cgetg(3,tx);
959: if (n>=0)
960: {
961: y[1]=lmul2n((GEN)x[1],n); y[2]=lcopy((GEN)x[2]);
962: }
963: else
964: {
965: y[2]=lmul2n((GEN)x[2],-n); y[1]=lcopy((GEN)x[1]);
966: }
967: return y;
968:
969: case t_PADIC:
970: l=avma; y=gmul2n(gun,n); tetpil=avma;
971: return gerepile(l,tetpil,gmul(y,x));
972: }
973: err(typeer,"gmul2n");
974: return NULL; /* not reached */
975: }
976:
977: /*******************************************************************/
978: /* */
979: /* INVERSE D'UN GEN */
980: /* */
981: /*******************************************************************/
982: GEN fix_rfrac_if_pol(GEN x, GEN y);
983:
984: GEN
985: ginv(GEN x)
986: {
987: long tx=typ(x),av,tetpil,s;
988: GEN z,y,p1,p2;
989:
990: switch(tx)
991: {
992: case t_INT:
993: if (is_pm1(x)) return icopy(x);
994: if (!signe(x)) err(gdiver2);
995: z=cgetg(3,t_FRAC);
996: z[1] = (signe(x)<0)? lnegi(gun): un;
997: z[2]=labsi(x); return z;
998:
999: case t_REAL:
1000: return divsr(1,x);
1001:
1002: case t_INTMOD: z=cgetg(3,t_INTMOD);
1003: icopyifstack(x[1],z[1]);
1004: z[2]=lmpinvmod((GEN)x[2],(GEN)x[1]); return z;
1005:
1006: case t_FRAC: case t_FRACN: z=cgetg(3,tx);
1007: s = signe(x[1]);
1008: if (!s) err(gdiver2);
1009: if (is_pm1(x[1]))
1010: return s>0? icopy((GEN)x[2]): negi((GEN)x[2]);
1011: z[1]=licopy((GEN)x[2]);
1012: z[2]=licopy((GEN)x[1]);
1013: if (s < 0)
1014: {
1015: setsigne(z[1],-signe(z[1]));
1016: setsigne(z[2],1);
1017: }
1018: return z;
1019:
1020: case t_COMPLEX: case t_QUAD:
1021: av=avma; p1=gnorm(x); p2=gconj(x); tetpil=avma;
1022: return gerepile(av,tetpil,gdiv(p2,p1));
1023:
1024: case t_PADIC: z = cgetg(5,t_PADIC);
1025: if (!signe(x[4])) err(gdiver2);
1026: z[1] = evalprecp(precp(x)) | evalvalp(-valp(x));
1027: icopyifstack(x[2], z[2]);
1028: z[3] = licopy((GEN)x[3]);
1029: z[4] = lmpinvmod((GEN)x[4],(GEN)z[3]); return z;
1030:
1031: case t_POLMOD: z=cgetg(3,t_POLMOD);
1032: copyifstack(x[1],z[1]);
1033: z[2]=linvmod((GEN)x[2],(GEN)x[1]); return z;
1034:
1035: case t_POL: case t_SER:
1036: return gdiv(gun,x);
1037:
1038: case t_RFRAC: case t_RFRACN:
1039: if (gcmp0((GEN) x[1])) err(gdiver2);
1040: p1 = fix_rfrac_if_pol((GEN)x[2],(GEN)x[1]);
1041: if (p1) return p1;
1042: z=cgetg(3,tx);
1043: z[1]=lcopy((GEN)x[2]);
1044: z[2]=lcopy((GEN)x[1]); return z;
1045:
1046: case t_QFR:
1047: {
1048: long k,l;
1049: l=signe(x[2]); setsigne(x[2],-l);
1050: k=signe(x[4]); setsigne(x[4],-k); z=redreal(x);
1051: setsigne(x[2],l); setsigne(x[4],k); return z;
1052: }
1053: case t_QFI:
1054: y=gcopy(x);
1055: if (!egalii((GEN)x[1],(GEN)x[2]) && !egalii((GEN)x[1],(GEN)x[3]))
1056: setsigne(y[2],-signe(y[2]));
1057: return y;
1058: case t_MAT:
1059: return (lg(x)==1)? cgetg(1,t_MAT): invmat(x);
1060: }
1061: err(typeer,"inverse");
1062: return NULL; /* not reached */
1063: }
1064:
1065: /*******************************************************************/
1066: /* */
1067: /* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
1068: /* */
1069: /*******************************************************************/
1070:
1071: /* Convert t_SER --> t_POL / t_RFRAC */
1072: static GEN
1073: gconvsp(GEN x, int flpile)
1074: {
1075: long v = varn(x), av,tetpil,i;
1076: GEN p1,y;
1077:
1078: if (gcmp0(x)) return zeropol(v);
1079: av=avma; y=dummycopy(x); settyp(y,t_POL);
1080: i=lg(x)-1; while (i>1 && gcmp0((GEN)y[i])) i--;
1081: setlgef(y,i+1);
1082: p1=gpuigs(polx[v],valp(x));
1083: tetpil=avma; p1=gmul(p1,y);
1084: return flpile? gerepile(av,tetpil,p1): p1;
1085: }
1086:
1087: GEN
1088: gsubst(GEN x, long v, GEN y)
1089: {
1090: long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
1091: long l,vx,vy,e,ex,ey,tetpil,av,i,j,k,jb,limite;
1092: GEN t,p1,p2,z;
1093:
1094: if (ty==t_MAT)
1095: {
1096: if (ly==1) return cgetg(1,t_MAT);
1097: if (ly != lg(y[1]))
1098: err(talker,"forbidden substitution by a non square matrix");
1099: } else if (is_graphicvec_t(ty))
1100: err(talker,"forbidden substitution by a vector");
1101:
1102: if (is_scalar_t(tx))
1103: {
1104: if (tx!=t_POLMOD || v <= varn(x[1]))
1105: {
1106: if (ty==t_MAT) return gscalmat(x,ly-1);
1107: return gcopy(x);
1108: }
1109: av=avma;
1110: p1=gsubst((GEN)x[1],v,y); vx=varn(p1);
1111: p2=gsubst((GEN)x[2],v,y); vy=gvar(p2);
1112: if (typ(p1)!=t_POL)
1113: err(talker,"forbidden substitution in a scalar type");
1114: if (vy>=vx)
1115: {
1116: tetpil=avma;
1117: return gerepile(av,tetpil,gmodulcp(p2,p1));
1118: }
1119: lx=lgef(p2); tetpil=avma; z=cgetg(lx,t_POL); z[1]=p2[1];
1120: for (i=2; i<lx; i++) z[i]=lmodulcp((GEN)p2[i],p1);
1121: return gerepile(av,tetpil, normalizepol_i(z,lx));
1122: }
1123:
1124: switch(tx)
1125: {
1126: case t_POL:
1127: l=lgef(x);
1128: if (l==2)
1129: return (ty==t_MAT)? gscalmat(gzero,ly-1): gzero;
1130:
1131: vx=varn(x);
1132: if (vx<v)
1133: {
1134: av=avma; p1=polx[vx]; z= gsubst((GEN)x[l-1],v,y);
1135: for (i=l-1; i>2; i--) z=gadd(gmul(z,p1),gsubst((GEN)x[i-1],v,y));
1136: return gerepileupto(av,z);
1137: }
1138: if (ty!=t_MAT)
1139: return (vx>v)? gcopy(x): poleval(x,y);
1140:
1141: if (vx>v) return gscalmat(x,ly-1);
1142: if (l==3) return gscalmat((GEN)x[2],ly-1);
1143: av=avma; z=(GEN)x[l-1];
1144: for (i=l-1; i>2; i--) z=gaddmat((GEN)x[i-1],gmul(z,y));
1145: return gerepileupto(av,z);
1146:
1147: case t_SER:
1148: vx=varn(x);
1149: if (vx > v)
1150: return (ty==t_MAT)? gscalmat(x,ly-1): gcopy(x);
1151: ex = valp(x);
1152: if (vx < v)
1153: {
1154: if (!signe(x)) return gcopy(x);
1155: /* a ameliorer */
1156: av=avma; setvalp(x,0); p1=gconvsp(x,0); setvalp(x,ex);
1157: p2=gsubst(p1,v,y); tetpil=avma; z=tayl(p2,vx,lx-2);
1158: if (ex)
1159: {
1160: p1=gpuigs(polx[vx],ex); tetpil=avma; z=gmul(z,p1);
1161: }
1162: return gerepile(av,tetpil,z);
1163: }
1164: switch(ty) /* here vx == v */
1165: {
1166: case t_SER:
1167: ey=valp(y); vy=varn(y);
1168: if (ey<1) return zeroser(vy,ey*(ex+lx-2));
1169: l=(lx-2)*ey+2;
1170: if (ex) { if (l>ly) l=ly; }
1171: else
1172: {
1173: if (gcmp0(y)) l=ey+2;
1174: else { if (l>ly) l=ey+ly; }
1175: }
1176: if (vy!=vx)
1177: {
1178: av=avma; z = zeroser(vy,0);
1179: for (i=lx-1; i>=2; i--)
1180: z = gadd((GEN)x[i], gmul(y,z));
1181: if (ex) z = gmul(z, gpuigs(y,ex));
1182: return gerepileupto(av,z);
1183: }
1184:
1185: av=avma; limite=stack_lim(av,1);
1186: t=cgetg(ly,t_SER);
1187: z = scalarser((GEN)x[2],varn(y),l-2);
1188: for (i=2; i<ly; i++) t[i]=y[i];
1189:
1190: for (i=3,jb=ey; jb<=l-2; i++,jb+=ey)
1191: {
1192: for (j=jb+2; j<l; j++)
1193: {
1194: p1=gmul((GEN)x[i],(GEN)t[j-jb]);
1195: z[j]=ladd((GEN)z[j],p1);
1196: }
1197: for (j=l-1-jb-ey; j>1; j--)
1198: {
1199: p1=gzero;
1200: for (k=2; k<j; k++)
1201: p1=gadd(p1,gmul((GEN)t[j-k+2],(GEN)y[k]));
1202: p2=gmul((GEN)t[2],(GEN)y[j]);
1203: t[j]=ladd(p1,p2);
1204: }
1205: if (low_stack(limite, stack_lim(av,1)))
1206: {
1207: GEN *gptr[2];
1208: if(DEBUGMEM>1) err(warnmem,"gsubst");
1209: gptr[0]=&z; gptr[1]=&t; gerepilemany(av,gptr,2);
1210: }
1211: }
1212: if (!ex) { tetpil=avma; return gerepile(av,tetpil,gcopy(z)); }
1213:
1214: if (l<ly) { setlg(y,l); p1=gpuigs(y,ex); setlg(y,ly); }
1215: else p1=gpuigs(y,ex);
1216: tetpil=avma; return gerepile(av,tetpil,gmul(z,p1));
1217:
1218: case t_POL: case t_RFRAC: case t_RFRACN:
1219: if (isexactzero(y)) return scalarser((GEN)x[2],v,lx-2);
1220: vy=gvar(y); e=gval(y,vy);
1221: if (e<=0)
1222: err(talker,"non positive valuation in a series substitution");
1223: av=avma; p1=gconvsp(x,0); p2=gsubst(p1,v,y); tetpil=avma;
1224: return gerepile(av,tetpil,tayl(p2,vy,e*(lx-2+ex)));
1225:
1226: default:
1227: err(talker,"non polynomial or series type substituted in a series");
1228: }
1229: break;
1230:
1231: case t_RFRAC: case t_RFRACN: av=avma;
1232: p1=gsubst((GEN)x[1],v,y);
1233: p2=gsubst((GEN)x[2],v,y); tetpil=avma;
1234: return gerepile(av,tetpil,gdiv(p1,p2));
1235:
1236: case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
1237: for (i=1; i<lx; i++) z[i]=lsubst((GEN)x[i],v,y);
1238: }
1239: return z;
1240: }
1241:
1242: /*******************************************************************/
1243: /* */
1244: /* SERIE RECIPROQUE D'UNE SERIE */
1245: /* */
1246: /*******************************************************************/
1247:
1248: GEN
1249: recip(GEN x)
1250: {
1251: long tetpil, av=avma, v=varn(x);
1252: GEN p1,p2,a,y,u;
1253:
1254: if (typ(x)!=t_SER) err(talker,"not a series in serreverse");
1255: if (valp(x)!=1) err(talker,"valuation not equal to 1 in serreverse");
1256:
1257: a=(GEN)x[2];
1258: if (gcmp1(a))
1259: {
1260: long i,j,k, lx=lg(x), lim=stack_lim(av,2);
1261:
1262: u=cgetg(lx,t_SER); y=cgetg(lx,t_SER);
1263: u[1]=y[1]=evalsigne(1) | evalvalp(1) | evalvarn(v);
1264: u[2]=un; u[3]=lmulsg(-2,(GEN)x[3]);
1265: y[2]=un; y[3]=lneg((GEN)x[3]);
1266: for (i=3; i<lx-1; )
1267: {
1268: for (j=3; j<i+1; j++)
1269: {
1270: p1=(GEN)u[j];
1271: for (k=j-1; k>2; k--)
1272: p1=gsub(p1,gmul((GEN)u[k],(GEN)x[j-k+2]));
1273: u[j]=lsub(p1,(GEN)x[j]);
1274: }
1275: p1=gmulsg(i,(GEN)x[i+1]);
1276: for (k=2; k<i; k++)
1277: {
1278: p2=gmul((GEN)x[k+1],(GEN)u[i-k+2]);
1279: p1=gadd(p1,gmulsg(k,p2));
1280: }
1281: i++; u[i]=lneg(p1); y[i]=ldivgs((GEN)u[i],i-1);
1282: if (low_stack(lim, stack_lim(av,2)))
1283: {
1284: GEN *gptr[2];
1285: if(DEBUGMEM>1) err(warnmem,"recip");
1286: for(k=i+1; k<lx; k++) u[k]=y[k]=zero; /* dummy */
1287: gptr[0]=&u; gptr[1]=&y; gerepilemany(av,gptr,2);
1288: }
1289: }
1290: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
1291: }
1292: y=gdiv(x,a); y[2]=un; y=recip(y);
1293: a=gdiv(polx[v],a); tetpil=avma;
1294: return gerepile(av,tetpil,gsubst(y,v,a));
1295: }
1296:
1297: /*******************************************************************/
1298: /* */
1299: /* DERIVATION ET INTEGRATION */
1300: /* */
1301: /*******************************************************************/
1302: GEN
1303: derivpol(GEN x)
1304: {
1305: long i,lx = lgef(x)-1;
1306: GEN y;
1307:
1308: if (lx<3) return gzero;
1309: y = cgetg(lx,t_POL);
1310: for (i=2; i<lx ; i++) y[i] = lmulsg(i-1,(GEN)x[i+1]);
1311: y[1] = x[1]; return normalizepol_i(y,i);
1312: }
1313:
1314: GEN
1315: derivser(GEN x)
1316: {
1317: long i,j,vx = varn(x), e = valp(x), lx = lg(x);
1318: GEN y;
1319: if (gcmp0(x)) return zeroser(vx,e-1);
1320: if (e)
1321: {
1322: y=cgetg(lx,t_SER); y[1] = evalsigne(1) | evalvalp(e-1) | evalvarn(vx);
1323: for (i=2; i<lx; i++) y[i]=lmulsg(i+e-2,(GEN)x[i]);
1324: return y;
1325: }
1326: i=3; while (i<lx && gcmp0((GEN)x[i])) i++;
1327: if (i==lx) return zeroser(vx,lx-3);
1328: lx--; if (lx<3) lx=3;
1329: lx = lx-i+3; y=cgetg(lx,t_SER);
1330: y[1]=evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
1331: for (j=2; j<lx; j++) y[j]=lmulsg(j+i-4,(GEN)x[i+j-2]);
1332: return y;
1333: }
1334:
1335: GEN
1336: deriv(GEN x, long v)
1337: {
1338: long lx,vx,tx,e,i,j,l,av,tetpil;
1339: GEN y,p1,p2;
1340:
1341: tx=typ(x); if (is_const_t(tx)) return gzero;
1342: if (v < 0) v = gvar(x);
1343: switch(tx)
1344: {
1345: case t_POLMOD:
1346: if (v<=varn(x[1])) return gzero;
1347: y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
1348: y[2]=lderiv((GEN)x[2],v); return y;
1349:
1350: case t_POL:
1351: vx=varn(x); if (vx>v) return gzero;
1352: if (vx<v)
1353: {
1354: lx = lgef(x);
1355: y = cgetg(lx,t_POL);
1356: for (i=2; i<lx; i++) y[i] = lderiv((GEN)x[i],v);
1357: y[1] = evalvarn(vx);
1358: return normalizepol_i(y,i);
1359: }
1360: return derivpol(x);
1361:
1362: case t_SER:
1363: vx=varn(x); if (vx>v) return gzero;
1364: if (vx<v)
1365: {
1366: if (!signe(x)) return gcopy(x);
1367: lx=lg(x); e=valp(x);
1368: l=avma;
1369: for (i=2; i<lx; i++)
1370: {
1371: if (!gcmp0(deriv((GEN)x[i],v))) break;
1372: avma=l;
1373: }
1374: if (i==lx) return ggrando(polx[vx],e+lx-2);
1375: y=cgetg(lx-i+2,t_SER);
1376: y[1] = evalsigne(1) | evalvalp(e+i-2) | evalvarn(vx);
1377: for (j=2; i<lx; j++,i++) y[j]=lderiv((GEN)x[i],v);
1378: return y;
1379: }
1380: return derivser(x);
1381:
1382: case t_RFRAC: case t_RFRACN: av=avma; y=cgetg(3,tx);
1383: y[2]=lsqr((GEN)x[2]); l=avma;
1384: p1=gmul((GEN)x[2],deriv((GEN)x[1],v));
1385: p2=gmul(gneg_i((GEN)x[1]),deriv((GEN)x[2],v));
1386: tetpil=avma; p1=gadd(p1,p2);
1387: if (tx==t_RFRACN) { y[1]=lpile(l,tetpil,p1); return y; }
1388: y[1]=(long)p1; return gerepileupto(av,gred_rfrac(y));
1389:
1390: case t_VEC: case t_COL: case t_MAT: lx=lg(x); y=cgetg(lx,tx);
1391: for (i=1; i<lx; i++) y[i]=lderiv((GEN)x[i],v);
1392: return y;
1393: }
1394: err(typeer,"deriv");
1395: return NULL; /* not reached */
1396: }
1397:
1398: /*******************************************************************/
1399: /* */
1400: /* INTEGRATION FORMELLE */
1401: /* */
1402: /*******************************************************************/
1403:
1404: static GEN
1405: triv_integ(GEN x, long v, long tx, long lx)
1406: {
1407: GEN y=cgetg(lx,tx);
1408: long i;
1409:
1410: y[1]=x[1];
1411: for (i=2; i<lx; i++) y[i]=linteg((GEN)x[i],v);
1412: return y;
1413: }
1414:
1415: GEN
1416: integ(GEN x, long v)
1417: {
1418: long lx,tx,e,i,j,vx,n,av=avma,tetpil;
1419: GEN y,p1;
1420:
1421: tx = typ(x);
1422: if (v < 0) v = gvar(x);
1423: if (is_scalar_t(tx))
1424: {
1425: if (tx == t_POLMOD && v>varn(x[1]))
1426: {
1427: y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
1428: y[2]=linteg((GEN)x[2],v); return y;
1429: }
1430: if (gcmp0(x)) return gzero;
1431:
1432: y=cgetg(4,t_POL); y[1] = evalsigne(1) | evallgef(4) | evalvarn(v);
1433: y[2]=zero; y[3]=lcopy(x); return y;
1434: }
1435:
1436: switch(tx)
1437: {
1438: case t_POL:
1439: vx=varn(x); lx=lgef(x);
1440: if (lx==2) return zeropol(min(v,vx));
1441: if (vx>v)
1442: {
1443: y=cgetg(4,t_POL);
1444: y[1] = signe(x)? evallgef(4) | evalvarn(v) | evalsigne(1)
1445: : evallgef(4) | evalvarn(v);
1446: y[2]=zero; y[3]=lcopy(x); return y;
1447: }
1448: if (vx<v) return triv_integ(x,v,tx,lx);
1449: y=cgetg(lx+1,tx); y[2]=zero;
1450: for (i=3; i<=lx; i++)
1451: y[i]=ldivgs((GEN)x[i-1],i-2);
1452: y[1] = signe(x)? evallgef(lx+1) | evalvarn(v) | evalsigne(1)
1453: : evallgef(lx+1) | evalvarn(v);
1454: return y;
1455:
1456: case t_SER:
1457: lx=lg(x); e=valp(x); vx=varn(x);
1458: if (!signe(x))
1459: {
1460: if (vx == v) e++; else if (vx < v) v = vx;
1461: return zeroser(v,e);
1462: }
1463: if (vx>v)
1464: {
1465: y=cgetg(4,t_POL);
1466: y[1] = evallgef(4) | evalvarn(v) | evalsigne(1);
1467: y[2]=zero; y[3]=lcopy(x); return y;
1468: }
1469: if (vx<v) return triv_integ(x,v,tx,lx);
1470: y=cgetg(lx,tx);
1471: for (i=2; i<lx; i++)
1472: {
1473: j=i+e-1;
1474: if (!j)
1475: {
1476: if (!gcmp0((GEN)x[i])) err(inter2);
1477: y[i]=zero;
1478: }
1479: else y[i] = ldivgs((GEN)x[i],j);
1480: }
1481: y[1]=x[1]+1; return y;
1482:
1483: case t_RFRAC: case t_RFRACN:
1484: vx = gvar(x);
1485: if (vx>v)
1486: {
1487: y=cgetg(4,t_POL);
1488: y[1] = signe(x[1])? evallgef(4) | evalvarn(v) | evalsigne(1)
1489: : evallgef(4) | evalvarn(v);
1490: y[2]=zero; y[3]=lcopy(x); return y;
1491: }
1492: if (vx<v)
1493: {
1494: p1=cgetg(v+2,t_VEC);
1495: for (i=0; i<vx; i++) p1[i+1] = lpolx[i];
1496: for (i=vx+1; i<v; i++) p1[i+1] = lpolx[i];
1497: p1[v+1] = lpolx[vx]; p1[vx+1] = lpolx[v];
1498: y=integ(changevar(x,p1),vx); tetpil=avma;
1499: return gerepile(av,tetpil,changevar(y,p1));
1500: }
1501:
1502: n=lgef(x[1])+lgef(x[2])-4;
1503: y=gdiv(gtrunc(gmul((GEN)x[2], integ(tayl(x,v,n),v))), (GEN)x[2]);
1504: if (!gegal(deriv(y,v),x)) err(inter2);
1505: if (typ(y)==t_RFRAC && lgef(y[1]) == lgef(y[2]))
1506: {
1507: GEN p2;
1508: tx=typ(y[1]); p1=is_scalar_t(tx)? (GEN)y[1]: leading_term(y[1]);
1509: tx=typ(y[2]); p2=is_scalar_t(tx)? (GEN)y[2]: leading_term(y[2]);
1510: y=gsub(y, gdiv(p1,p2));
1511: }
1512: return gerepileupto(av,y);
1513:
1514: case t_VEC: case t_COL: case t_MAT:
1515: lx=lg(x); y=cgetg(lx,tx);
1516: for (i=1; i<lg(x); i++) y[i]=linteg((GEN)x[i],v);
1517: return y;
1518: }
1519: err(typeer,"integ");
1520: return NULL; /* not reached */
1521: }
1522:
1523: /*******************************************************************/
1524: /* */
1525: /* PARTIES ENTIERES */
1526: /* */
1527: /*******************************************************************/
1528:
1529: GEN
1530: gfloor(GEN x)
1531: {
1532: GEN y;
1533: long i,lx, tx = typ(x);
1534:
1535: switch(tx)
1536: {
1537: case t_INT: case t_POL:
1538: return gcopy(x);
1539:
1540: case t_REAL:
1541: return mpent(x);
1542:
1543: case t_FRAC: case t_FRACN:
1544: return truedvmdii((GEN)x[1],(GEN)x[2],NULL);
1545:
1546: case t_RFRAC: case t_RFRACN:
1547: return gdeuc((GEN)x[1],(GEN)x[2]);
1548:
1549: case t_VEC: case t_COL: case t_MAT:
1550: lx=lg(x); y=cgetg(lx,tx);
1551: for (i=1; i<lx; i++) y[i]=lfloor((GEN)x[i]);
1552: return y;
1553: }
1554: err(typeer,"gfloor");
1555: return NULL; /* not reached */
1556: }
1557:
1558: GEN
1559: gfrac(GEN x)
1560: {
1561: long av = avma,tetpil;
1562: GEN p1 = gneg_i(gfloor(x));
1563:
1564: tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
1565: }
1566:
1567: GEN
1568: gceil(GEN x)
1569: {
1570: GEN y,p1;
1571: long i,lx,tx=typ(x),av,tetpil;
1572:
1573: switch(tx)
1574: {
1575: case t_INT: case t_POL:
1576: return gcopy(x);
1577:
1578: case t_REAL:
1579: av=avma; y=mpent(x);
1580: if (!gegal(x,y))
1581: {
1582: tetpil=avma; return gerepile(av,tetpil,addsi(1,y));
1583: }
1584: return y;
1585:
1586: case t_FRAC: case t_FRACN:
1587: av=avma; y=dvmdii((GEN)x[1],(GEN)x[2],&p1);
1588: if (p1 != gzero && gsigne(x)>0)
1589: {
1590: cgiv(p1); tetpil=avma;
1591: return gerepile(av,tetpil,addsi(1,y));
1592: }
1593: return y;
1594:
1595: case t_RFRAC: case t_RFRACN:
1596: return gdeuc((GEN)x[1],(GEN)x[2]);
1597:
1598: case t_VEC: case t_COL: case t_MAT:
1599: lx=lg(x); y=cgetg(lx,tx);
1600: for (i=1; i<lx; i++) y[i]=lceil((GEN)x[i]);
1601: return y;
1602: }
1603: err(typeer,"gceil");
1604: return NULL; /* not reached */
1605: }
1606:
1607: GEN
1608: round0(GEN x, GEN *pte)
1609: {
1610: if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
1611: return ground(x);
1612: }
1613:
1614: GEN
1615: ground(GEN x)
1616: {
1617: GEN y,p1;
1618: long i,lx,tx=typ(x),av,tetpil;
1619:
1620: switch(tx)
1621: {
1622: case t_INT: case t_INTMOD: case t_QUAD:
1623: return gcopy(x);
1624:
1625: case t_REAL:
1626: {
1627: long ex, s = signe(x);
1628: if (s==0 || (ex=expo(x)) < -1) return gzero;
1629: if (ex < 0) return s>0? gun: negi(gun);
1630: av=avma; p1 = realun(3 + (ex>>TWOPOTBITS_IN_LONG));
1631: setexpo(p1,-1); /* p1 = 0.5 */
1632: p1 = addrr(x,p1); tetpil = avma;
1633: return gerepile(av,tetpil,mpent(p1));
1634: }
1635: case t_FRAC: case t_FRACN:
1636: av=avma; p1 = addii(shifti((GEN)x[2], -1), (GEN)x[1]);
1637: return gerepileuptoint(av, truedvmdii(p1, (GEN)x[2], NULL));
1638:
1639: case t_POLMOD: y=cgetg(3,t_POLMOD);
1640: copyifstack(x[1],y[1]);
1641: y[2]=lround((GEN)x[2]); return y;
1642:
1643: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
1644: case t_VEC: case t_COL: case t_MAT:
1645: lx = (tx==t_POL)? lgef(x): lg(x);
1646: y=cgetg(lx,tx);
1647: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
1648: for ( ; i<lx; i++) y[i]=lround((GEN)x[i]);
1649: if (tx==t_POL) return normalizepol_i(y, lx);
1650: if (tx==t_SER) return normalize(y);
1651: return y;
1652: }
1653: err(typeer,"ground");
1654: return NULL; /* not reached */
1655: }
1656:
1657: /* e = number of error bits on integral part */
1658: GEN
1659: grndtoi(GEN x, long *e)
1660: {
1661: GEN y,p1;
1662: long i,tx=typ(x), lx,av,ex,e1;
1663:
1664: *e = -HIGHEXPOBIT;
1665: switch(tx)
1666: {
1667: case t_INT: case t_INTMOD: case t_QUAD:
1668: case t_FRAC: case t_FRACN:
1669: return ground(x);
1670:
1671: case t_REAL:
1672: av=avma; p1=gadd(ghalf,x); ex=expo(p1);
1673: if (ex<0)
1674: {
1675: if (signe(p1)>=0) { *e=expo(x); avma=av; return gzero; }
1676: *e=expo(addsr(1,x)); avma=av; return negi(gun);
1677: }
1678: lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
1679: settyp(p1,t_INT); setlgefint(p1,lx);
1680: y=shifti(p1,e1); if (signe(x)<0) y=addsi(-1,y);
1681: y = gerepileupto(av,y);
1682:
1683: if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
1684: *e=e1; return y;
1685:
1686: case t_POLMOD: y=cgetg(3,t_POLMOD);
1687: copyifstack(x[1],y[1]);
1688: y[2]=lrndtoi((GEN)x[2],e); return y;
1689:
1690: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
1691: case t_VEC: case t_COL: case t_MAT:
1692: lx=(tx==t_POL)? lgef(x): lg(x);
1693: y=cgetg(lx,tx);
1694: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
1695: for ( ; i<lx; i++)
1696: {
1697: y[i]=lrndtoi((GEN)x[i],&e1);
1698: if (e1>*e) *e=e1;
1699: }
1700: return y;
1701: }
1702: err(typeer,"grndtoi");
1703: return NULL; /* not reached */
1704: }
1705:
1706: /* e = number of error bits on integral part */
1707: GEN
1708: gcvtoi(GEN x, long *e)
1709: {
1710: long tx=typ(x), lx,i,ex,av,e1;
1711: GEN y;
1712:
1713: *e = -HIGHEXPOBIT;
1714: if (tx == t_REAL)
1715: {
1716: long x0,x1;
1717: ex=expo(x); if (ex<0) { *e=ex; return gzero; }
1718: lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
1719: x0=x[0]; x1=x[1]; settyp(x,t_INT); setlgefint(x,lx);
1720: y=shifti(x,e1); x[0]=x0; x[1]=x1;
1721: if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
1722: *e=e1; return y;
1723: }
1724: if (is_matvec_t(tx))
1725: {
1726: lx=lg(x); y=cgetg(lx,tx);
1727: for (i=1; i<lx; i++)
1728: {
1729: y[i]=lcvtoi((GEN)x[i],&e1);
1730: if (e1>*e) *e=e1;
1731: }
1732: return y;
1733: }
1734: return gtrunc(x);
1735: }
1736:
1737: GEN
1738: gtrunc(GEN x)
1739: {
1740: long tx=typ(x),av,tetpil,i,v;
1741: GEN y;
1742:
1743: switch(tx)
1744: {
1745: case t_INT: case t_POL:
1746: return gcopy(x);
1747:
1748: case t_REAL:
1749: return mptrunc(x);
1750:
1751: case t_FRAC: case t_FRACN:
1752: return divii((GEN)x[1],(GEN)x[2]);
1753:
1754: case t_PADIC:
1755: if (!signe(x[4])) return gzero;
1756: v = valp(x);
1757: if (!v) return gcopy((GEN)x[4]);
1758: if (v>0)
1759: { /* here p^v is an integer */
1760: av=avma; y=gpuigs((GEN)x[2],v); tetpil=avma;
1761: return gerepile(av,tetpil, mulii(y,(GEN)x[4]));
1762: }
1763: y=cgetg(3,t_FRAC);
1764: y[1]=licopy((GEN)x[4]);
1765: y[2]=lpuigs((GEN)x[2],-v);
1766: return y;
1767:
1768: case t_RFRAC: case t_RFRACN:
1769: return gdeuc((GEN)x[1],(GEN)x[2]);
1770:
1771: case t_SER:
1772: return gconvsp(x,1);
1773:
1774: case t_VEC: case t_COL: case t_MAT:
1775: {
1776: long lx=lg(x);
1777:
1778: y=cgetg(lx,tx);
1779: for (i=1; i<lx; i++) y[i]=ltrunc((GEN)x[i]);
1780: return y;
1781: }
1782: }
1783: err(typeer,"gtrunc");
1784: return NULL; /* not reached */
1785: }
1786:
1787: GEN
1788: trunc0(GEN x, GEN *pte)
1789: {
1790: if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
1791: return gtrunc(x);
1792: }
1793:
1794: /*******************************************************************/
1795: /* */
1796: /* CONVERSIONS GEN --> POLYNOMIALS & SERIES */
1797: /* */
1798: /*******************************************************************/
1799:
1800: GEN
1801: zeropol(long v)
1802: {
1803: GEN x = cgetg(2,t_POL);
1804: x[1] = evallgef(2) | evalvarn(v); return x;
1805: }
1806:
1807: GEN
1808: scalarpol(GEN x, long v)
1809: {
1810: GEN y=cgetg(3,t_POL);
1811: y[1] = gcmp0(x)? evallgef(3) | evalvarn(v)
1812: : evallgef(3) | evalvarn(v) | evalsigne(1);
1813: y[2]=lcopy(x); return y;
1814: }
1815:
1816: static GEN
1817: gtopoly0(GEN x, long v, int reverse)
1818: {
1819: long tx=typ(x),lx,i,j;
1820: GEN y;
1821:
1822: if (v<0) v = 0;
1823: if (isexactzero(x)) return zeropol(v);
1824: if (is_scalar_t(tx)) return scalarpol(x,v);
1825:
1826: switch(tx)
1827: {
1828: case t_POL:
1829: y=gcopy(x); break;
1830: case t_SER:
1831: y=gconvsp(x,1); break;
1832: case t_RFRAC: case t_RFRACN:
1833: y=gdeuc((GEN)x[1],(GEN)x[2]); break;
1834: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
1835: lx=lg(x);
1836: if (reverse)
1837: {
1838: while (lx-- && isexactzero((GEN)x[lx]));
1839: i=lx+2; y=cgetg(i,t_POL);
1840: y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
1841: for (j=2; j<i; j++) y[j]=lcopy((GEN)x[j-1]);
1842: }
1843: else
1844: {
1845: i=1; j=lx; while (lx-- && isexactzero((GEN)x[i++]));
1846: i=lx+2; y=cgetg(i,t_POL);
1847: y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
1848: lx = j-1;
1849: for (j=2; j<i; j++) y[j]=lcopy((GEN)x[lx--]);
1850: }
1851: break;
1852: default: err(typeer,"gtopoly");
1853: }
1854: setvarn(y,v); return y;
1855: }
1856:
1857: GEN
1858: gtopolyrev(GEN x, long v) { return gtopoly0(x,v,1); }
1859:
1860: GEN
1861: gtopoly(GEN x, long v) { return gtopoly0(x,v,0); }
1862:
1863: GEN
1864: zeroser(long v, long val)
1865: {
1866: GEN x = cgetg(2,t_SER);
1867: x[1]=evalvalp(val) | evalvarn(v); return x;
1868: }
1869:
1870: GEN
1871: scalarser(GEN x, long v, long prec)
1872: {
1873: GEN y=cgetg(prec+2,t_SER);
1874: long i;
1875:
1876: y[1]=evalsigne(1) | evalvalp(0) | evalvarn(v);
1877: y[2]=lcopy(x); for (i=3; i<=prec+1; i++) y[i]=zero;
1878: return y;
1879: }
1880:
1881: GEN
1882: gtoser(GEN x, long v)
1883: {
1884: long tx=typ(x),lx,i,j,l,av,tetpil;
1885: GEN y,p1,p2;
1886:
1887: if (v<0) v = 0;
1888: if (tx==t_SER) { y=gcopy(x); setvarn(y,v); return y; }
1889: if (isexactzero(x)) return zeroser(v,precdl);
1890: if (is_scalar_t(tx)) return scalarser(x,v,precdl);
1891:
1892: switch(tx)
1893: {
1894: case t_POL:
1895: lx=lgef(x); i=2; while (i<lx && gcmp0((GEN)x[i])) i++;
1896: l=lx-i; if (precdl>l) l=precdl;
1897: y=cgetg(l+2,t_SER);
1898: y[1] = evalsigne(1) | evalvalp(i-2) | evalvarn(v);
1899: for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
1900: for ( ; j<=l+1; j++) y[j]=zero;
1901: break;
1902:
1903: case t_RFRAC: case t_RFRACN:
1904: av=avma; p1=gtoser((GEN)x[1],v); p2=gtoser((GEN)x[2],v);
1905: tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
1906:
1907: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
1908: lx=lg(x); i=1; while (i<lx && isexactzero((GEN)x[i])) i++;
1909: y = cgetg(lx-i+2,t_SER);
1910: y[1] = evalsigne(1) | evalvalp(i-1) | evalvarn(v);
1911: for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
1912: break;
1913:
1914: default: err(typeer,"gtoser");
1915: }
1916: return y;
1917: }
1918:
1919: GEN
1920: gtovec(GEN x)
1921: {
1922: long tx,lx,i;
1923: GEN y;
1924:
1925: if (!x) return cgetg(1,t_VEC);
1926: tx = typ(x);
1927: if (is_scalar_t(tx) || is_rfrac_t(tx))
1928: {
1929: y=cgetg(2,t_VEC); y[1]=lcopy(x);
1930: return y;
1931: }
1932: if (is_graphicvec_t(tx))
1933: {
1934: lx=lg(x); y=cgetg(lx,t_VEC);
1935: for (i=1; i<lx; i++) y[i]=lcopy((GEN)x[i]);
1936: return y;
1937: }
1938: if (tx==t_POL)
1939: {
1940: lx=lgef(x); y=cgetg(lx-1,t_VEC);
1941: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[lx-i]);
1942: return y;
1943: }
1944: if (tx==t_LIST)
1945: {
1946: lx=lgef(x); y=cgetg(lx-1,t_VEC); x++;
1947: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
1948: return y;
1949: }
1950: if (!signe(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
1951: lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
1952: for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
1953: return y;
1954: }
1955:
1956: GEN
1957: compo(GEN x, long n)
1958: {
1959: long l,tx=typ(x);
1960:
1961: if (tx==t_POL && n+1 >= lgef(x)) return gzero;
1962: if (tx==t_SER && !signe(x)) return gzero;
1963: if (!is_recursive_t(tx))
1964: err(talker, "this object doesn't have components (is a leaf)");
1965: l=lontyp[tx]+n-1;
1966: if (n<1 || l>=lg(x))
1967: err(talker,"nonexistent component");
1968: return gcopy((GEN)x[l]);
1969: }
1970:
1971: /* with respect to the main variable if v<0, with respect to the variable v
1972: otherwise. v ignored if x is not a polynomial/series. */
1973:
1974: GEN
1975: polcoeff0(GEN x, long n, long v)
1976: {
1977: long tx=typ(x),lx,ex,w,av,tetpil;
1978: GEN xinit;
1979:
1980: if (is_scalar_t(tx)) return n? gzero: gcopy(x);
1981:
1982: switch(tx)
1983: {
1984: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
1985: if (n<1 || n>=lg(x)) break;
1986: return gcopy((GEN)x[n]);
1987:
1988: case t_POL:
1989: if (n<0) return gzero;
1990: w=varn(x);
1991: if (v < 0 || v == w)
1992: return (n>=lgef(x)-2)? gzero: gcopy((GEN)x[n+2]);
1993: if (v < w) return n? gzero: gcopy(x);
1994: av=avma; xinit=x;
1995: x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
1996: if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
1997: if (typ(x) == t_POL)
1998: {
1999: if (n>=lgef(x)-2) { avma=av; return gzero; }
2000: x = (GEN)x[n+2];
2001: }
2002: else
2003: x = polcoeff0(x, n, 0);
2004: tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
2005:
2006: case t_SER:
2007: w=varn(x);
2008: if (v < 0 || v == w)
2009: {
2010: if (!signe(x)) return gzero;
2011: lx=lg(x); ex=valp(x); if (n<ex) return gzero;
2012: if (n>=ex+lx-2) break;
2013: return gcopy((GEN)x[n-ex+2]);
2014: }
2015: if (v < w) return n? gzero: gcopy(x);
2016: av=avma; xinit=x;
2017: x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
2018: if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
2019: if (gcmp0(x)) { avma=av; return gzero; }
2020: if (typ(x) == t_SER)
2021: {
2022: lx=lg(x); ex=valp(x); if (n<ex) return gzero;
2023: if (n>=ex+lx-2) break;
2024: x = (GEN)x[n-ex+2];
2025: }
2026: else
2027: x = polcoeff0(x, n, 0);
2028: tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
2029:
2030: case t_RFRAC: case t_RFRACN:
2031: w = precdl; av = avma;
2032: if (v<0) v = gvar(x);
2033: ex = ggval((GEN)x[2], polx[v]);
2034: precdl = n + ex + 1; x = gtoser(x, v); precdl = w;
2035: return gerepileupto(av, polcoeff0(x, n, v));
2036: }
2037: err(talker,"nonexistent component in truecoeff");
2038: return NULL; /* not reached */
2039: }
2040:
2041: GEN
2042: truecoeff(GEN x, long n)
2043: {
2044: return polcoeff0(x,n,-1);
2045: }
2046:
2047: GEN
2048: denom(GEN x)
2049: {
2050: long lx,i,av,tetpil;
2051: GEN s,t;
2052:
2053: switch(typ(x))
2054: {
2055: case t_INT: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
2056: return gun;
2057:
2058: case t_FRAC: case t_FRACN:
2059: return absi((GEN)x[2]);
2060:
2061: case t_COMPLEX:
2062: av=avma; t=denom((GEN)x[1]); s=denom((GEN)x[2]); tetpil=avma;
2063: return gerepile(av,tetpil,glcm(s,t));
2064:
2065: case t_QUAD:
2066: av=avma; t=denom((GEN)x[2]); s=denom((GEN)x[3]); tetpil=avma;
2067: return gerepile(av,tetpil,glcm(s,t));
2068:
2069: case t_POLMOD:
2070: return denom((GEN)x[2]);
2071:
2072: case t_RFRAC: case t_RFRACN:
2073: return gcopy((GEN)x[2]);
2074:
2075: case t_POL:
2076: return polun[varn(x)];
2077:
2078: case t_VEC: case t_COL: case t_MAT:
2079: lx=lg(x); if (lx==1) return gun;
2080: av = tetpil = avma; s = denom((GEN)x[1]);
2081: for (i=2; i<lx; i++)
2082: {
2083: t = denom((GEN)x[i]);
2084: /* t != gun est volontaire */
2085: if (t != gun) { tetpil=avma; s=glcm(s,t); }
2086: }
2087: return gerepile(av,tetpil,s);
2088: }
2089: err(typeer,"denom");
2090: return NULL; /* not reached */
2091: }
2092:
2093: GEN
2094: numer(GEN x)
2095: {
2096: long av,tetpil;
2097: GEN s;
2098:
2099: switch(typ(x))
2100: {
2101: case t_INT: case t_REAL: case t_INTMOD:
2102: case t_PADIC: case t_POL: case t_SER:
2103: return gcopy(x);
2104:
2105: case t_FRAC: case t_FRACN:
2106: return (signe(x[2])>0)? gcopy((GEN)x[1]): gneg((GEN)x[1]);
2107:
2108: case t_POLMOD:
2109: av=avma; s=numer((GEN)x[2]); tetpil=avma;
2110: return gerepile(av,tetpil,gmodulcp(s,(GEN)x[1]));
2111:
2112: case t_RFRAC: case t_RFRACN:
2113: return gcopy((GEN)x[1]);
2114:
2115: case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
2116: av=avma; s=denom(x); tetpil=avma;
2117: return gerepile(av,tetpil,gmul(s,x));
2118: }
2119: err(typeer,"numer");
2120: return NULL; /* not reached */
2121: }
2122:
2123: /* Lift only intmods if v does not occur in x, lift with respect to main
2124: * variable of x if v < 0, with respect to variable v otherwise.
2125: */
2126: GEN
2127: lift0(GEN x, long v)
2128: {
2129: long lx,tx=typ(x),i;
2130: GEN y;
2131:
2132: switch(tx)
2133: {
2134: case t_INT: case t_REAL:
2135: return gcopy(x);
2136:
2137: case t_INTMOD:
2138: return gcopy((GEN)x[2]);
2139:
2140: case t_POLMOD:
2141: if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
2142: y = cgetg(3,tx);
2143: y[1] = (long)lift0((GEN)x[1],v);
2144: y[2] = (long)lift0((GEN)x[2],v); return y;
2145:
2146: case t_SER:
2147: if (!signe(x)) return gcopy(x);
2148: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
2149: for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
2150: return y;
2151:
2152: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
2153: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
2154: lx=lg(x); y=cgetg(lx,tx);
2155: for (i=1; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
2156: return y;
2157:
2158: case t_POL:
2159: y=cgetg(lx=lgef(x),tx); y[1]=x[1];
2160: for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
2161: return y;
2162:
2163: case t_QUAD:
2164: y=cgetg(4,tx); copyifstack(x[1],y[1]);
2165: for (i=2; i<4; i++) y[i]=(long)lift0((GEN)x[i],v);
2166: return y;
2167: }
2168: err(typeer,"lift");
2169: return NULL; /* not reached */
2170: }
2171:
2172: GEN
2173: lift(GEN x)
2174: {
2175: return lift0(x,-1);
2176: }
2177:
2178: /* same as lift, without copy. May DESTROY x. For internal use only.
2179: Conventions on v as for lift. */
2180: GEN
2181: lift_intern0(GEN x, long v)
2182: {
2183: long i,lx,tx=typ(x);
2184:
2185: switch(tx)
2186: {
2187: case t_INT: case t_REAL:
2188: return x;
2189:
2190: case t_INTMOD:
2191: return (GEN)x[2];
2192:
2193: case t_POLMOD:
2194: if (v < 0 || v == varn((GEN)x[1])) return (GEN)x[2];
2195: x[1]=(long)lift_intern0((GEN)x[1],v);
2196: x[2]=(long)lift_intern0((GEN)x[2],v);
2197: return x;
2198:
2199: case t_SER: if (!signe(x)) return x; /* fall through */
2200: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_QUAD: case t_POL:
2201: case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
2202: lx = (tx==t_POL)? lgef(x): lg(x);
2203: for (i = lx-1; i>=lontyp[tx]; i--)
2204: x[i] = (long) lift_intern0((GEN)x[i],v);
2205: return x;
2206: }
2207: err(typeer,"lift_intern");
2208: return NULL; /* not reached */
2209: }
2210:
2211: /* memes conventions pour v que lift */
2212: GEN
2213: centerlift0(GEN x, long v)
2214: {
2215: long lx,tx=typ(x),i,av;
2216: GEN y;
2217:
2218: switch(tx)
2219: {
2220: case t_INT:
2221: return icopy(x);
2222:
2223: case t_INTMOD:
2224: av=avma; i=cmpii(shifti((GEN)x[2],1),(GEN)x[1]); avma=av;
2225: return (i>0)? subii((GEN)x[2],(GEN)x[1]): icopy((GEN)x[2]);
2226:
2227: case t_POLMOD:
2228: if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
2229: y=cgetg(3,tx);
2230: y[1]=(long)centerlift0((GEN)x[1],v); y[2]=(long)centerlift0((GEN)x[2],v);
2231: return y;
2232:
2233: case t_SER:
2234: if (!signe(x)) return gcopy(x);
2235: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
2236: for (i=2; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
2237: return y;
2238:
2239: case t_POL:
2240: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
2241: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
2242: lx = (tx==t_POL)? lgef(x): lg(x);
2243: y=cgetg(lx,tx); y[1]=x[1];
2244: for (i=lontyp[tx]; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
2245: return y;
2246:
2247: case t_QUAD:
2248: y=cgetg(4,tx); copyifstack(x[1],y[1]);
2249: for (i=2; i<4; i++) y[i]=(long)centerlift0((GEN)x[i],v);
2250: return y;
2251: }
2252: err(typeer,"centerlift");
2253: return NULL; /* not reached */
2254: }
2255:
2256: GEN
2257: centerlift(GEN x)
2258: {
2259: return centerlift0(x,-1);
2260: }
2261:
2262: /*******************************************************************/
2263: /* */
2264: /* PARTIES REELLE ET IMAGINAIRES */
2265: /* */
2266: /*******************************************************************/
2267:
2268: static GEN
2269: op_ReIm(GEN f(GEN), GEN x)
2270: {
2271: long lx,i,j,av,tetpil, tx = typ(x);
2272: GEN p1,p2,r2,i2,z;
2273:
2274: switch(tx)
2275: {
2276: case t_POL:
2277: lx=lgef(x); av=avma;
2278: for (i=lx-1; i>=2; i--)
2279: if (!gcmp0(f((GEN)x[i]))) break;
2280: avma=av; if (i==1) return zeropol(varn(x));
2281:
2282: z=cgetg(i+1,tx); z[1]=evalsigne(1)|evallgef(1+i)|evalvarn(varn(x));
2283: for (j=2; j<=i; j++) z[j] = (long)f((GEN)x[j]);
2284: return z;
2285:
2286: case t_SER:
2287: if (gcmp0(x)) { z=cgetg(2,tx); z[1]=x[1]; return z; }
2288: lx=lg(x); av=avma;
2289: for (i=2; i<lx; i++)
2290: if (!gcmp0(f((GEN)x[i]))) break;
2291: avma=av; if (i==lx) return zeroser(varn(x),lx-2+valp(x));
2292:
2293: z=cgetg(lx-i+2,tx); z[1]=x[1]; setvalp(z, valp(x)+i-2);
2294: for (j=2; i<lx; j++,i++) z[j] = (long) f((GEN)x[i]);
2295: return z;
2296:
2297: case t_RFRAC: case t_RFRACN:
2298: {
2299: av=avma; r2=greal((GEN)x[2]); i2=gimag((GEN)x[2]);
2300: if (f==greal)
2301: p1 = gadd(gmul(greal((GEN)x[1]),r2),
2302: gmul(gimag((GEN)x[1]),i2));
2303: else
2304: p1 = gsub(gmul(gimag((GEN)x[1]),r2),
2305: gmul(greal((GEN)x[1]),i2));
2306: p2=gadd(gsqr(r2), gsqr(i2));
2307: tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
2308: }
2309:
2310: case t_VEC: case t_COL: case t_MAT:
2311: lx=lg(x); z=cgetg(lx,tx);
2312: for (i=1; i<lx; i++) z[i] = (long) f((GEN)x[i]);
2313: return z;
2314: }
2315: err(typeer,"greal/gimag");
2316: return NULL; /* not reached */
2317: }
2318:
2319: GEN
2320: greal(GEN x)
2321: {
2322: switch(typ(x))
2323: {
2324: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
2325: return gcopy(x);
2326:
2327: case t_COMPLEX:
2328: return gcopy((GEN)x[1]);
2329:
2330: case t_QUAD:
2331: return gcopy((GEN)x[2]);
2332: }
2333: return op_ReIm(greal,x);
2334: }
2335:
2336: GEN
2337: gimag(GEN x)
2338: {
2339: switch(typ(x))
2340: {
2341: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
2342: return gzero;
2343:
2344: case t_COMPLEX:
2345: return gcopy((GEN)x[2]);
2346:
2347: case t_QUAD:
2348: return gcopy((GEN)x[3]);
2349: }
2350: return op_ReIm(gimag,x);
2351: }
2352:
2353: /*******************************************************************/
2354: /* */
2355: /* LOGICAL OPERATIONS */
2356: /* */
2357: /*******************************************************************/
2358:
2359: GEN
2360: glt(GEN x, GEN y) { return gcmp(x,y)<0? gun: gzero; }
2361:
2362: GEN
2363: gle(GEN x, GEN y) { return gcmp(x,y)<=0? gun: gzero; }
2364:
2365: GEN
2366: gge(GEN x, GEN y) { return gcmp(x,y)>=0? gun: gzero; }
2367:
2368: GEN
2369: ggt(GEN x, GEN y) { return gcmp(x,y)>0? gun: gzero; }
2370:
2371: GEN
2372: geq(GEN x, GEN y) { return gegal(x,y)? gun: gzero; }
2373:
2374: GEN
2375: gne(GEN x, GEN y) { return gegal(x,y)? gzero: gun; }
2376:
2377: GEN
2378: gand(GEN x, GEN y) { return gcmp0(x)? gzero: (gcmp0(y)? gzero: gun); }
2379:
2380: GEN
2381: gor(GEN x, GEN y) { return gcmp0(x)? (gcmp0(y)? gzero: gun): gun; }
2382:
2383: GEN
2384: gnot(GEN x) { return gcmp0(x)? gun: gzero; }
2385:
2386: /*******************************************************************/
2387: /* */
2388: /* FORMAL SIMPLIFICATIONS */
2389: /* */
2390: /*******************************************************************/
2391:
2392: GEN
2393: geval(GEN x)
2394: {
2395: long av,tetpil,lx,i, tx = typ(x);
2396: GEN y,z;
2397:
2398: if (is_const_t(tx)) return gcopy(x);
2399: if (is_graphicvec_t(tx) || tx == t_RFRACN)
2400: {
2401: lx=lg(x); y=cgetg(lx, tx);
2402: for (i=1; i<lx; i++) y[i] = (long)geval((GEN)x[i]);
2403: return y;
2404: }
2405:
2406: switch(tx)
2407: {
2408: case t_STR:
2409: return flisexpr(GSTR(x));
2410:
2411: case t_POLMOD: y=cgetg(3,tx);
2412: y[1]=(long)geval((GEN)x[1]);
2413: av=avma; z=geval((GEN)x[2]); tetpil=avma;
2414: y[2]=lpile(av,tetpil,gmod(z,(GEN)y[1]));
2415: return y;
2416:
2417: case t_POL:
2418: lx=lgef(x); if (lx==2) return gzero;
2419: {
2420: entree *ep = varentries[varn(x)];
2421: if (!ep) return gcopy(x);
2422: z = (GEN)ep->value;
2423: }
2424: y=gzero; av=avma;
2425: for (i=lx-1; i>1; i--)
2426: {
2427: tetpil=avma;
2428: y = gadd(geval((GEN)x[i]), gmul(z,y));
2429: }
2430: return gerepile(av,tetpil,y);
2431:
2432: case t_SER:
2433: err(impl, "evaluation of a power series");
2434:
2435: case t_RFRAC:
2436: return gdiv(geval((GEN)x[1]),geval((GEN)x[2]));
2437: }
2438: err(typeer,"geval");
2439: return NULL; /* not reached */
2440: }
2441:
2442: GEN
2443: simplify(GEN x)
2444: {
2445: long tx=typ(x),av,tetpil,i,lx;
2446: GEN p1,y;
2447:
2448: switch(tx)
2449: {
2450: case t_INT: case t_REAL: case t_FRAC:
2451: return is_universal_constant(x)? x: gcopy(x);
2452:
2453: case t_INTMOD: case t_PADIC: case t_QFR: case t_QFI:
2454: case t_LIST: case t_STR: case t_VECSMALL:
2455: return gcopy(x);
2456:
2457: case t_FRACN:
2458: return gred(x);
2459:
2460: case t_COMPLEX:
2461: if (is_universal_constant(x)) return x;
2462: if (isexactzero((GEN)x[2])) return simplify((GEN)x[1]);
2463: y=cgetg(3,tx);
2464: y[1]=(long)simplify((GEN)x[1]);
2465: y[2]=(long)simplify((GEN)x[2]); return y;
2466:
2467: case t_QUAD:
2468: if (isexactzero((GEN)x[3])) return simplify((GEN)x[2]);
2469: y=cgetg(4,tx);
2470: y[1]=lcopy((GEN)x[1]);
2471: y[2]=(long)simplify((GEN)x[2]);
2472: y[3]=(long)simplify((GEN)x[3]); return y;
2473:
2474: case t_POLMOD: y=cgetg(3,tx);
2475: y[1]=(long)simplify((GEN)x[1]);
2476: av=avma; p1=gmod((GEN)x[2],(GEN)y[1]); tetpil=avma;
2477: y[2]=lpile(av,tetpil,simplify(p1)); return y;
2478:
2479: case t_POL:
2480: lx=lgef(x); if (lx==2) return gzero;
2481: if (lx==3) return simplify((GEN)x[2]);
2482: y=cgetg(lx,tx); y[1]=x[1];
2483: for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
2484: return y;
2485:
2486: case t_SER:
2487: if (!signe(x)) return gcopy(x);
2488: lx=lg(x);
2489: y=cgetg(lx,tx); y[1]=x[1];
2490: for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
2491: return y;
2492:
2493: case t_RFRAC: y=cgetg(3,tx);
2494: y[1]=(long)simplify((GEN)x[1]);
2495: y[2]=(long)simplify((GEN)x[2]); return y;
2496:
2497: case t_RFRACN:
2498: av=avma; y=gred_rfrac(x); tetpil=avma;
2499: return gerepile(av,tetpil,simplify(y));
2500:
2501: case t_VEC: case t_COL: case t_MAT:
2502: lx=lg(x); y=cgetg(lx,tx);
2503: for (i=1; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
2504: return y;
2505: }
2506: err(typeer,"simplify");
2507: return NULL; /* not reached */
2508: }
2509:
2510: /*******************************************************************/
2511: /* */
2512: /* EVALUATION OF SOME SIMPLE FORMAL OBJECTS */
2513: /* */
2514: /*******************************************************************/
2515: static GEN
2516: qfeval0_i(GEN q, GEN x, long n)
2517: {
2518: long i,j,av=avma;
2519: GEN res=gzero;
2520:
2521: for (i=2;i<n;i++)
2522: for (j=1;j<i;j++)
2523: res = gadd(res, gmul(gcoeff(q,i,j), mulii((GEN)x[i],(GEN)x[j])) );
2524: res=gshift(res,1);
2525: for (i=1;i<n;i++)
2526: res = gadd(res, gmul(gcoeff(q,i,i), sqri((GEN)x[i])) );
2527: return gerepileupto(av,res);
2528: }
2529:
2530: #if 0
2531: static GEN
2532: qfeval0(GEN q, GEN x, long n)
2533: {
2534: long i,j,av=avma;
2535: GEN res=gzero;
2536:
2537: for (i=2;i<n;i++)
2538: for (j=1;j<i;j++)
2539: res = gadd(res, gmul(gcoeff(q,i,j), gmul((GEN)x[i],(GEN)x[j])) );
2540: res=gshift(res,1);
2541: for (i=1;i<n;i++)
2542: res = gadd(res, gmul(gcoeff(q,i,i), gsqr((GEN)x[i])) );
2543: return gerepileupto(av,res);
2544: }
2545: #else
2546: static GEN
2547: qfeval0(GEN q, GEN x, long n)
2548: {
2549: long i,j,av=avma;
2550: GEN res = gmul(gcoeff(q,1,1), gsqr((GEN)x[1]));
2551:
2552: for (i=2; i<n; i++)
2553: {
2554: GEN l = (GEN)q[i];
2555: GEN sx = gmul((GEN)l[1], (GEN)x[1]);
2556: for (j=2; j<i; j++)
2557: sx = gadd(sx, gmul((GEN)l[j],(GEN)x[j]));
2558: sx = gadd(gshift(sx,1), gmul((GEN)l[i],(GEN)x[i]));
2559:
2560: res = gadd(res, gmul((GEN)x[i], sx));
2561: }
2562: return gerepileupto(av,res);
2563: }
2564: #endif
2565:
2566: /* We assume q is a real symetric matrix */
2567: GEN
2568: qfeval(GEN q, GEN x)
2569: {
2570: long n=lg(q);
2571:
2572: if (n==1)
2573: {
2574: if (typ(q) != t_MAT || lg(x) != 1)
2575: err(talker,"invalid data in qfeval");
2576: return gzero;
2577: }
2578: if (typ(q) != t_MAT || lg(q[1]) != n)
2579: err(talker,"invalid quadratic form in qfeval");
2580: if (typ(x) != t_COL || lg(x) != n)
2581: err(talker,"invalid vector in qfeval");
2582:
2583: return qfeval0(q,x,n);
2584: }
2585:
2586: /* the Horner-type evaluation (mul x 2/3) would force us to use gmul and not
2587: * mulii (more than 4 x slower for small entries). Not worth it.
2588: */
2589: static GEN
2590: qfbeval0_i(GEN q, GEN x, GEN y, long n)
2591: {
2592: long i,j,av=avma;
2593: GEN res = gmul(gcoeff(q,1,1), mulii((GEN)x[1],(GEN)y[1]));
2594:
2595: for (i=2;i<n;i++)
2596: {
2597: for (j=1;j<i;j++)
2598: {
2599: GEN p1 = addii(mulii((GEN)x[i],(GEN)y[j]), mulii((GEN)x[j],(GEN)y[i]));
2600: res = gadd(res, gmul(gcoeff(q,i,j),p1));
2601: }
2602: res = gadd(res, gmul(gcoeff(q,i,i), mulii((GEN)x[i],(GEN)y[i])));
2603: }
2604: return gerepileupto(av,res);
2605: }
2606:
2607: #if 0
2608: static GEN
2609: qfbeval0(GEN q, GEN x, GEN y, long n)
2610: {
2611: long i,j,av=avma;
2612: GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1],(GEN)y[1]));
2613:
2614: for (i=2;i<n;i++)
2615: {
2616: for (j=1;j<i;j++)
2617: {
2618: GEN p1 = gadd(gmul((GEN)x[i],(GEN)y[j]), gmul((GEN)x[j],(GEN)y[i]));
2619: res = gadd(res, gmul(gcoeff(q,i,j),p1));
2620: }
2621: res = gadd(res, gmul(gcoeff(q,i,i), gmul((GEN)x[i],(GEN)y[i])));
2622: }
2623: return gerepileupto(av,res);
2624: }
2625: #else
2626: static GEN
2627: qfbeval0(GEN q, GEN x, GEN y, long n)
2628: {
2629: long i,j,av=avma;
2630: GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1], (GEN)y[1]));
2631:
2632: for (i=2; i<n; i++)
2633: {
2634: GEN l = (GEN)q[i];
2635: GEN sx = gmul((GEN)l[1], (GEN)y[1]);
2636: GEN sy = gmul((GEN)l[1], (GEN)x[1]);
2637: for (j=2; j<i; j++)
2638: {
2639: sx = gadd(sx, gmul((GEN)l[j],(GEN)y[j]));
2640: sy = gadd(sy, gmul((GEN)l[j],(GEN)x[j]));
2641: }
2642: sx = gadd(sx, gmul((GEN)l[i],(GEN)y[i]));
2643:
2644: sx = gmul((GEN)x[i], sx);
2645: sy = gmul((GEN)y[i], sy);
2646: res = gadd(res, gadd(sx,sy));
2647: }
2648: return gerepileupto(av,res);
2649: }
2650: #endif
2651:
2652: /* We assume q is a real symetric matrix */
2653: GEN
2654: qfbeval(GEN q, GEN x, GEN y)
2655: {
2656: long n=lg(q);
2657:
2658: if (n==1)
2659: {
2660: if (typ(q) != t_MAT || lg(x) != 1 || lg(y) != 1)
2661: err(talker,"invalid data in qfbeval");
2662: return gzero;
2663: }
2664: if (typ(q) != t_MAT || lg(q[1]) != n)
2665: err(talker,"invalid quadratic form in qfbeval");
2666: if (typ(x) != t_COL || lg(x) != n || typ(y) != t_COL || lg(y) != n)
2667: err(talker,"invalid vector in qfbeval");
2668:
2669: return qfbeval0(q,x,y,n);
2670: }
2671:
2672: /* yield X = M'.q.M, assuming q is symetric.
2673: * X_ij are X_ji identical, not copies
2674: * if flag is set, M has integer entries
2675: */
2676: GEN
2677: qf_base_change(GEN q, GEN M, int flag)
2678: {
2679: long i,j, k = lg(M), n=lg(q);
2680: GEN res = cgetg(k,t_MAT);
2681: GEN (*qf)(GEN,GEN,long) = flag? qfeval0_i: qfeval0;
2682: GEN (*qfb)(GEN,GEN,GEN,long) = flag? qfbeval0_i: qfbeval0;
2683:
2684: if (n==1)
2685: {
2686: if (typ(q) != t_MAT || k != 1)
2687: err(talker,"invalid data in qf_base_change");
2688: return res;
2689: }
2690: if (typ(M) != t_MAT || k == 1 || lg(M[1]) != n)
2691: err(talker,"invalid base change matrix in qf_base_change");
2692:
2693: for (i=1;i<k;i++)
2694: {
2695: res[i] = lgetg(k,t_COL);
2696: coeff(res,i,i) = (long) qf(q,(GEN)M[i],n);
2697: }
2698: for (i=2;i<k;i++)
2699: for (j=1;j<i;j++)
2700: coeff(res,i,j)=coeff(res,j,i) = (long) qfb(q,(GEN)M[i],(GEN)M[j],n);
2701: return res;
2702: }
2703:
2704: /* compute M'.M */
2705: GEN
2706: gram_matrix(GEN M)
2707: {
2708: long n=lg(M),i,j,k, av;
2709: GEN res = cgetg(n,t_MAT),p1;
2710:
2711: if (n==1)
2712: {
2713: if (typ(M) != t_MAT)
2714: err(talker,"invalid data in gram_matrix");
2715: return res;
2716: }
2717: if (typ(M) != t_MAT || lg(M[1]) != n)
2718: err(talker,"not a square matrix in gram_matrix");
2719:
2720: for (i=1;i<n;i++) res[i]=lgetg(n,t_COL);
2721: av=avma;
2722: for (i=1;i<n;i++)
2723: {
2724: p1 = gzero;
2725: for (k=1;k<n;k++)
2726: p1 = gadd(p1, gsqr(gcoeff(M,k,i)));
2727: coeff(res,i,i) = (long) gerepileupto(av,p1);
2728: av=avma;
2729: }
2730: for (i=2;i<n;i++)
2731: for (j=1;j<i;j++)
2732: {
2733: p1=gzero;
2734: for (k=1;k<n;k++)
2735: p1 = gadd(p1, gmul(gcoeff(M,k,i),gcoeff(M,k,j)));
2736: coeff(res,i,j)=coeff(res,j,i) = lpileupto(av,p1);
2737: av=avma;
2738: }
2739: return res;
2740: }
2741:
2742: /* return Re(x * y), x and y scalars */
2743: static GEN
2744: mul_real(GEN x, GEN y)
2745: {
2746: if (typ(x) == t_COMPLEX)
2747: {
2748: if (typ(y) == t_COMPLEX)
2749: {
2750: long av=avma, tetpil;
2751: GEN p1 = gmul((GEN)x[1], (GEN) y[1]);
2752: GEN p2 = gneg(gmul((GEN)x[2], (GEN) y[2]));
2753: tetpil=avma; return gerepile(av,tetpil,gadd(p1,p2));
2754: }
2755: x = (GEN)x[1];
2756: }
2757: else if (typ(y) == t_COMPLEX) y = (GEN)y[1];
2758: return gmul(x,y);
2759: }
2760:
2761: /* Compute Re(x * y), x and y matrices of compatible dimensions
2762: * assume lx, ly > 1, and scalar entries */
2763: GEN
2764: mulmat_real(GEN x, GEN y)
2765: {
2766: long av,tetpil,i,j,k, lx = lg(x), ly = lg(y), l = lg(x[1]);
2767: GEN p1,p2, z = cgetg(ly,t_MAT);
2768:
2769: for (j=1; j<ly; j++)
2770: {
2771: z[j] = lgetg(l,t_COL);
2772: for (i=1; i<l; i++)
2773: {
2774: p1=gzero; av=avma;
2775: for (k=1; k<lx; k++)
2776: {
2777: p2=mul_real(gcoeff(x,i,k),gcoeff(y,k,j));
2778: tetpil=avma; p1=gadd(p1,p2);
2779: }
2780: coeff(z,i,j)=lpile(av,tetpil,p1);
2781: }
2782: }
2783: return z;
2784: }
2785:
2786: static GEN
2787: hqfeval0(GEN q, GEN x, long n)
2788: {
2789: long i,j, av=avma;
2790: GEN res=gzero;
2791:
2792: for (i=2;i<n;i++)
2793: for (j=1;j<i;j++)
2794: {
2795: GEN p1 = gmul((GEN)x[i],gconj((GEN)x[j]));
2796: res = gadd(res, mul_real(gcoeff(q,i,j),p1));
2797: }
2798: res=gshift(res,1);
2799: for (i=1;i<n;i++)
2800: res = gadd(res, gmul(gcoeff(q,i,i), gnorm((GEN)x[i])) );
2801: return gerepileupto(av,res);
2802: }
2803:
2804: /* We assume q is a hermitian complex matrix */
2805: GEN
2806: hqfeval(GEN q, GEN x)
2807: {
2808: long n=lg(q);
2809:
2810: if (n==1)
2811: {
2812: if (typ(q) != t_MAT || lg(x) != 1)
2813: err(talker,"invalid data in hqfeval");
2814: return gzero;
2815: }
2816: if (typ(q) != t_MAT || lg(q[1]) != n)
2817: err(talker,"invalid quadratic form in hqfeval");
2818: if (typ(x) != t_COL || lg(x) != n)
2819: err(talker,"invalid vector in hqfeval");
2820:
2821: return hqfeval0(q,x,n);
2822: }
2823:
2824: GEN
2825: poleval(GEN x, GEN y)
2826: {
2827: long av,tetpil,i,j,imin,tx=typ(x);
2828: GEN p1,p2,p3,r,s;
2829:
2830: if (is_scalar_t(tx)) return gcopy(x);
2831: switch(tx)
2832: {
2833: case t_POL:
2834: i=lgef(x)-1; imin=2; break;
2835:
2836: case t_RFRAC: case t_RFRACN: av=avma;
2837: p1=poleval((GEN)x[1],y);
2838: p2=poleval((GEN)x[2],y); tetpil=avma;
2839: return gerepile(av,tetpil,gdiv(p1,p2));
2840:
2841: case t_VEC: case t_COL:
2842: i=lg(x)-1; imin=1; break;
2843: default: err(typeer,"poleval");
2844: }
2845: if (i<=imin)
2846: return (i==imin)? gcopy((GEN)x[imin]): gzero;
2847:
2848: av=avma; p1=(GEN)x[i]; i--;
2849: if (typ(y)!=t_COMPLEX)
2850: {
2851: #if 0 /* standard Horner's rule */
2852: for ( ; i>=imin; i--)
2853: p1 = gadd(gmul(p1,y),(GEN)x[i]);
2854: #endif
2855: /* specific attention to sparse polynomials */
2856: for ( ; i>=imin; i=j-1)
2857: {
2858: for (j=i; gcmp0((GEN)x[j]); j--)
2859: if (j==imin)
2860: {
2861: if (i!=j) y = gpuigs(y,i-j+1);
2862: tetpil=avma; return gerepile(av,tetpil,gmul(p1,y));
2863: }
2864: r = (i==j)? y: gpuigs(y,i-j+1);
2865: p1 = gadd(gmul(p1,r), (GEN)x[j]);
2866: }
2867: return gerepileupto(av,p1);
2868: }
2869:
2870: p2=(GEN)x[i]; i--; r=gtrace(y); s=gneg_i(gnorm(y));
2871: for ( ; i>=imin; i--)
2872: {
2873: p3=gadd(p2,gmul(r,p1));
2874: p2=gadd((GEN)x[i],gmul(s,p1)); p1=p3;
2875: }
2876: p1=gmul(y,p1); tetpil=avma;
2877: return gerepile(av,tetpil,gadd(p1,p2));
2878: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen3.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