![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen2.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|
Return to gen2.c CVS log | Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath |
1.1 noro 1: /* $Id: gen2.c,v 1.33 2001/10/01 12:11:30 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: /** (second part) **/
20: /** **/
21: /********************************************************************/
22: #include "pari.h"
23:
24: /*******************************************************************/
25: /* */
26: /* OPERATIONS BY VALUE */
27: /* f is a pointer to the function called. */
28: /* result is gaffect-ed to last parameter */
29: /* */
30: /*******************************************************************/
31:
32: void
33: gop1z(GEN (*f)(GEN), GEN x, GEN y)
34: {
35: long av=avma; gaffect(f(x),y); avma=av;
36: }
37:
38: void
39: gop2z(GEN (*f)(GEN, GEN), GEN x, GEN y, GEN z)
40: {
41: long av=avma; gaffect(f(x,y),z); avma=av;
42: }
43:
44: void
45: gops2gsz(GEN (*f)(GEN, long), GEN x, long s, GEN z)
46: {
47: long av=avma; gaffect(f(x,s),z); avma=av;
48: }
49:
50: void
51: gops2sgz(GEN (*f)(long, GEN), long s, GEN y, GEN z)
52: {
53: long av=avma; gaffect(f(s,y),z); avma=av;
54: }
55:
56: void
57: gops2ssz(GEN (*f)(long, long), long s, long y, GEN z)
58: {
59: long av=avma; gaffect(f(s,y),z); avma=av;
60: }
61:
62: /*******************************************************************/
63: /* */
64: /* OPERATIONS USING SMALL INTEGERS */
65: /* */
66: /*******************************************************************/
67:
68: /* small int prototype */
69: static long court_p[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
70:
71: GEN
72: gopsg2(GEN (*f)(GEN, GEN), long s, GEN y)
73: {
74: affsi(s,court_p); return f(court_p,y);
75: }
76:
77: GEN
78: gopgs2(GEN (*f)(GEN, GEN), GEN y, long s)
79: {
80: affsi(s,court_p); return f(y,court_p);
81: }
82:
83: long
84: opgs2(int (*f)(GEN, GEN), GEN y, long s)
85: {
86: affsi(s,court_p); return f(y,court_p);
87: }
88:
89: void
90: gopsg2z(GEN (*f)(GEN, GEN), long s, GEN y, GEN z)
91: {
92: long av=avma;
93: affsi(s,court_p); gaffect(f(court_p,y),z); avma=av;
94: }
95:
96: void
97: gopgs2z(GEN (*f)(GEN, GEN), GEN y, long s, GEN z)
98: {
99: long av=avma;
100: affsi(s,court_p); gaffect(f(y,court_p),z); avma=av;
101: }
102:
103: /*******************************************************************/
104: /* */
105: /* CREATION OF A P-ADIC GEN */
106: /* */
107: /*******************************************************************/
108:
109: /* y[4] not filled */
110: static GEN
111: cgetp2(GEN x, long v)
112: {
113: GEN y = cgetg(5,t_PADIC);
114: y[1] = evalprecp(precp(x)) | evalvalp(v);
115: icopyifstack(x[2], y[2]);
116: y[3] = licopy((GEN)x[3]); return y;
117: }
118:
119: GEN
120: cgetp(GEN x)
121: {
122: GEN y = cgetp2(x, 0);
123: y[4] = lgeti(lgefint(x[3])); return y;
124: }
125:
126: GEN
127: pureimag(GEN x)
128: {
129: GEN y = cgetg(3,t_COMPLEX);
130: y[1] = zero;
131: y[2] = (long)x; return y;
132: }
133:
134: /*******************************************************************/
135: /* */
136: /* SIZES */
137: /* */
138: /*******************************************************************/
139:
140: long
141: glength(GEN x)
142: {
143: switch(typ(x))
144: {
145: case t_INT: return lgefint(x)-2;
146: case t_POL: case t_LIST: return lgef(x)-2;
147: case t_REAL: return signe(x)? lg(x)-2: 0;
148: case t_STR: return strlen(GSTR(x));
149: }
150: return lg(x)-lontyp[typ(x)];
151: }
152:
153: GEN
154: matsize(GEN x)
155: {
156: GEN y=cgetg(3,t_VEC);
157:
158: switch(typ(x))
159: {
160: case t_VEC:
161: y[1]=un; y[2]=lstoi(lg(x)-1); break;
162: case t_COL:
163: y[1]=lstoi(lg(x)-1); y[2]=un; break;
164: case t_MAT:
165: y[2]=lstoi(lg(x)-1);
166: y[1]=(lg(x)==1)? zero: lstoi(lg(x[1])-1); break;
167: default: err(typeer,"matsize");
168: }
169: return y;
170: }
171:
172: long
173: taille(GEN x)
174: {
175: long i,n,lx, tx = typ(x);
176: n = lg(x);
177: if (is_recursive_t(tx))
178: {
179: lx = (tx==t_POL || tx==t_LIST)? lgef(x): n;
180: for (i=lontyp[tx]; i<lx; i++) n += taille((GEN)x[i]);
181: }
182: return n;
183: }
184:
185: long
186: taille2(GEN x) { return taille(x)<<TWOPOTBYTES_IN_LONG; }
187:
188: /*******************************************************************/
189: /* */
190: /* GREFFE */
191: /* Greffe d'une serie sur un polynome */
192: /* */
193: /*******************************************************************/
194:
195: GEN
196: greffe(GEN x, long l, long use_stack)
197: {
198: long i,e,k,vx;
199: GEN y;
200:
201: if (typ(x)!=t_POL) err(notpoler,"greffe");
202: if (use_stack) y = cgetg(l,t_SER);
203: else
204: {
205: y = (GEN) gpmalloc(l*sizeof(long));
206: y[0] = evaltyp(t_SER) | evallg(l);
207: }
208: if (gcmp0(x))
209: {
210: y[1]=evalvalp(l-2) | evalvarn(varn(x));
211: i=2; while(i<l) { y[i]=x[2]; i++; }
212: return y;
213: }
214: vx=varn(x); e=gval(x,vx);
215: y[1]=evalsigne(1) | evalvalp(e) | evalvarn(vx);
216: k=lgef(x)-e-1; i=l-1;
217: if (k<l)
218: while (i>k) { y[i]=zero; i--; }
219: while (i>=2) { y[i]=x[i+e]; i--; }
220: return y;
221: }
222:
223: /*******************************************************************/
224: /* */
225: /* CONVERSION GEN --> long */
226: /* */
227: /*******************************************************************/
228:
229: long
230: gtolong(GEN x)
231: {
232: long y,tx=typ(x),av=avma;
233:
234: switch(tx)
235: {
236: case t_INT:
237: return itos(x);
238: case t_REAL: case t_FRAC: case t_FRACN:
239: y=itos(ground(x)); avma=av; return y;
240: case t_COMPLEX:
241: if (gcmp0((GEN)x[2])) return gtolong((GEN)x[1]); break;
242: case t_QUAD:
243: if (gcmp0((GEN)x[3])) return gtolong((GEN)x[2]); break;
244: }
245: err(typeer,"gtolong");
246: return 0; /* not reached */
247: }
248:
249: /*******************************************************************/
250: /* */
251: /* COMPARISONS */
252: /* */
253: /*******************************************************************/
254:
255: /* returns 1 whenever x = 0, and 0 otherwise */
256: int
257: gcmp0(GEN x)
258: {
259: switch(typ(x))
260: {
261: case t_INT: case t_REAL: case t_POL: case t_SER:
262: return !signe(x);
263:
264: case t_INTMOD: case t_POLMOD:
265: return gcmp0((GEN)x[2]);
266:
267: case t_FRAC: case t_FRACN:
268: return !signe(x[1]);
269:
270: case t_COMPLEX:
271: /* is 0 iff norm(x) would be 0 (can happen with Re(x) and Im(x) != 0
272: * only if Re(x) and Im(x) are of type t_REAL). See mp.c:addrr().
273: */
274: if (gcmp0((GEN)x[1]))
275: {
276: if (gcmp0((GEN)x[2])) return 1;
277: if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
278: return (expo(x[1])>expo(x[2]));
279: }
280: if (gcmp0((GEN)x[2]))
281: {
282: if (typ(x[1])!=t_REAL || typ(x[2])!=t_REAL) return 0;
283: return (expo(x[2])>expo(x[1]));
284: }
285: return 0;
286:
287: case t_PADIC:
288: return !signe(x[4]);
289:
290: case t_QUAD:
291: return gcmp0((GEN)x[2]) && gcmp0((GEN)x[3]);
292:
293: case t_RFRAC: case t_RFRACN:
294: return gcmp0((GEN)x[1]);
295:
296: case t_VEC: case t_COL: case t_MAT:
297: {
298: long i;
299: for (i=lg(x)-1; i; i--)
300: if (!gcmp0((GEN)x[i])) return 0;
301: return 1;
302: }
303: }
304: return 0;
305: }
306:
307: /* returns 1 whenever x = 1, 0 otherwise */
308: int
309: gcmp1(GEN x)
310: {
311: switch(typ(x))
312: {
313: case t_INT:
314: return is_pm1(x) && signe(x)==1;
315:
316: case t_REAL:
317: if (signe(x) > 0 && expo(x)==0 && (ulong)x[2]==HIGHBIT)
318: {
319: long i,lx = lg(x);
320: for (i=3; i<lx; i++)
321: if (x[i]) return 0;
322: return 1;
323: }
324: return 0;
325:
326: case t_INTMOD: case t_POLMOD:
327: return gcmp1((GEN)x[2]);
328:
329: case t_FRAC:
330: return gcmp1((GEN)x[1]) && gcmp1((GEN)x[2]);
331:
332: case t_FRACN:
333: return egalii((GEN)x[1],(GEN)x[2]);
334:
335: case t_COMPLEX:
336: return gcmp1((GEN)x[1]) && gcmp0((GEN)x[2]);
337:
338: case t_PADIC:
339: return !valp(x) && gcmp1((GEN)x[4]);
340:
341: case t_QUAD:
342: return gcmp1((GEN)x[2]) && gcmp0((GEN)x[3]);
343:
344: case t_POL:
345: return lgef(x)==3 && gcmp1((GEN)x[2]);
346: }
347: return 0;
348: }
349:
350: /* returns 1 whenever the x = -1, 0 otherwise */
351: int
352: gcmp_1(GEN x)
353: {
354: long l,y;
355: GEN p1;
356:
357: switch(typ(x))
358: {
359: case t_INT:
360: return is_pm1(x) && signe(x)== -1;
361:
362: case t_REAL:
363: if (signe(x) < 0 && expo(x)==0 && (ulong)x[2]==HIGHBIT)
364: {
365: long i,lx = lg(x);
366: for (i=3; i<lx; i++)
367: if (x[i]) return 0;
368: return 1;
369: }
370: return 0;
371:
372: case t_INTMOD:
373: l=avma; y=egalii(addsi(1,(GEN)x[2]), (GEN)x[1]); avma=l; return y;
374:
375: case t_FRAC: case t_FRACN:
376: l = signe(x[1]);
377: return l && l == -signe(x[2]) && !absi_cmp((GEN)x[1],(GEN)x[2]);
378:
379: case t_COMPLEX:
380: return gcmp_1((GEN)x[1]) && gcmp0((GEN)x[2]);
381:
382: case t_QUAD:
383: return gcmp_1((GEN)x[2]) && gcmp0((GEN)x[3]);
384:
385: case t_PADIC:
386: l=avma; y=gegal(addsi(1,(GEN)x[4]), (GEN)x[3]); avma=l; return y;
387:
388: case t_POLMOD:
389: l=avma; p1=gadd(gun,(GEN)x[2]);
390: y = signe(p1) && !gegal(p1,(GEN)x[1]); avma=l; return !y;
391:
392: case t_POL:
393: return lgef(x)==3 && gcmp_1((GEN)x[2]);
394: }
395: return 0;
396: }
397:
398: /* returns the sign of x - y when it makes sense. 0 otherwise */
399: int
400: gcmp(GEN x, GEN y)
401: {
402: long tx,ty,f,av;
403:
404: tx=typ(x); ty=typ(y);
405: if (is_intreal_t(tx))
406: { if (is_intreal_t(ty)) return mpcmp(x,y); }
407: else
408: {
409: if (tx==t_STR)
410: {
411: if (ty != t_STR) return 1;
412: return strcmp(GSTR(x),GSTR(y));
413: }
414: if (!is_frac_t(tx)) err(typeer,"comparison");
415: }
416: if (ty == t_STR) return -1;
417: if (!is_intreal_t(ty) && !is_frac_t(ty)) err(typeer,"comparison");
418: av=avma; y=gneg_i(y); f=gsigne(gadd(x,y)); avma=av; return f;
419: }
420:
421: /* as gcmp for vector/matrices, using lexicographic ordering on components */
422: int
423: lexcmp(GEN x, GEN y)
424: {
425: const long tx=typ(x), ty=typ(y);
426: long lx,ly,l,fl,i;
427:
428: ly=lg(y);
429: if (!is_matvec_t(tx))
430: {
431: if (!is_matvec_t(ty)) return gcmp(x,y);
432: if (ly==1) return 1;
433: fl = lexcmp(x,(GEN)y[1]);
434: if (fl) return fl;
435: return (ly>2)? -1:0;
436: }
437:
438: lx=lg(x);
439: if (!is_matvec_t(ty))
440: {
441: if (lx==1) return -1;
442: fl = lexcmp(y,(GEN)x[1]);
443: if (fl) return -fl;
444: return (lx>2)? 1:0;
445: }
446:
447: /* x and y are matvec_t */
448: if (ly==1) return (lx==1)?0:1;
449: if (lx==1) return -1;
450: if (ty==t_MAT)
451: {
452: if (tx != t_MAT)
453: {
454: fl = lexcmp(x,(GEN)y[1]);
455: if (fl) return fl;
456: return (ly>2)?-1:0;
457: }
458: }
459: else if (tx==t_MAT)
460: {
461: fl = lexcmp(y,(GEN)x[1]);
462: if (fl) return -fl;
463: return (ly>2)?1:0;
464: }
465:
466: /* tx = ty = t_MAT, or x and y are both vect_t */
467: l=min(lx,ly);
468: for (i=1; i<l; i++)
469: {
470: fl = lexcmp((GEN)x[i],(GEN)y[i]);
471: if (fl) return fl;
472: }
473: return (ly != lx)? -1 : 0;
474: }
475:
476: /*****************************************************************/
477: /* */
478: /* EQUALITY */
479: /* returns 1 if x == y, 0 otherwise */
480: /* */
481: /*****************************************************************/
482:
483: /* x,y t_POL */
484: static int
485: polegal(GEN x, GEN y)
486: {
487: long i, lx = lgef(x);
488:
489: if (x[1] != y[1] && (lgef(y) != lx || lx > 3)) return 0;
490: for (i = 2; i < lx; i++)
491: if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
492: return 1;
493: }
494:
495: #define MASK(x) (((ulong)(x)) & (TYPBITS | LGBITS))
496: /* typ(x) = typ(y) = t_VEC/COL/MAT */
497: static int
498: vecegal(GEN x, GEN y)
499: {
500: long i;
501:
502: if (MASK(x[0]) != MASK(y[0])) return 0;
503:
504: i = lg(x)-1;
505: if (typ(x) != t_MAT)
506: {
507: for ( ; i; i--)
508: if (! gegal((GEN)x[i],(GEN)y[i]) ) return 0;
509: return 1;
510: }
511: for ( ; i; i--)
512: if (! vecegal((GEN)x[i],(GEN)y[i]) ) return 0;
513: return 1;
514: }
515:
516: int
517: gegal(GEN x, GEN y)
518: {
519: ulong av,tx;
520: long i;
521:
522: if (x == y) return 1;
523: tx = typ(x);
524: if (tx==typ(y))
525: switch(tx)
526: {
527: case t_INT:
528: return egalii(x,y);
529:
530: case t_POL:
531: return polegal(x,y);
532:
533: case t_COMPLEX:
534: return gegal((GEN)x[1],(GEN)y[1]) && gegal((GEN)x[2],(GEN)y[2]);
535:
536: case t_INTMOD: case t_POLMOD:
537: return gegal((GEN)x[2],(GEN)y[2])
538: && (x[1]==y[1] || gegal((GEN)x[1],(GEN)y[1]));
539:
540: case t_QFR:
541: if (!gegal((GEN)x[4],(GEN)y[4])) return 0; /* fall through */
542: case t_QUAD: case t_QFI:
543: return gegal((GEN)x[1],(GEN)y[1])
544: && gegal((GEN)x[2],(GEN)y[2])
545: && gegal((GEN)x[3],(GEN)y[3]);
546:
547: case t_FRAC:
548: return gegal((GEN)x[1], (GEN)y[1])
549: && gegal((GEN)x[2], (GEN)y[2]);
550:
551: case t_FRACN: case t_RFRAC: case t_RFRACN:
552: av=avma; i=gegal(gmul((GEN)x[1],(GEN)y[2]),gmul((GEN)x[2],(GEN)y[1]));
553: avma=av; return i;
554:
555: case t_STR:
556: return !strcmp(GSTR(x),GSTR(y));
557:
558: case t_VEC: case t_COL: case t_MAT:
559: return vecegal(x,y);
560: case t_VECSMALL:
561: if (MASK(x[0]) != MASK(y[0])) return 0;
562: for (i = lg(x)-1; i; i--)
563: if (x[i] != y[i]) return 0;
564: return 1;
565: }
566: {
567: jmp_buf env;
568: VOLATILE long av = avma;
569: void *c;
570: if (setjmp(env)) i = 0;
571: else
572: {
573: c = err_catch(-1, env, NULL);
574: i = gcmp0(gadd(x, gneg_i(y)));
575: }
576: err_leave(&c);
577: avma = av;
578: }
579: return i;
580: }
581: #undef MASK
582:
583: /*******************************************************************/
584: /* */
585: /* VALUATION */
586: /* p is either an int or a polynomial. */
587: /* returns the largest exponent of p dividing x when this makes */
588: /* sense : error for types real, integermod and polymod if p does */
589: /* not divide the modulus, q-adic if q!=p. */
590: /* */
591: /*******************************************************************/
592:
593: static long
594: minval(GEN x, GEN p, long first, long lx)
595: {
596: long i,k, val = VERYBIGINT;
597: for (i=first; i<lx; i++)
598: {
599: k = ggval((GEN)x[i],p);
600: if (k < val) val = k;
601: }
602: return val;
603: }
604:
605: long
606: ggval(GEN x, GEN p)
607: {
608: long tx=typ(x), tp=typ(p), av, limit,vx,v,i,val;
609: GEN p1,p2;
610:
611: if (isexactzero(x)) return VERYBIGINT;
612: if (is_const_t(tx) && tp==t_POL) return 0;
613: switch(tx)
614: {
615: case t_INT:
616: if (tp != t_INT) break;
617: av=avma; val = pvaluation(x,p, &p1);
618: avma=av; return val;
619:
620: case t_INTMOD:
621: av=avma;
622: p1=cgeti(lgefint(x[1]));
623: p2=cgeti(lgefint(x[2]));
624: if (tp!=t_INT || !mpdivis((GEN)x[1],p,p1)) break;
625: if (!mpdivis((GEN)x[2],p,p2)) { avma=av; return 0; }
626: val=1; while (mpdivis(p1,p,p1) && mpdivis(p2,p,p2)) val++;
627: avma=av; return val;
628:
629: case t_PADIC:
630: if (tp!=t_INT || !gegal(p,(GEN)x[2])) break;
631: return valp(x);
632:
633: case t_POLMOD:
634: if (tp==t_INT) return ggval((GEN)x[2],p);
635: av = avma;
636: if (tp!=t_POL || !poldivis((GEN)x[1],p,&p1)) break;
637: if (!poldivis((GEN)x[2],p,&p2)) { avma=av; return 0; }
638: val=1; while (poldivis(p1,p,&p1)&&poldivis(p2,p,&p2)) val++;
639: avma = av; return val;
640:
641: case t_POL:
642: if (tp==t_POL)
643: {
644: v=varn(p); vx=varn(x);
645: if (vx == v)
646: {
647: if ((p>=(GEN)polx && p <= (GEN)(polx+MAXVARN)) || ismonome(p))
648: {
649: i=2; while (isexactzero((GEN)x[i])) i++;
650: return i-2;
651: }
652: av = avma; limit=stack_lim(av,1);
653: for (val=0; ; val++)
654: {
655: if (!poldivis(x,p,&x)) break;
656: if (low_stack(limit, stack_lim(av,1)))
657: {
658: if(DEBUGMEM>1) err(warnmem,"ggval");
659: x = gerepilecopy(av, x);
660: }
661: }
662: avma = av; return val;
663: }
664: if (vx > v) return 0;
665: }
666: else if (tp!=t_INT) break;
667: i=2; while (isexactzero((GEN)x[i])) i++;
668: return minval(x,p,i,lgef(x));
669:
670: case t_SER:
671: if (tp!=t_POL && tp!=t_SER && tp!=t_INT) break;
672: v=gvar(p); vx=varn(x);
673: if (vx==v) return (long)(valp(x)/ggval(p,polx[v]));
674: if (vx>v) return 0;
675: return minval(x,p,2,lg(x));
676:
677: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
678: return ggval((GEN)x[1],p) - ggval((GEN)x[2],p);
679:
680: case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
681: return minval(x,p,1,lg(x));
682: }
683: err(talker,"forbidden or conflicting type in gval");
684: return 0; /* not reached */
685: }
686:
687: /* x is non-zero */
688: long
689: svaluation(ulong x, ulong p, long *py)
690: {
691: ulong v = 0;
692: for(;;)
693: {
694: if (x%p) { *py = x; return v; }
695: v++; x/=p;
696: }
697: }
698:
699: /* x is a non-zero integer */
700: long
701: pvaluation(GEN x, GEN p, GEN *py)
702: {
703: long av,v;
704: GEN p1,p2;
705:
706: if (egalii(p,gdeux))
707: {
708: v = vali(x);
709: if (py) *py = shifti(x, -v);
710: return v;
711: }
712: if (is_pm1(p))
713: {
714: v = (signe(p) < 0 && signe(x) < 0);
715: if (py) { *py = v? negi(x): icopy(x); }
716: return v;
717: }
718: if (!is_bigint(x))
719: {
720: long y;
721: if (!is_bigint(p))
722: {
723: v = svaluation(x[2],p[2], &y);
724: if (signe(x) < 0) y = -y;
725: if (py) *py = stoi(y);
726: }
727: else
728: {
729: v = 0;
730: if (py) *py = icopy(x);
731: }
732: return v;
733: }
734: av = avma; v = 0; (void)new_chunk(lgefint(x));
735: for(;;)
736: {
737: p1 = dvmdii(x,p,&p2);
738: if (p2 != gzero) { avma=av; if (py) *py = icopy(x); return v; }
739: v++; x = p1;
740: }
741: }
742:
743: /*******************************************************************/
744: /* */
745: /* NEGATION: Create -x */
746: /* */
747: /*******************************************************************/
748:
749: GEN
750: gneg(GEN x)
751: {
752: long tx=typ(x),lx,i;
753: GEN y;
754:
755: if (gcmp0(x)) return gcopy(x);
756: switch(tx)
757: {
758: case t_INT: case t_REAL:
759: return mpneg(x);
760:
761: case t_INTMOD: y=cgetg(3,t_INTMOD);
762: icopyifstack(x[1],y[1]);
763: y[2] = lsubii((GEN)y[1],(GEN)x[2]);
764: break;
765:
766: case t_POLMOD: y=cgetg(3,t_POLMOD);
767: copyifstack(x[1],y[1]);
768: y[2]=lneg((GEN)x[2]); break;
769:
770: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
771: y=cgetg(3,tx);
772: y[1]=lneg((GEN)x[1]);
773: y[2]=lcopy((GEN)x[2]); break;
774:
775: case t_PADIC:
776: y=cgetp2(x,valp(x));
777: y[4]=lsubii((GEN)x[3],(GEN)x[4]);
778: break;
779:
780: case t_QUAD:
781: y=cgetg(4,t_QUAD); copyifstack(x[1],y[1]);
782: y[2]=lneg((GEN)x[2]);
783: y[3]=lneg((GEN)x[3]); break;
784:
785: case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
786: lx=lg(x); y=cgetg(lx,tx);
787: for (i=1; i<lx; i++) y[i]=lneg((GEN)x[i]);
788: break;
789:
790: case t_POL:
791: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
792: for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
793: break;
794:
795: case t_SER:
796: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
797: for (i=2; i<lx; i++) y[i]=lneg((GEN)x[i]);
798: break;
799:
800: default:
801: err(typeer,"negation");
802: return NULL; /* not reached */
803: }
804: return y;
805: }
806:
807: GEN
808: gneg_i(GEN x)
809: {
810: long tx=typ(x),lx,i;
811: GEN y;
812:
813: if (gcmp0(x)) return x;
814: switch(tx)
815: {
816: case t_INT: case t_REAL:
817: return mpneg(x);
818:
819: case t_INTMOD: y=cgetg(3,t_INTMOD); y[1]=x[1];
820: y[2] = lsubii((GEN)y[1],(GEN)x[2]);
821: break;
822:
823: case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
824: y=cgetg(3,tx); y[2]=x[2];
825: y[1]=(long)gneg_i((GEN)x[1]); break;
826:
827: case t_PADIC: y = cgetg(5,t_PADIC); y[2]=x[2]; y[3]=x[3];
828: y[1] = evalprecp(precp(x)) | evalvalp(valp(x));
829: y[4]=lsubii((GEN)x[3],(GEN)x[4]); break;
830:
831: case t_POLMOD: y=cgetg(3,t_POLMOD); y[1]=x[1];
832: y[2]=(long)gneg_i((GEN)x[2]); break;
833:
834: case t_QUAD: y=cgetg(4,t_QUAD); y[1]=x[1];
835: y[2]=(long)gneg_i((GEN)x[2]);
836: y[3]=(long)gneg_i((GEN)x[3]); break;
837:
838: case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
839: lx=lg(x); y=cgetg(lx,tx);
840: for (i=1; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
841: break;
842:
843: case t_POL: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
844: for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
845: break;
846:
847: case t_SER: lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
848: for (i=2; i<lx; i++) y[i]=(long)gneg_i((GEN)x[i]);
849: break;
850:
851: default:
852: err(typeer,"negation");
853: return NULL; /* not reached */
854: }
855: return y;
856: }
857:
858: /******************************************************************/
859: /* */
860: /* ABSOLUTE VALUE */
861: /* Create abs(x) if x is integer, real, fraction or complex. */
862: /* Error otherwise. */
863: /* */
864: /******************************************************************/
865:
866: GEN
867: gabs(GEN x, long prec)
868: {
869: long tx=typ(x),lx,i,l,tetpil;
870: GEN y,p1;
871:
872: switch(tx)
873: {
874: case t_INT: case t_REAL:
875: return mpabs(x);
876:
877: case t_FRAC: case t_FRACN: y=cgetg(lg(x),tx);
878: y[1]=labsi((GEN)x[1]);
879: y[2]=labsi((GEN)x[2]); return y;
880:
881: case t_COMPLEX:
882: l=avma; p1=gnorm(x);
883: switch(typ(p1))
884: {
885: case t_INT:
886: if (!carrecomplet(p1, &y)) break;
887: return gerepileupto(l, y);
888: case t_FRAC:
889: case t_FRACN:
890: {
891: GEN a,b;
892: if (!carrecomplet((GEN)p1[1], &a)) break;
893: if (!carrecomplet((GEN)p1[2], &b)) break;
894: return gerepileupto(l, gdiv(a,b));
895: }
896: }
897: tetpil=avma;
898: return gerepile(l,tetpil,gsqrt(p1,prec));
899:
900: case t_QUAD:
901: l=avma; p1=gmul(x, realun(prec)); tetpil=avma;
902: return gerepile(l,tetpil,gabs(p1,prec));
903:
904: case t_POL:
905: lx=lgef(x); if (lx<=2) return gcopy(x);
906: p1=(GEN)x[lx-1];
907: switch(typ(p1))
908: {
909: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
910: if (gsigne(p1)<0) return gneg(x);
911: }
912: return gcopy(x);
913:
914: case t_SER:
915: if (valp(x) || !signe(x) || lg(x)<3)
916: err(talker,"abs is not analytic at 0");
917: if (gsigne((GEN) x[2])<0) return gneg(x);
918: return gcopy(x);
919:
920: case t_VEC: case t_COL: case t_MAT:
921: lx=lg(x); y=cgetg(lx,tx);
922: for (i=1; i<lx; i++)
923: y[i]=(long)gabs((GEN)x[i],prec);
924: return y;
925: }
926: err(typeer,"gabs");
927: return NULL; /* not reached */
928: }
929:
930: GEN
931: gmax(GEN x, GEN y)
932: {
933: if (gcmp(x,y)>0) y=x; return gcopy(y);
934: }
935:
936: GEN
937: gmin(GEN x, GEN y)
938: {
939: if (gcmp(x,y)<0) y=x; return gcopy(y);
940: }
941:
942: GEN
943: vecmax(GEN x)
944: {
945: long tx=typ(x),lx,lx2,i,j;
946: GEN *p1,s;
947:
948: if (!is_matvec_t(tx)) return gcopy(x);
949: lx=lg(x); if (lx==1) return stoi(-VERYBIGINT);
950: if (tx!=t_MAT)
951: {
952: s=(GEN)x[1];
953: for (i=2; i<lx; i++)
954: if (gcmp((GEN)x[i],s) > 0) s=(GEN)x[i];
955: }
956: else
957: {
958: lx2 = lg(x[1]);
959: if (lx2==1) return stoi(-VERYBIGINT);
960: s=gmael(x,1,1); i=2;
961: for (j=1; j<lx; j++)
962: {
963: p1 = (GEN *) x[j];
964: for (; i<lx2; i++)
965: if (gcmp(p1[i],s) > 0) s=p1[i];
966: i=1;
967: }
968: }
969: return gcopy(s);
970: }
971:
972: GEN
973: vecmin(GEN x)
974: {
975: long tx=typ(x),lx,lx2,i,j;
976: GEN *p1,s;
977:
978: if (!is_matvec_t(tx)) return gcopy(x);
979: lx=lg(x); if (lx==1) return stoi(VERYBIGINT);
980: if (tx!=t_MAT)
981: {
982: s=(GEN)x[1];
983: for (i=2; i<lx; i++)
984: if (gcmp((GEN)x[i],s) < 0) s=(GEN)x[i];
985: }
986: else
987: {
988: lx2 = lg(x[1]);
989: if (lx2==1) return stoi(VERYBIGINT);
990: s=gmael(x,1,1); i=2;
991: for (j=1; j<lx; j++)
992: {
993: p1 = (GEN *) x[j];
994: for (; i<lx2; i++)
995: if (gcmp(p1[i],s) < 0) s=p1[i];
996: i=1;
997: }
998: }
999: return gcopy(s);
1000: }
1001:
1002: /*******************************************************************/
1003: /* */
1004: /* AFFECT long --> GEN */
1005: /* affect long s to GEN x. Useful for initialization. */
1006: /* */
1007: /*******************************************************************/
1008:
1009: static void
1010: padicaff0(GEN x)
1011: {
1012: if (signe(x[4]))
1013: {
1014: setvalp(x, valp(x)|precp(x));
1015: affsi(0,(GEN)x[4]);
1016: }
1017: }
1018:
1019: void
1020: gaffsg(long s, GEN x)
1021: {
1022: long l,i,v;
1023: GEN p;
1024:
1025: switch(typ(x))
1026: {
1027: case t_INT:
1028: affsi(s,x); break;
1029:
1030: case t_REAL:
1031: affsr(s,x); break;
1032:
1033: case t_INTMOD:
1034: modsiz(s,(GEN)x[1],(GEN)x[2]); break;
1035:
1036: case t_FRAC: case t_FRACN:
1037: affsi(s,(GEN)x[1]); affsi(1,(GEN)x[2]); break;
1038:
1039: case t_COMPLEX:
1040: gaffsg(s,(GEN)x[1]); gaffsg(0,(GEN)x[2]); break;
1041:
1042: case t_PADIC:
1043: if (!s) { padicaff0(x); break; }
1044: p = (GEN)x[2];
1045: v = (is_bigint(p))? 0: svaluation(s,p[2], &i);
1046: setvalp(x,v); modsiz(i,(GEN)x[3],(GEN)x[4]);
1047: break;
1048:
1049: case t_QUAD:
1050: gaffsg(s,(GEN)x[2]); gaffsg(0,(GEN)x[3]); break;
1051:
1052: case t_POLMOD:
1053: gaffsg(s,(GEN)x[2]); break;
1054:
1055: case t_POL:
1056: v=varn(x);
1057: if (!s) x[1]=evallgef(2) | evalvarn(v);
1058: else
1059: {
1060: x[1]=evalsigne(1) | evallgef(3) | evalvarn(v);
1061: gaffsg(s,(GEN)x[2]);
1062: }
1063: break;
1064:
1065: case t_SER:
1066: v=varn(x); gaffsg(s,(GEN)x[2]); l=lg(x);
1067: if (!s)
1068: x[1] = evalvalp(l-2) | evalvarn(v);
1069: else
1070: x[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
1071: for (i=3; i<l; i++) gaffsg(0,(GEN)x[i]);
1072: break;
1073:
1074: case t_RFRAC: case t_RFRACN:
1075: gaffsg(s,(GEN)x[1]); gaffsg(1,(GEN)x[2]); break;
1076:
1077: case t_VEC: case t_COL: case t_MAT:
1078: if (lg(x)!=2) err(operi,"",t_INT,typ(x));
1079: gaffsg(s,(GEN)x[1]); break;
1080:
1081: default: err(operf,"",t_INT,typ(x));
1082: }
1083: }
1084:
1085: /*******************************************************************/
1086: /* */
1087: /* GENERIC AFFECTATION */
1088: /* Affect the content of x to y, whenever possible */
1089: /* */
1090: /*******************************************************************/
1091: void
1092: gaffect(GEN x, GEN y)
1093: {
1094: long i,j,k,l,v,vy,lx,ly,tx,ty,av;
1095: GEN p1,num,den;
1096:
1097:
1098: tx=typ(x); ty=typ(y); lx=lg(x); ly=lg(y);
1099: if (is_scalar_t(tx))
1100: {
1101: if (is_scalar_t(ty))
1102: {
1103: switch(tx)
1104: {
1105: case t_INT:
1106: switch(ty)
1107: {
1108: case t_INT:
1109: /* y = gzero, gnil, gun or gdeux */
1110: if (is_universal_constant(y))
1111: {
1112: if (y==gzero) err(overwriter,"gaffect (gzero)");
1113: if (y==gun) err(overwriter,"gaffect (gun)");
1114: if (y==gdeux) err(overwriter,"gaffect (gdeux)");
1115: err(overwriter,"gaffect (gnil)");
1116: }
1117: affii(x,y); break;
1118:
1119: case t_REAL:
1120: if (y == gpi) err(overwriter,"gaffect (gpi)");
1121: if (y==geuler) err(overwriter,"gaffect (geuler)");
1122: affir(x,y); break;
1123:
1124: case t_INTMOD:
1125: modiiz(x,(GEN)y[1],(GEN)y[2]); break;
1126:
1127: case t_FRAC:
1128: if (y == ghalf) err(overwriter,"gaffect (ghalf)");
1129: case t_FRACN:
1130: affii(x,(GEN)y[1]); affsi(1,(GEN)y[2]); break;
1131:
1132: case t_COMPLEX:
1133: if (y == gi) err(overwriter,"gaffect (gi)");
1134: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
1135:
1136: case t_PADIC:
1137: if (!signe(x)) { padicaff0(y); break; }
1138: l=avma;
1139: setvalp(y, pvaluation(x,(GEN)y[2],&p1));
1140: modiiz(p1,(GEN)y[3],(GEN)y[4]);
1141: avma=l; break;
1142:
1143: case t_QUAD:
1144: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
1145:
1146: case t_POLMOD:
1147: gaffect(x,(GEN)y[2]); break;
1148:
1149: default: err(operf,"",tx,ty);
1150: }
1151: break;
1152:
1153: case t_REAL:
1154: switch(ty)
1155: {
1156: case t_REAL:
1157: affrr(x,y); break;
1158:
1159: case t_COMPLEX:
1160: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
1161:
1162: case t_POLMOD:
1163: gaffect(x,(GEN)y[2]); break;
1164:
1165: case t_INT: case t_INTMOD: case t_FRAC:
1166: case t_FRACN: case t_PADIC: case t_QUAD:
1167: err(operf,"",tx,ty);
1168:
1169: default: err(operf,"",tx,ty);
1170: }
1171: break;
1172:
1173: case t_INTMOD:
1174: switch(ty)
1175: {
1176: case t_INTMOD:
1177: if (!divise((GEN)x[1],(GEN)y[1]))
1178: err(operi,"",tx,ty);
1179: modiiz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;
1180:
1181: case t_POLMOD:
1182: gaffect(x,(GEN)y[2]); break;
1183:
1184: default: err(operf,"",tx,ty);
1185: }
1186: break;
1187:
1188: case t_FRAC:
1189: switch(ty)
1190: {
1191: case t_INT:
1192: if (! mpdivis((GEN)x[1],(GEN)x[2],y))
1193: err(operi,"",tx,ty);
1194: break;
1195:
1196: case t_REAL:
1197: diviiz((GEN)x[1],(GEN)x[2],y); break;
1198:
1199: case t_INTMOD:
1200: av=avma; p1=mpinvmod((GEN)x[2],(GEN)y[1]);
1201: modiiz(mulii((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
1202: avma=av; break;
1203:
1204: case t_FRAC: case t_FRACN:
1205: affii((GEN)x[1],(GEN)y[1]);
1206: affii((GEN)x[2],(GEN)y[2]);
1207: break;
1208:
1209: case t_COMPLEX:
1210: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
1211:
1212: case t_PADIC:
1213: if (!signe(x[1])) { padicaff0(y); break; }
1214: av=avma; num=(GEN)x[1]; den=(GEN)x[2];
1215: v = pvaluation(num, (GEN) y[2], &num);
1216: if (!v) v = -pvaluation(den,(GEN)y[2],&den);
1217: p1=mulii(num,mpinvmod(den,(GEN)y[3]));
1218: modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
1219: setvalp(y,v); break;
1220:
1221: case t_QUAD:
1222: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
1223: case t_POLMOD:
1224: gaffect(x,(GEN)y[2]); break;
1225: default: err(operf,"",tx,ty);
1226: }
1227: break;
1228:
1229: case t_FRACN:
1230: switch(ty)
1231: {
1232: case t_INT:
1233: if (! mpdivis((GEN)x[1],(GEN)x[2],y)) err(operi,"",tx,ty);
1234: break;
1235:
1236: case t_REAL:
1237: diviiz((GEN)x[1],(GEN)x[2],y); break;
1238: case t_INTMOD:
1239: av=avma; gaffect(gred(x),y); avma=av; break;
1240: case t_FRAC:
1241: gredz(x,y); break;
1242: case t_FRACN:
1243: affii((GEN)x[1],(GEN)y[1]);
1244: affii((GEN)x[2],(GEN)y[2]); break;
1245: case t_COMPLEX:
1246: gaffect(x,(GEN)y[1]); gaffsg(0,(GEN)y[2]); break;
1247: case t_PADIC:
1248: if (!signe(x[1])) { padicaff0(y); break; }
1249: av=avma; num=(GEN)x[1]; den=(GEN)x[2];
1250: v = pvaluation(num,(GEN)y[2],&num)
1251: - pvaluation(den,(GEN)y[2],&den);
1252: p1=mulii(num,mpinvmod(den,(GEN)y[3]));
1253: modiiz(p1,(GEN)y[3],(GEN)y[4]); avma=av;
1254: setvalp(y,v); break;
1255:
1256: case t_QUAD:
1257: gaffect(x,(GEN)y[2]); gaffsg(0,(GEN)y[3]); break;
1258: case t_POLMOD:
1259: gaffect(x,(GEN)y[2]); break;
1260: default: err(operf,"",tx,ty);
1261: }
1262: break;
1263:
1264: case t_COMPLEX:
1265: switch(ty)
1266: {
1267: case t_INT: case t_REAL: case t_INTMOD:
1268: case t_FRAC: case t_FRACN: case t_PADIC: case t_QUAD:
1269: if (!gcmp0((GEN)x[2])) err(operi,"",tx,ty);
1270: gaffect((GEN)x[1],y); break;
1271:
1272: case t_COMPLEX:
1273: gaffect((GEN)x[1],(GEN)y[1]);
1274: gaffect((GEN)x[2],(GEN)y[2]);
1275: break;
1276:
1277: case t_POLMOD:
1278: gaffect(x,(GEN)y[2]); break;
1279:
1280: default: err(operf,"",tx,ty);
1281: }
1282: break;
1283:
1284: case t_PADIC:
1285: switch(ty)
1286: {
1287: case t_INTMOD:
1288: if (valp(x)<0) err(operi,"",tx,ty);
1289: av=avma;
1290: v = pvaluation((GEN)y[1],(GEN)x[2],&p1);
1291: k = signe(x[4])? (precp(x)+valp(x)): valp(x)+1;
1292: if (!gcmp1(p1) || v > k) err(operi,"",tx,ty);
1293: modiiz(gtrunc(x),(GEN)y[1],(GEN)y[2]); avma=av; break;
1294:
1295: case t_PADIC:
1296: if (!egalii((GEN)x[2],(GEN)y[2])) err(operi,"",tx,ty);
1297: modiiz((GEN)x[4],(GEN)y[3],(GEN)y[4]);
1298: setvalp(y,valp(x)); break;
1299:
1300: case t_POLMOD:
1301: gaffect(x,(GEN)y[2]); break;
1302:
1303: default: err(operf,"",tx,ty);
1304: }
1305: break;
1306:
1307: case t_QUAD:
1308: switch(ty)
1309: {
1310: case t_INT: case t_INTMOD: case t_FRAC:
1311: case t_FRACN: case t_PADIC:
1312: if (!gcmp0((GEN)x[3])) err(operi,"",tx,ty);
1313: gaffect((GEN)x[2],y); break;
1314:
1315: case t_REAL:
1316: av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av; break;
1317: case t_COMPLEX:
1318: ly=precision(y);
1319: if (ly)
1320: {
1321: av=avma; p1=co8(x,ly); gaffect(p1,y); avma=av;
1322: }
1323: else
1324: {
1325: if (!gcmp0((GEN)x[3])) err(operi,"",tx,ty);
1326: gaffect((GEN)x[2],y);
1327: }
1328: break;
1329:
1330: case t_QUAD:
1331: if (! gegal((GEN)x[1],(GEN)y[1])) err(operi,"",tx,ty);
1332: affii((GEN)x[2],(GEN)y[2]);
1333: affii((GEN)x[3],(GEN)y[3]);
1334: break;
1335: case t_POLMOD:
1336: gaffect(x,(GEN)y[2]); break;
1337: default: err(operf,"",tx,ty);
1338: }
1339: break;
1340:
1341: case t_POLMOD:
1342: if (ty!=t_POLMOD) err(operf,"",tx,ty);
1343: if (! gdivise((GEN)x[1],(GEN)y[1])) err(operi,"",tx,ty);
1344: gmodz((GEN)x[2],(GEN)y[1],(GEN)y[2]); break;
1345:
1346: default: err(operf,"",tx,ty);
1347: }
1348: return;
1349: }
1350:
1351: /* here y is not scalar */
1352: switch(ty)
1353: {
1354: case t_POL:
1355: v=varn(y);
1356: if (y==polun[v] || y==polx[v])
1357: err(talker,"trying to overwrite a universal polynomial");
1358: gaffect(x,(GEN)y[2]);
1359: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
1360: if (gcmp0(x)) y[1]=evallgef(2)+evalvarn(v);
1361: else y[1]=evalsigne(1)+evallgef(3)+evalvarn(v);
1362: break;
1363:
1364: case t_SER:
1365: v=varn(y); gaffect(x,(GEN)y[2]);
1366: if (gcmp0(x))
1367: y[1] = evalvalp(ly-2) | evalvarn(v);
1368: else
1369: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
1370: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
1371: break;
1372:
1373: case t_RFRAC: case t_RFRACN:
1374: gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;
1375:
1376: default: err(operf,"",tx,ty);
1377: }
1378: return;
1379: }
1380:
1381: if (is_const_t(ty)) {
1382: #define initial_value(ep) ((ep)+1)
1383: if (tx == t_POL) {
1384: entree *var = varentries[ordvar[varn(x)]];
1385: /* Is a bound expression, thus should not loop: */
1386: if (var && var->value != (void*)initial_value(var)) {
1387: gaffect(geval(x),y);
1388: return;
1389: }
1390: } else if (is_rfrac_t(tx)) {
1391: GEN num = (GEN)x[1], den = (GEN)x[2];
1392: entree *varnum = varentries[ordvar[varn(num)]];
1393: entree *varden = varentries[ordvar[varn(den)]];
1394:
1395: /* Are bound expressions, thus should not loop: */
1396: if (varnum && varnum->value != (void*)initial_value(varnum)
1397: && varden && varden->value != (void*)initial_value(varden)) {
1398: gaffect(geval(x),y);
1399: return;
1400: }
1401: }
1402: #undef initial_value
1403: err(operf,"",tx,ty);
1404: }
1405:
1406: switch(tx)
1407: {
1408: case t_POL:
1409: v=varn(x);
1410: switch(ty)
1411: {
1412: case t_POL:
1413: vy=varn(y); if (vy>v) err(operf,"",tx,ty);
1414: if (vy==v)
1415: {
1416: l=lgef(x); if (l>ly) err(operi,"",tx,ty);
1417: y[1]=x[1]; for (i=2; i<l; i++) gaffect((GEN)x[i],(GEN)y[i]);
1418: }
1419: else
1420: {
1421: gaffect(x,(GEN)y[2]);
1422: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
1423: if (signe(x)) y[1]=evalsigne(1)+evallgef(3)+evalvarn(vy);
1424: else y[1]=evallgef(2)+evalvarn(vy);
1425: }
1426: break;
1427:
1428: case t_SER:
1429: vy=varn(y); if (vy>v) err(operf,"",tx,ty);
1430: if (!signe(x)) { gaffsg(0,y); return; }
1431: if (vy==v)
1432: {
1433: i=gval(x,v); y[1]=evalvarn(v) | evalvalp(i) | evalsigne(1);
1434: k=lgef(x)-i; if (k>ly) k=ly;
1435: for (j=2; j<k; j++) gaffect((GEN)x[i+j],(GEN)y[j]);
1436: for ( ; j<ly; j++) gaffsg(0,(GEN)y[j]);
1437: }
1438: else
1439: {
1440: gaffect(x,(GEN)y[2]);
1441: if (!signe(x))
1442: y[1] = evalvalp(ly-2) | evalvarn(vy);
1443: else
1444: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
1445: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
1446: }
1447: break;
1448:
1449: case t_POLMOD:
1450: gmodz(x,(GEN)y[1],(GEN)y[2]); break;
1451:
1452: case t_RFRAC: case t_RFRACN:
1453: gaffect(x,(GEN)y[1]); gaffsg(1,(GEN)y[2]); break;
1454:
1455: default: err(operf,"",tx,ty);
1456: }
1457: break;
1458:
1459: case t_SER:
1460: if (ty!=t_SER) err(operf,"",tx,ty);
1461: v=varn(x); vy=varn(y); if (vy>v) err(operf,"",tx,ty);
1462: if (vy==v)
1463: {
1464: y[1]=x[1]; k=lx; if (k>ly) k=ly;
1465: for (i=2; i< k; i++) gaffect((GEN)x[i],(GEN)y[i]);
1466: for ( ; i<ly; i++) gaffsg(0,(GEN)y[i]);
1467: }
1468: else
1469: {
1470: gaffect(x,(GEN)y[2]);
1471: if (!signe(x))
1472: y[1] = evalvalp(ly-2) | evalvarn(vy);
1473: else
1474: y[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
1475: for (i=3; i<ly; i++) gaffsg(0,(GEN)y[i]);
1476: }
1477: break;
1478:
1479: case t_RFRAC: case t_RFRACN:
1480: switch(ty)
1481: {
1482: case t_POL: case t_VEC: case t_COL: case t_MAT:
1483: err(operf,"",tx,ty);
1484:
1485: case t_POLMOD:
1486: av=avma; p1=ginvmod((GEN)x[2],(GEN)y[1]);
1487: gmodz(gmul((GEN)x[1],p1),(GEN)y[1],(GEN)y[2]);
1488: avma=av; break;
1489:
1490: case t_SER:
1491: gdivz((GEN)x[1],(GEN)x[2],y); break;
1492:
1493: case t_RFRAC:
1494: if (tx==t_RFRACN) gredz(x,y);
1495: else
1496: {
1497: gaffect((GEN)x[1],(GEN)y[1]);
1498: gaffect((GEN)x[2],(GEN)y[2]);
1499: }
1500: break;
1501:
1502: case t_RFRACN:
1503: gaffect((GEN)x[1],(GEN)y[1]);
1504: gaffect((GEN)x[2],(GEN)y[2]); break;
1505:
1506: default: err(operf,"",tx,ty);
1507: }
1508: break;
1509:
1510: case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
1511: if (ty != tx || ly != lx) err(operi,"",tx,ty);
1512: for (i=1; i<lx; i++) gaffect((GEN)x[i],(GEN)y[i]);
1513: break;
1514:
1515: default: err(operf,"",tx,ty);
1516: }
1517: }
1518:
1519: /*******************************************************************/
1520: /* */
1521: /* CONVERSION QUAD --> REAL, COMPLEX OR P-ADIC */
1522: /* */
1523: /*******************************************************************/
1524:
1525: GEN
1526: co8(GEN x, long prec)
1527: {
1528: long av=avma,tetpil;
1529: GEN p1, p = (GEN) x[1];
1530:
1531: p1 = subii(sqri((GEN)p[3]), shifti((GEN)p[2],2));
1532: if (signe(p1) > 0)
1533: {
1534: p1 = subri(gsqrt(p1,prec), (GEN)p[3]);
1535: setexpo(p1, expo(p1)-1);
1536: }
1537: else
1538: {
1539: p1 = gsub(gsqrt(p1,prec), (GEN)p[3]);
1540: p1[1] = lmul2n((GEN)p1[1],-1);
1541: setexpo(p1[2], expo(p1[2])-1);
1542: }/* p1 = (-b + sqrt(D)) / 2 */
1543: p1 = gmul((GEN)x[3],p1); tetpil=avma;
1544: return gerepile(av,tetpil,gadd((GEN)x[2],p1));
1545: }
1546:
1547: GEN
1548: cvtop(GEN x, GEN p, long l)
1549: {
1550: GEN p1,p2,p3;
1551: long av,tetpil,n;
1552:
1553: if (typ(p)!=t_INT)
1554: err(talker,"not an integer modulus in cvtop or gcvtop");
1555: if (gcmp0(x)) return ggrandocp(p,l);
1556: switch(typ(x))
1557: {
1558: case t_INT:
1559: return gadd(x,ggrandocp(p,ggval(x,p)+l));
1560:
1561: case t_INTMOD:
1562: n=ggval((GEN)x[1],p); if (n>l) n=l;
1563: return gadd((GEN)x[2],ggrandocp(p,n));
1564:
1565: case t_FRAC: case t_FRACN:
1566: n = ggval((GEN)x[1],p) - ggval((GEN)x[2],p);
1567: return gadd(x,ggrandocp(p,n+l));
1568:
1569: case t_COMPLEX:
1570: av=avma; p1=gsqrt(gaddgs(ggrandocp(p,l),-1),0);
1571: p1=gmul(p1,(GEN)x[2]); tetpil=avma;
1572: return gerepile(av,tetpil,gadd(p1,(GEN)x[1]));
1573:
1574: case t_PADIC:
1575: return gprec(x,l);
1576:
1577: case t_QUAD:
1578: av=avma; p1=(GEN)x[1]; p3=gmul2n((GEN)p1[3],-1);
1579: p2=gsub(gsqr(p3),(GEN)p1[2]);
1580: switch(typ(p2))
1581: {
1582: case t_INT:
1583: n=ggval(p2,p); p2=gadd(p2,ggrandocp(p,n+l)); break;
1584:
1585: case t_INTMOD: case t_PADIC:
1586: break;
1587:
1588: case t_FRAC: case t_FRACN:
1589: n = ggval((GEN)p2[1],p) - ggval((GEN)p2[2],p);
1590: p2=gadd(p2,ggrandocp(p,n+l)); break;
1591:
1592: default: err(operi,"",t_QUAD,t_QUAD);
1593: }
1594: p2=gsqrt(p2,0); p1=gmul((GEN)x[3],gsub(p2,p3)); tetpil=avma;
1595: return gerepile(av,tetpil,gadd((GEN)x[2],p1));
1596:
1597: }
1598: err(typeer,"cvtop");
1599: return NULL; /* not reached */
1600: }
1601:
1602: GEN
1603: gcvtop(GEN x, GEN p, long r)
1604: {
1605: long i,lx, tx=typ(x);
1606: GEN y;
1607:
1608: if (is_const_t(tx)) return cvtop(x,p,r);
1609: switch(tx)
1610: {
1611: case t_POL:
1612: lx=lgef(x); y=cgetg(lx,t_POL); y[1]=x[1];
1613: for (i=2; i<lx; i++)
1614: y[i]=(long)gcvtop((GEN)x[i],p,r);
1615: break;
1616:
1617: case t_SER:
1618: lx=lg(x); y=cgetg(lx,t_SER); y[1]=x[1];
1619: for (i=2; i<lx; i++)
1620: y[i]=(long)gcvtop((GEN)x[i],p,r);
1621: break;
1622:
1623: case t_POLMOD: case t_RFRAC: case t_RFRACN:
1624: case t_VEC: case t_COL: case t_MAT:
1625: lx=lg(x); y=cgetg(lx,tx);
1626: for (i=1; i<lx; i++)
1627: y[i]=(long)gcvtop((GEN)x[i],p,r);
1628: break;
1629:
1630: default: err(typeer,"gcvtop");
1631: return NULL; /* not reached */
1632: }
1633: return y;
1634: }
1635:
1636: long
1637: gexpo(GEN x)
1638: {
1639: long tx=typ(x),lx,e,i,y,av;
1640:
1641: switch(tx)
1642: {
1643: case t_INT:
1644: return expi(x);
1645:
1646: case t_FRAC: case t_FRACN:
1647: if (!signe(x[1])) return -HIGHEXPOBIT;
1648: return expi((GEN)x[1]) - expi((GEN)x[2]);
1649:
1650: case t_REAL:
1651: return expo(x);
1652:
1653: case t_COMPLEX:
1654: e = gexpo((GEN)x[1]);
1655: y = gexpo((GEN)x[2]); return max(e,y);
1656:
1657: case t_QUAD:
1658: av=avma; e = gexpo(co8(x,3)); avma=av; return e;
1659:
1660: case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
1661: lx=(tx==t_POL)? lgef(x): lg(x);
1662: y = -HIGHEXPOBIT;
1663: for (i=lontyp[tx]; i<lx; i++) { e=gexpo((GEN)x[i]); if (e>y) y=e; }
1664: return y;
1665: }
1666: err(typeer,"gexpo");
1667: return 0; /* not reached */
1668: }
1669:
1670: long
1671: sizedigit(GEN x)
1672: {
1673: return gcmp0(x)? 0: (long) ((gexpo(x)+1) * L2SL10) + 1;
1674: }
1675:
1676: /* Normalize series x in place.
1677: * Assumption: x,x[2],...,x[lg(x)-1] have been created in that order.
1678: * All intermediate objects will be destroyed.
1679: */
1680: GEN
1681: normalize(GEN x)
1682: {
1683: long i,j, lx = lg(x);
1684:
1685: if (typ(x)!=t_SER) err(typeer,"normalize");
1686: if (lx==2) { setsigne(x,0); avma = (long) x; return x; }
1687: if (! isexactzero((GEN)x[2])) { setsigne(x,1); return x; }
1688:
1689: for (i=3; i<lx; i++)
1690: if (! isexactzero((GEN)x[i]))
1691: {
1692: long tetpil = avma;
1693: GEN p1 = cgetg(lx-i+2,t_SER);
1694: p1[1] = evalsigne(1) | evalvalp(valp(x)+i-2) | evalvarn(varn(x));
1695: j=i; i=2; while (j<lx) p1[i++] = lcopy((GEN)x[j++]);
1696: return gerepile((long) (x+lx),tetpil,p1);
1697: }
1698: avma = (long) (x+lx); return zeroser(varn(x),lx-2);
1699: }
1700:
1701: GEN
1702: normalizepol_i(GEN x, long lx)
1703: {
1704: long i;
1705: for (i=lx-1; i>1; i--)
1706: if (! isexactzero((GEN)x[i])) break;
1707: setlgef(x,i+1);
1708: for (; i>1; i--)
1709: if (! gcmp0((GEN)x[i]) ) { setsigne(x,1); return x; }
1710: setsigne(x,0); return x;
1711: }
1712:
1713: /* Normalize polynomial x in place. See preceding comment */
1714: GEN
1715: normalizepol(GEN x)
1716: {
1717: if (typ(x)!=t_POL) err(typeer,"normalizepol");
1718: return normalizepol_i(x, lgef(x));
1719: }
1720:
1721: int
1722: gsigne(GEN x)
1723: {
1724: switch(typ(x))
1725: {
1726: case t_INT: case t_REAL:
1727: return signe(x);
1728:
1729: case t_FRAC: case t_FRACN:
1730: return (signe(x[2])>0) ? signe(x[1]) : -signe(x[1]);
1731: }
1732: err(typeer,"gsigne");
1733: return 0; /* not reached */
1734: }
1735:
1736: /*******************************************************************/
1737: /* */
1738: /* LISTS */
1739: /* */
1740: /*******************************************************************/
1741:
1742: GEN
1743: listcreate(long n)
1744: {
1745: GEN list;
1746:
1747: if (n<0) err(talker,"negative length in listcreate");
1748: n += 2;
1749: if (n & ~LGEFBITS) err(talker,"list too long (max = %ld)",LGEFBITS-2);
1750: list=cgetg(n,t_LIST); list[1]=evallgef(2);
1751: return list;
1752: }
1753:
1754: static void
1755: listaffect(GEN list, long index, GEN object)
1756: {
1757: if (index<lgef(list) && isclone(list[index]))
1758: gunclone((GEN)list[index]);
1759: list[index]=lclone(object);
1760: }
1761:
1762: void
1763: listkill(GEN list)
1764: {
1765: long lx,i;
1766:
1767: if (typ(list)!=t_LIST) err(typeer,"listkill");
1768: lx=lgef(list);
1769: for (i=2; i<lx; i++)
1770: if (isclone(list[i])) gunclone((GEN)list[i]);
1771: list[1]=evallgef(2); return;
1772: }
1773:
1774: GEN
1775: listput(GEN list, GEN object, long index)
1776: {
1777: long lx=lgef(list);
1778:
1779: if (typ(list)!=t_LIST) err(typeer,"listput");
1780: if (index < 0) err(talker,"negative index (%ld) in listput", index);
1781: if (!index || index >= lx-1)
1782: {
1783: index = lx-1; lx++;
1784: if (lx > lg(list))
1785: err(talker,"no more room in this list (size %ld)", lg(list)-2);
1786: }
1787: listaffect(list,index+1,object);
1788: list[1]=evallgef(lx);
1789: return (GEN)list[index+1];
1790: }
1791:
1792: GEN
1793: listinsert(GEN list, GEN object, long index)
1794: {
1795: long lx = lgef(list), i;
1796:
1797: if (typ(list)!=t_LIST) err(typeer,"listinsert");
1798: if (index <= 0 || index > lx-1) err(talker,"bad index in listinsert");
1799: lx++; if (lx > lg(list)) err(talker,"no more room in this list");
1800:
1801: for (i=lx-2; i > index; i--) list[i+1]=list[i];
1802: list[index+1] = lclone(object);
1803: list[1] = evallgef(lx); return (GEN)list[index+1];
1804: }
1805:
1806: GEN
1807: gtolist(GEN x)
1808: {
1809: long tx,lx,i;
1810: GEN list;
1811:
1812: if (!x) { list = cgetg(2, t_LIST); list[1] = evallgef(2); return list; }
1813:
1814: tx = typ(x);
1815: lx = (tx==t_LIST)? lgef(x): lg(x);
1816: switch(tx)
1817: {
1818: case t_VEC: case t_COL:
1819: lx++; x--; /* fall through */
1820: case t_LIST:
1821: list = cgetg(lx,t_LIST);
1822: for (i=2; i<lx; i++) list[i] = lclone((GEN)x[i]);
1823: break;
1824: default: err(typeer,"gtolist");
1825: return NULL; /* not reached */
1826: }
1827: list[1] = evallgef(lx); return list;
1828: }
1829:
1830: GEN
1831: listconcat(GEN list1, GEN list2)
1832: {
1833: long i,l1,lx;
1834: GEN list;
1835:
1836: if (typ(list1)!=t_LIST || typ(list2)!=t_LIST)
1837: err(typeer,"listconcat");
1838: l1=lgef(list1)-2; lx=l1+lgef(list2);
1839: list = listcreate(lx-2);
1840: for (i=2; i<=l1+1; i++) listaffect(list,i,(GEN)list1[i]);
1841: for ( ; i<lx; i++) listaffect(list,i,(GEN)list2[i-l1]);
1842: list[1]=evallgef(lx); return list;
1843: }
1844:
1845: GEN
1846: listsort(GEN list, long flag)
1847: {
1848: long i,av=avma, c=list[1], lx = lgef(list)-1;
1849: GEN perm,vec,l;
1850:
1851: if (typ(list) != t_LIST) err(typeer,"listsort");
1852: vec = list+1;
1853: vec[0] = evaltyp(t_VEC) | evallg(lx);
1854: perm = sindexsort(vec);
1855: list[1] = c; l = cgetg(lx,t_VEC);
1856: for (i=1; i<lx; i++) l[i] = vec[perm[i]];
1857: if (flag)
1858: {
1859: c=1; vec[1] = l[1];
1860: for (i=2; i<lx; i++)
1861: if (!gegal((GEN)l[i], (GEN)vec[c]))
1862: vec[++c] = l[i];
1863: else
1864: if (isclone(l[i])) gunclone((GEN)l[i]);
1865: setlgef(list,c+2);
1866: }
1867: else
1868: for (i=1; i<lx; i++) vec[i] = l[i];
1869:
1870: avma=av; return list;
1871: }