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