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