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