1 /* Copyright (C) 2000 The PARI group.
2
3 This file is part of the PARI/GP package.
4
5 PARI/GP is free software; you can redistribute it and/or modify it under the
6 terms of the GNU General Public License as published by the Free Software
7 Foundation; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 /********************************************************************/
16 /** **/
17 /** GENERIC OPERATIONS **/
18 /** (second part) **/
19 /** **/
20 /********************************************************************/
21 #include "pari.h"
22 #include "paripriv.h"
23
24 /*********************************************************************/
25 /** **/
26 /** MAP FUNCTIONS WITH GIVEN PROTOTYPES **/
27 /** **/
28 /*********************************************************************/
29 GEN
map_proto_G(GEN (* f)(GEN),GEN x)30 map_proto_G(GEN (*f)(GEN), GEN x)
31 {
32 if (is_matvec_t(typ(x)))
33 {
34 long lx, i;
35 GEN y = cgetg_copy(x, &lx);
36 for (i=1; i<lx; i++) gel(y,i) = map_proto_G(f, gel(x,i));
37 return y;
38 }
39 return f(x);
40 }
41
42 GEN
map_proto_lG(long (* f)(GEN),GEN x)43 map_proto_lG(long (*f)(GEN), GEN x)
44 {
45 if (is_matvec_t(typ(x)))
46 {
47 long lx, i;
48 GEN y = cgetg_copy(x, &lx);
49 for (i=1; i<lx; i++) gel(y,i) = map_proto_lG(f, gel(x,i));
50 return y;
51 }
52 return stoi(f(x));
53 }
54
55 GEN
map_proto_lGL(long (* f)(GEN,long),GEN x,long y)56 map_proto_lGL(long (*f)(GEN,long), GEN x, long y)
57 {
58 if (is_matvec_t(typ(x)))
59 {
60 long l, i;
61 GEN t = cgetg_copy(x, &l);
62 for (i=1; i<l; i++) gel(t,i) = map_proto_lGL(f,gel(x,i),y);
63 return t;
64 }
65 return stoi(f(x,y));
66 }
67
68 static GEN
_domul(void * data,GEN x,GEN y)69 _domul(void *data, GEN x, GEN y)
70 {
71 GEN (*mul)(GEN,GEN)=(GEN (*)(GEN,GEN)) data;
72 return mul(x,y);
73 }
74
75 GEN
gassoc_proto(GEN (* f)(GEN,GEN),GEN x,GEN y)76 gassoc_proto(GEN (*f)(GEN,GEN), GEN x, GEN y)
77 {
78 if (!y)
79 {
80 pari_sp av = avma;
81 switch(typ(x))
82 {
83 case t_LIST:
84 x = list_data(x); if (!x) return gen_1;
85 case t_VEC:
86 case t_COL: break;
87 default: pari_err_TYPE("association",x);
88 }
89 return gerepileupto(av, gen_product(x, (void *)f, _domul));
90
91 }
92 return f(x,y);
93 }
94 /*******************************************************************/
95 /* */
96 /* CREATION OF A P-ADIC GEN */
97 /* */
98 /*******************************************************************/
99 GEN
cgetp(GEN x)100 cgetp(GEN x)
101 {
102 GEN y = cgetg(5,t_PADIC);
103 y[1] = (x[1]&PRECPBITS) | _evalvalp(0);
104 gel(y,2) = icopy(gel(x,2));
105 gel(y,3) = icopy(gel(x,3));
106 gel(y,4) = cgeti(lgefint(gel(x,3))); return y;
107 }
108
109 /*******************************************************************/
110 /* */
111 /* SIZES */
112 /* */
113 /*******************************************************************/
114
115 long
glength(GEN x)116 glength(GEN x)
117 {
118 long tx = typ(x);
119 switch(tx)
120 {
121 case t_INT: return lgefint(x)-2;
122 case t_LIST: {
123 GEN L = list_data(x);
124 return L? lg(L)-1: 0;
125 }
126 case t_REAL: return signe(x)? lg(x)-2: 0;
127 case t_STR: return strlen( GSTR(x) );
128 case t_VECSMALL: return lg(x)-1;
129 }
130 return lg(x) - lontyp[tx];
131 }
132
133 GEN
matsize(GEN x)134 matsize(GEN x)
135 {
136 long L = lg(x) - 1;
137 switch(typ(x))
138 {
139 case t_VEC: return mkvec2s(1, L);
140 case t_COL: return mkvec2s(L, 1);
141 case t_MAT: return mkvec2s(L? nbrows(x): 0, L);
142 }
143 pari_err_TYPE("matsize",x);
144 return NULL; /* LCOV_EXCL_LINE */
145 }
146
147 /*******************************************************************/
148 /* */
149 /* CONVERSION GEN --> long */
150 /* */
151 /*******************************************************************/
152
153 long
gtolong(GEN x)154 gtolong(GEN x)
155 {
156 switch(typ(x))
157 {
158 case t_INT:
159 return itos(x);
160 case t_REAL:
161 return (long)(rtodbl(x) + 0.5);
162 case t_FRAC:
163 { pari_sp av = avma; return gc_long(av, itos(ground(x))); }
164 case t_COMPLEX:
165 if (gequal0(gel(x,2))) return gtolong(gel(x,1)); break;
166 case t_QUAD:
167 if (gequal0(gel(x,3))) return gtolong(gel(x,2)); break;
168 }
169 pari_err_TYPE("gtolong",x);
170 return 0; /* LCOV_EXCL_LINE */
171 }
172
173 /*******************************************************************/
174 /* */
175 /* COMPARISONS */
176 /* */
177 /*******************************************************************/
178 int
isexactzero(GEN g)179 isexactzero(GEN g)
180 {
181 long i, lx;
182 switch (typ(g))
183 {
184 case t_INT:
185 return !signe(g);
186 case t_INTMOD:
187 return !signe(gel(g,2));
188 case t_COMPLEX:
189 return isexactzero(gel(g,1)) && isexactzero(gel(g,2));
190 case t_FFELT:
191 return FF_equal0(g);
192 case t_QUAD:
193 return isexactzero(gel(g,2)) && isexactzero(gel(g,3));
194 case t_POLMOD:
195 return isexactzero(gel(g,2));
196 case t_POL:
197 lx = lg(g); /* cater for Mod(0,2)*x^0 */
198 return lx == 2 || (lx == 3 && isexactzero(gel(g,2)));
199 case t_RFRAC:
200 return isexactzero(gel(g,1)); /* may occur: Mod(0,2)/x */
201 case t_VEC: case t_COL: case t_MAT:
202 for (i=lg(g)-1; i; i--)
203 if (!isexactzero(gel(g,i))) return 0;
204 return 1;
205 }
206 return 0;
207 }
208 GEN
gisexactzero(GEN g)209 gisexactzero(GEN g)
210 {
211 long i, lx;
212 GEN a, b;
213 switch (typ(g))
214 {
215 case t_INT:
216 return !signe(g)? g: NULL;
217 case t_INTMOD:
218 return !signe(gel(g,2))? g: NULL;
219 case t_COMPLEX:
220 a = gisexactzero(gel(g,1)); if (!a) return NULL;
221 b = gisexactzero(gel(g,2)); if (!b) return NULL;
222 return ggcd(a,b);
223 case t_FFELT:
224 return FF_equal0(g)? g: NULL;
225 case t_QUAD:
226 a = gisexactzero(gel(g,2)); if (!a) return NULL;
227 b = gisexactzero(gel(g,3)); if (!b) return NULL;
228 return ggcd(a,b);
229 case t_POLMOD:
230 return gisexactzero(gel(g,2));
231 case t_POL:
232 lx = lg(g); /* cater for Mod(0,2)*x^0 */
233 if (lx == 2) return gen_0;
234 if (lx == 3) return gisexactzero(gel(g,2));
235 return NULL;
236 case t_RFRAC:
237 return gisexactzero(gel(g,1)); /* may occur: Mod(0,2)/x */
238 case t_VEC: case t_COL: case t_MAT:
239 a = gen_0;
240 for (i=lg(g)-1; i; i--)
241 {
242 b = gisexactzero(gel(g,i));
243 if (!b) return NULL;
244 a = ggcd(a, b);
245 }
246 return a;
247 }
248 return NULL;
249 }
250
251 int
isrationalzero(GEN g)252 isrationalzero(GEN g)
253 {
254 long i;
255 switch (typ(g))
256 {
257 case t_INT:
258 return !signe(g);
259 case t_COMPLEX:
260 return isintzero(gel(g,1)) && isintzero(gel(g,2));
261 case t_QUAD:
262 return isintzero(gel(g,2)) && isintzero(gel(g,3));
263 case t_POLMOD:
264 return isrationalzero(gel(g,2));
265 case t_POL: return lg(g) == 2;
266 case t_VEC: case t_COL: case t_MAT:
267 for (i=lg(g)-1; i; i--)
268 if (!isrationalzero(gel(g,i))) return 0;
269 return 1;
270 }
271 return 0;
272 }
273
274 int
gequal0(GEN x)275 gequal0(GEN x)
276 {
277 switch(typ(x))
278 {
279 case t_INT: case t_REAL: case t_POL: case t_SER:
280 return !signe(x);
281
282 case t_INTMOD:
283 return !signe(gel(x,2));
284
285 case t_FFELT:
286 return FF_equal0(x);
287
288 case t_COMPLEX:
289 /* is 0 iff norm(x) would be 0 (can happen with Re(x) and Im(x) != 0
290 * only if Re(x) and Im(x) are of type t_REAL). See mp.c:addrr().
291 */
292 if (gequal0(gel(x,1)))
293 {
294 if (gequal0(gel(x,2))) return 1;
295 if (typ(gel(x,1))!=t_REAL || typ(gel(x,2))!=t_REAL) return 0;
296 return (expo(gel(x,1))>=expo(gel(x,2)));
297 }
298 if (gequal0(gel(x,2)))
299 {
300 if (typ(gel(x,1))!=t_REAL || typ(gel(x,2))!=t_REAL) return 0;
301 return (expo(gel(x,2))>=expo(gel(x,1)));
302 }
303 return 0;
304
305 case t_PADIC:
306 return !signe(gel(x,4));
307
308 case t_QUAD:
309 return gequal0(gel(x,2)) && gequal0(gel(x,3));
310
311 case t_POLMOD:
312 return gequal0(gel(x,2));
313
314 case t_RFRAC:
315 return gequal0(gel(x,1));
316
317 case t_VEC: case t_COL: case t_MAT:
318 {
319 long i;
320 for (i=lg(x)-1; i; i--)
321 if (!gequal0(gel(x,i))) return 0;
322 return 1;
323 }
324 }
325 return 0;
326 }
327
328 /* x a t_POL or t_SER, considered as having valuation v; let X(t) = t^(-v) x(t)
329 * return 1 (true) if coeff(X,i) = 0 for all i != 0 and test(coeff(X, 0))
330 * is true. Return 0 (false) otherwise, or if x == 0 */
331 static int
is_monomial_test(GEN x,long v,int (* test)(GEN))332 is_monomial_test(GEN x, long v, int(*test)(GEN))
333 {
334 long d, i, l;
335 if (!signe(x)) return (typ(x) == t_SER && v <= 0);
336 if (v > 0) return 0;
337 l = lg(x); d = 2-v;
338 if (l <= d) return 0;
339 /* 2 <= d < l */
340 if (!test(gel(x,d))) return 0;
341 for (i = 2; i < d; i++)
342 if (!gequal0(gel(x,i))) return 0;
343 for (i = d+1; i < l; i++)
344 if (!gequal0(gel(x,i))) return 0;
345 return 1;
346 }
347 static int
col_test(GEN x,int (* test)(GEN))348 col_test(GEN x, int(*test)(GEN))
349 {
350 long i, l = lg(x);
351 if (l == 1 || !test(gel(x,1))) return 0;
352 for (i = 2; i < l; i++)
353 if (!gequal0(gel(x,i))) return 0;
354 return 1;
355 }
356 static int
mat_test(GEN x,int (* test)(GEN))357 mat_test(GEN x, int(*test)(GEN))
358 {
359 long i, j, l = lg(x);
360 if (l == 1) return 1;
361 if (l != lgcols(x)) return 0;
362 for (i = 1; i < l; i++)
363 for (j = 1; j < l; j++)
364 if (i == j) {
365 if (!test(gcoeff(x,i,i))) return 0;
366 } else {
367 if (!gequal0(gcoeff(x,i,j))) return 0;
368 }
369 return 1;
370 }
371
372 /* returns 1 whenever x = 1, and 0 otherwise */
373 int
gequal1(GEN x)374 gequal1(GEN x)
375 {
376 switch(typ(x))
377 {
378 case t_INT:
379 return equali1(x);
380
381 case t_REAL:
382 {
383 long s = signe(x);
384 if (!s) return expo(x) >= 0;
385 return s > 0 ? absrnz_equal1(x): 0;
386 }
387 case t_INTMOD:
388 return is_pm1(gel(x,2)) || is_pm1(gel(x,1));
389 case t_POLMOD:
390 return gequal1(gel(x,2)) || gequal1(gel(x,1));
391
392 case t_FFELT:
393 return FF_equal1(x);
394
395 case t_FRAC:
396 return 0;
397
398 case t_COMPLEX:
399 return gequal1(gel(x,1)) && gequal0(gel(x,2));
400
401 case t_PADIC:
402 return !valp(x) && gequal1(gel(x,4));
403
404 case t_QUAD:
405 return gequal1(gel(x,2)) && gequal0(gel(x,3));
406
407 case t_POL: return is_monomial_test(x, 0, &gequal1);
408 case t_SER: return is_monomial_test(x, valp(x), &gequal1);
409
410 case t_RFRAC: return gequal(gel(x,1), gel(x,2));
411 case t_COL: return col_test(x, &gequal1);
412 case t_MAT: return mat_test(x, &gequal1);
413 }
414 return 0;
415 }
416
417 /* returns 1 whenever the x = -1, 0 otherwise */
418 int
gequalm1(GEN x)419 gequalm1(GEN x)
420 {
421 pari_sp av;
422 GEN p1;
423
424 switch(typ(x))
425 {
426 case t_INT:
427 return equalim1(x);
428
429 case t_REAL:
430 {
431 long s = signe(x);
432 if (!s) return expo(x) >= 0;
433 return s < 0 ? absrnz_equal1(x): 0;
434 }
435 case t_INTMOD:
436 av = avma; return gc_bool(av, equalii(addui(1,gel(x,2)), gel(x,1)));
437
438 case t_FRAC:
439 return 0;
440
441 case t_FFELT:
442 return FF_equalm1(x);
443
444 case t_COMPLEX:
445 return gequalm1(gel(x,1)) && gequal0(gel(x,2));
446
447 case t_QUAD:
448 return gequalm1(gel(x,2)) && gequal0(gel(x,3));
449
450 case t_PADIC:
451 av = avma; return gc_bool(av, equalii(addui(1,gel(x,4)), gel(x,3)));
452
453 case t_POLMOD:
454 av = avma; p1 = gaddgs(gel(x,2), 1);
455 return gc_bool(av, gequal0(p1) || gequal(p1,gel(x,1)));
456
457 case t_POL: return is_monomial_test(x, 0, &gequalm1);
458 case t_SER: return is_monomial_test(x, valp(x), &gequalm1);
459
460 case t_RFRAC:
461 av = avma; return gc_bool(av, gequal(gel(x,1), gneg_i(gel(x,2))));
462 case t_COL: return col_test(x, &gequalm1);
463 case t_MAT: return mat_test(x, &gequalm1);
464 }
465 return 0;
466 }
467
468 int
gequalX(GEN x)469 gequalX(GEN x) { return typ(x) == t_POL && lg(x) == 4
470 && isintzero(gel(x,2)) && isint1(gel(x,3)); }
471
472 static int
cmp_str(const char * x,const char * y)473 cmp_str(const char *x, const char *y)
474 {
475 int f = strcmp(x, y);
476 return f > 0? 1
477 : f? -1: 0;
478 }
479
480 static int
cmp_universal_rec(GEN x,GEN y,long i0)481 cmp_universal_rec(GEN x, GEN y, long i0)
482 {
483 long i, lx = lg(x), ly = lg(y);
484 if (lx < ly) return -1;
485 if (lx > ly) return 1;
486 for (i = i0; i < lx; i++)
487 {
488 int f = cmp_universal(gel(x,i), gel(y,i));
489 if (f) return f;
490 }
491 return 0;
492 }
493 /* Universal "meaningless" comparison function. Transitive, returns 0 iff
494 * gidentical(x,y) */
495 int
cmp_universal(GEN x,GEN y)496 cmp_universal(GEN x, GEN y)
497 {
498 long lx, ly, i, tx = typ(x), ty = typ(y);
499
500 if (tx < ty) return -1;
501 if (ty < tx) return 1;
502 switch(tx)
503 {
504 case t_INT: return cmpii(x,y);
505 case t_STR: return cmp_str(GSTR(x),GSTR(y));
506 case t_REAL:
507 case t_VECSMALL:
508 lx = lg(x);
509 ly = lg(y);
510 if (lx < ly) return -1;
511 if (lx > ly) return 1;
512 for (i = 1; i < lx; i++)
513 {
514 if (x[i] < y[i]) return -1;
515 if (x[i] > y[i]) return 1;
516 }
517 return 0;
518
519 case t_POL:
520 case t_SER:
521 case t_FFELT:
522 case t_CLOSURE:
523 if (x[1] < y[1]) return -1;
524 if (x[1] > y[1]) return 1;
525 return cmp_universal_rec(x, y, 2);
526
527 case t_LIST:
528 {
529 long tx = list_typ(x), ty = list_typ(y);
530 GEN vx, vy;
531 pari_sp av;
532 if (tx < ty) return -1;
533 if (tx > ty) return 1;
534 vx = list_data(x);
535 vy = list_data(y);
536 if (!vx) return vy? -1: 0;
537 if (!vy) return 1;
538 av = avma;
539 if (tx == t_LIST_MAP)
540 {
541 vx = maptomat_shallow(x);
542 vy = maptomat_shallow(y);
543 }
544 return gc_int(av, cmp_universal_rec(vx, vy, 1));
545 }
546 default:
547 return cmp_universal_rec(x, y, lontyp[tx]);
548 }
549 }
550
551 static int
cmpfrac(GEN x,GEN y)552 cmpfrac(GEN x, GEN y)
553 {
554 pari_sp av = avma;
555 GEN a = gel(x,1), b = gel(x,2);
556 GEN c = gel(y,1), d = gel(y,2);
557 return gc_bool(av, cmpii(mulii(a, d), mulii(b, c)));
558 }
559 static int
cmpifrac(GEN a,GEN y)560 cmpifrac(GEN a, GEN y)
561 {
562 pari_sp av = avma;
563 GEN c = gel(y,1), d = gel(y,2);
564 return gc_int(av, cmpii(mulii(a, d), c));
565 }
566 static int
cmprfrac(GEN a,GEN y)567 cmprfrac(GEN a, GEN y)
568 {
569 pari_sp av = avma;
570 GEN c = gel(y,1), d = gel(y,2);
571 return gc_int(av, cmpri(mulri(a, d), c));
572 }
573 static int
cmpgen(GEN x,GEN y)574 cmpgen(GEN x, GEN y)
575 {
576 pari_sp av = avma;
577 return gc_int(av, gsigne(gsub(x,y)));
578 }
579
580 /* returns the sign of x - y when it makes sense. 0 otherwise */
581 int
gcmp(GEN x,GEN y)582 gcmp(GEN x, GEN y)
583 {
584 long tx = typ(x), ty = typ(y);
585
586 if (tx == ty) /* generic case */
587 switch(tx)
588 {
589 case t_INT: return cmpii(x, y);
590 case t_REAL: return cmprr(x, y);
591 case t_FRAC: return cmpfrac(x, y);
592 case t_QUAD: return cmpgen(x, y);
593 case t_STR: return cmp_str(GSTR(x), GSTR(y));
594 case t_INFINITY:
595 {
596 long sx = inf_get_sign(x), sy = inf_get_sign(y);
597 if (sx < sy) return -1;
598 if (sx > sy) return 1;
599 return 0;
600 }
601 }
602 if (ty == t_INFINITY) return -inf_get_sign(y);
603 switch(tx)
604 {
605 case t_INT:
606 switch(ty)
607 {
608 case t_REAL: return cmpir(x, y);
609 case t_FRAC: return cmpifrac(x, y);
610 case t_QUAD: return cmpgen(x, y);
611 }
612 break;
613 case t_REAL:
614 switch(ty)
615 {
616 case t_INT: return cmpri(x, y);
617 case t_FRAC: return cmprfrac(x, y);
618 case t_QUAD: return cmpgen(x, y);
619 }
620 break;
621 case t_FRAC:
622 switch(ty)
623 {
624 case t_INT: return -cmpifrac(y, x);
625 case t_REAL: return -cmprfrac(y, x);
626 case t_QUAD: return cmpgen(x, y);
627 }
628 break;
629 case t_QUAD:
630 return cmpgen(x, y);
631 case t_INFINITY: return inf_get_sign(x);
632 }
633 pari_err_TYPE2("comparison",x,y);
634 return 0;/*LCOV_EXCL_LINE*/
635 }
636
637 int
gcmpsg(long s,GEN y)638 gcmpsg(long s, GEN y)
639 {
640 switch(typ(y))
641 {
642 case t_INT: return cmpsi(s,y);
643 case t_REAL: return cmpsr(s,y);
644 case t_FRAC: {
645 pari_sp av = avma;
646 return gc_int(av, cmpii(mulsi(s,gel(y,2)), gel(y,1)));
647 }
648 case t_QUAD: {
649 pari_sp av = avma;
650 return gc_int(av, gsigne(gsubsg(s, y)));
651 }
652 case t_INFINITY: return -inf_get_sign(y);
653 }
654 pari_err_TYPE2("comparison",stoi(s),y);
655 return 0; /* LCOV_EXCL_LINE */
656 }
657
658 static long
roughtype(GEN x)659 roughtype(GEN x)
660 {
661 switch(typ(x))
662 {
663 case t_MAT: return t_MAT;
664 case t_VEC: case t_COL: return t_VEC;
665 case t_VECSMALL: return t_VECSMALL;
666 default: return t_INT;
667 }
668 }
669
670 static int lexcmpsg(long x, GEN y);
lexcmpgs(GEN x,long y)671 static int lexcmpgs(GEN x, long y) { return -lexcmpsg(y,x); }
672 /* lexcmp(stoi(x),y), y t_VEC/t_COL/t_MAT */
673 static int
lexcmp_s_matvec(long x,GEN y)674 lexcmp_s_matvec(long x, GEN y)
675 {
676 int fl;
677 if (lg(y)==1) return 1;
678 fl = lexcmpsg(x,gel(y,1));
679 if (fl) return fl;
680 return -1;
681 }
682 /* x a scalar, y a t_VEC/t_COL/t_MAT */
683 static int
lexcmp_scal_matvec(GEN x,GEN y)684 lexcmp_scal_matvec(GEN x, GEN y)
685 {
686 int fl;
687 if (lg(y)==1) return 1;
688 fl = lexcmp(x,gel(y,1));
689 if (fl) return fl;
690 return -1;
691 }
692 /* x a scalar, y a t_VECSMALL */
693 static int
lexcmp_scal_vecsmall(GEN x,GEN y)694 lexcmp_scal_vecsmall(GEN x, GEN y)
695 {
696 int fl;
697 if (lg(y)==1) return 1;
698 fl = lexcmpgs(x, y[1]);
699 if (fl) return fl;
700 return -1;
701 }
702
703 /* tx = ty = t_MAT, or x and y are both vect_t */
704 static int
lexcmp_similar(GEN x,GEN y)705 lexcmp_similar(GEN x, GEN y)
706 {
707 long i, lx = lg(x), ly = lg(y), l = minss(lx,ly);
708 for (i=1; i<l; i++)
709 {
710 int fl = lexcmp(gel(x,i),gel(y,i));
711 if (fl) return fl;
712 }
713 if (lx == ly) return 0;
714 return (lx < ly)? -1 : 1;
715 }
716 /* x a t_VECSMALL, y a t_VEC/t_COL ~ lexcmp_similar */
717 static int
lexcmp_vecsmall_vec(GEN x,GEN y)718 lexcmp_vecsmall_vec(GEN x, GEN y)
719 {
720 long i, lx = lg(x), ly = lg(y), l = minss(lx,ly);
721 for (i=1; i<l; i++)
722 {
723 int fl = lexcmpsg(x[i], gel(y,i));
724 if (fl) return fl;
725 }
726 if (lx == ly) return 0;
727 return (lx < ly)? -1 : 1;
728 }
729
730 /* x t_VEC/t_COL, y t_MAT */
731 static int
lexcmp_vec_mat(GEN x,GEN y)732 lexcmp_vec_mat(GEN x, GEN y)
733 {
734 int fl;
735 if (lg(x)==1) return -1;
736 if (lg(y)==1) return 1;
737 fl = lexcmp_similar(x,gel(y,1));
738 if (fl) return fl;
739 return -1;
740 }
741 /* x t_VECSMALl, y t_MAT ~ lexcmp_vec_mat */
742 static int
lexcmp_vecsmall_mat(GEN x,GEN y)743 lexcmp_vecsmall_mat(GEN x, GEN y)
744 {
745 int fl;
746 if (lg(x)==1) return -1;
747 if (lg(y)==1) return 1;
748 fl = lexcmp_vecsmall_vec(x, gel(y,1));
749 if (fl) return fl;
750 return -1;
751 }
752
753 /* x a t_VECSMALL, not y */
754 static int
lexcmp_vecsmall_other(GEN x,GEN y,long ty)755 lexcmp_vecsmall_other(GEN x, GEN y, long ty)
756 {
757 switch(ty)
758 {
759 case t_MAT: return lexcmp_vecsmall_mat(x, y);
760 case t_VEC: return lexcmp_vecsmall_vec(x, y);
761 default: return -lexcmp_scal_vecsmall(y, x); /*y scalar*/
762 }
763 }
764
765 /* lexcmp(stoi(s), y) */
766 static int
lexcmpsg(long x,GEN y)767 lexcmpsg(long x, GEN y)
768 {
769 switch(roughtype(y))
770 {
771 case t_MAT:
772 case t_VEC:
773 return lexcmp_s_matvec(x,y);
774 case t_VECSMALL: /* ~ lexcmp_scal_matvec */
775 if (lg(y)==1) return 1;
776 return (x > y[1])? 1: -1;
777 default: return gcmpsg(x,y);
778 }
779 }
780
781 /* as gcmp for vector/matrices, using lexicographic ordering on components */
782 static int
lexcmp_i(GEN x,GEN y)783 lexcmp_i(GEN x, GEN y)
784 {
785 const long tx = roughtype(x), ty = roughtype(y);
786 if (tx == ty)
787 switch(tx)
788 {
789 case t_MAT:
790 case t_VEC: return lexcmp_similar(x,y);
791 case t_VECSMALL: return vecsmall_lexcmp(x,y);
792 default: return gcmp(x,y);
793 }
794 if (tx == t_VECSMALL) return lexcmp_vecsmall_other(x,y,ty);
795 if (ty == t_VECSMALL) return -lexcmp_vecsmall_other(y,x,tx);
796
797 if (tx == t_INT) return lexcmp_scal_matvec(x,y); /*scalar*/
798 if (ty == t_INT) return -lexcmp_scal_matvec(y,x);
799
800 if (ty==t_MAT) return lexcmp_vec_mat(x,y);
801 return -lexcmp_vec_mat(y,x); /*tx==t_MAT*/
802 }
803 int
lexcmp(GEN x,GEN y)804 lexcmp(GEN x, GEN y)
805 {
806 pari_sp av = avma;
807 if (typ(x) == t_COMPLEX)
808 {
809 x = mkvec2(gel(x,1), gel(x,2));
810 if (typ(y) == t_COMPLEX) y = mkvec2(gel(y,1), gel(y,2));
811 else y = mkvec2(y, gen_0);
812 }
813 else if (typ(y) == t_COMPLEX)
814 {
815 x = mkvec2(x, gen_0);
816 y = mkvec2(gel(y,1), gel(y,2));
817 }
818 return gc_int(av, lexcmp_i(x, y));
819 }
820
821 /*****************************************************************/
822 /* */
823 /* EQUALITY */
824 /* returns 1 if x == y, 0 otherwise */
825 /* */
826 /*****************************************************************/
827 /* x,y t_POL */
828 static int
polidentical(GEN x,GEN y)829 polidentical(GEN x, GEN y)
830 {
831 long lx;
832 if (x[1] != y[1]) return 0;
833 lx = lg(x); if (lg(y) != lg(x)) return 0;
834 for (lx--; lx >= 2; lx--) if (!gidentical(gel(x,lx), gel(y,lx))) return 0;
835 return 1;
836 }
837 /* x,y t_SER */
838 static int
seridentical(GEN x,GEN y)839 seridentical(GEN x, GEN y) { return polidentical(x,y); }
840 /* typ(x) = typ(y) = t_VEC/COL/MAT */
841 static int
vecidentical(GEN x,GEN y)842 vecidentical(GEN x, GEN y)
843 {
844 long i;
845 if ((x[0] ^ y[0]) & (TYPBITS|LGBITS)) return 0;
846 for (i = lg(x)-1; i; i--)
847 if (! gidentical(gel(x,i),gel(y,i)) ) return 0;
848 return 1;
849 }
850 static int
identicalrr(GEN x,GEN y)851 identicalrr(GEN x, GEN y)
852 {
853 long i, lx = lg(x);
854 if (lg(y) != lx) return 0;
855 if (x[1] != y[1]) return 0;
856 i=2; while (i<lx && x[i]==y[i]) i++;
857 return (i == lx);
858 }
859
860 static int
closure_identical(GEN x,GEN y)861 closure_identical(GEN x, GEN y)
862 {
863 if (lg(x)!=lg(y) || x[1]!=y[1]) return 0;
864 if (!gidentical(gel(x,2),gel(y,2)) || !gidentical(gel(x,3),gel(y,3))
865 || !gidentical(gel(x,4),gel(y,4))) return 0;
866 if (lg(x)<8) return 1;
867 return gidentical(gel(x,7),gel(y,7));
868 }
869
870 static int
list_cmp(GEN x,GEN y,int cmp (GEN x,GEN y))871 list_cmp(GEN x, GEN y, int cmp(GEN x, GEN y))
872 {
873 int t = list_typ(x);
874 GEN vx, vy;
875 long lvx, lvy;
876 if (list_typ(y)!=t) return 0;
877 vx = list_data(x);
878 vy = list_data(y);
879 lvx = vx ? lg(vx): 1;
880 lvy = vy ? lg(vy): 1;
881 if (lvx==1 && lvy==1) return 1;
882 if (lvx != lvy) return 0;
883 switch (t)
884 {
885 case t_LIST_MAP:
886 {
887 pari_sp av = avma;
888 GEN mx = maptomat_shallow(x), my = maptomat_shallow(y);
889 int ret = gidentical(gel(mx, 1), gel(my, 1)) && cmp(gel(mx, 2), gel(my, 2));
890 return gc_bool(av, ret);
891 }
892 default:
893 return cmp(vx, vy);
894 }
895 }
896
897 int
gidentical(GEN x,GEN y)898 gidentical(GEN x, GEN y)
899 {
900 long tx;
901
902 if (x == y) return 1;
903 tx = typ(x); if (typ(y) != tx) return 0;
904 switch(tx)
905 {
906 case t_INT:
907 return equalii(x,y);
908
909 case t_REAL:
910 return identicalrr(x,y);
911
912 case t_FRAC: case t_INTMOD:
913 return equalii(gel(x,2), gel(y,2)) && equalii(gel(x,1), gel(y,1));
914
915 case t_COMPLEX:
916 return gidentical(gel(x,2),gel(y,2)) && gidentical(gel(x,1),gel(y,1));
917 case t_PADIC:
918 return valp(x) == valp(y)
919 && equalii(gel(x,2),gel(y,2))
920 && equalii(gel(x,3),gel(y,3))
921 && equalii(gel(x,4),gel(y,4));
922 case t_POLMOD:
923 return gidentical(gel(x,2),gel(y,2)) && polidentical(gel(x,1),gel(y,1));
924 case t_POL:
925 return polidentical(x,y);
926 case t_SER:
927 return seridentical(x,y);
928 case t_FFELT:
929 return FF_equal(x,y);
930
931 case t_QFR:
932 if (!identicalrr(gel(x,4),gel(y,4))) return 0; /* fall through */
933 case t_QFI:
934 return equalii(gel(x,1),gel(y,1))
935 && equalii(gel(x,2),gel(y,2))
936 && equalii(gel(x,3),gel(y,3));
937
938 case t_QUAD:
939 return ZX_equal(gel(x,1),gel(y,1))
940 && gidentical(gel(x,2),gel(y,2))
941 && gidentical(gel(x,3),gel(y,3));
942
943 case t_RFRAC:
944 return gidentical(gel(x,1),gel(y,1)) && gidentical(gel(x,2),gel(y,2));
945
946 case t_STR:
947 return !strcmp(GSTR(x),GSTR(y));
948 case t_VEC: case t_COL: case t_MAT:
949 return vecidentical(x,y);
950 case t_VECSMALL:
951 return zv_equal(x,y);
952 case t_CLOSURE:
953 return closure_identical(x,y);
954 case t_LIST:
955 return list_cmp(x, y, gidentical);
956 case t_INFINITY: return gidentical(gel(x,1),gel(y,1));
957 }
958 return 0;
959 }
960 /* x,y t_POL in the same variable */
961 static int
polequal(GEN x,GEN y)962 polequal(GEN x, GEN y)
963 {
964 long lx, ly;
965 if (signe(x) != signe(y)) return 0;
966 lx = lg(x); ly = lg(y);
967 while (lx > ly) if (!gequal0(gel(x,--lx))) return 0;
968 while (ly > lx) if (!gequal0(gel(y,--ly))) return 0;
969 for (lx--; lx >= 2; lx--) if (!gequal(gel(x,lx), gel(y,lx))) return 0;
970 return 1;
971 }
972
973 /* x,y t_SER in the same variable */
974 static int
serequal(GEN x,GEN y)975 serequal(GEN x, GEN y)
976 {
977 long LX, LY, lx, ly, vx, vy;
978 if (!signe(x) && !signe(y)) return 1;
979 lx = lg(x); vx = valp(x); LX = lx + vx;
980 ly = lg(y); vy = valp(y); LY = ly + vy;
981 if (LX > LY) lx = LY - vx; else ly = LX - vy;
982 while (lx >= 3 && ly >= 3)
983 if (!gequal(gel(x,--lx), gel(y,--ly))) return 0;
984 while(--ly >= 2) if (!gequal0(gel(y,ly))) return 0;
985 while(--lx >= 2) if (!gequal0(gel(x,lx))) return 0;
986 return 1;
987 }
988
989 /* typ(x) = typ(y) = t_VEC/COL/MAT */
990 static int
vecequal(GEN x,GEN y)991 vecequal(GEN x, GEN y)
992 {
993 long i;
994 if ((x[0] ^ y[0]) & (TYPBITS|LGBITS)) return 0;
995 for (i = lg(x)-1; i; i--)
996 if (! gequal(gel(x,i),gel(y,i)) ) return 0;
997 return 1;
998 }
999
1000 static int
gequal_try(GEN x,GEN y)1001 gequal_try(GEN x, GEN y)
1002 {
1003 int i;
1004 pari_CATCH(CATCH_ALL) {
1005 GEN E = pari_err_last();
1006 switch(err_get_num(E))
1007 {
1008 case e_STACK: case e_MEM: case e_ALARM:
1009 pari_err(0, E); /* rethrow */
1010 }
1011 return 0;
1012 } pari_TRY {
1013 i = gequal0(gadd(x, gneg_i(y)));
1014 } pari_ENDCATCH;
1015 return i;
1016 }
1017
1018 int
gequal(GEN x,GEN y)1019 gequal(GEN x, GEN y)
1020 {
1021 pari_sp av;
1022 long tx, ty;
1023 long i;
1024
1025 if (x == y) return 1;
1026 tx = typ(x);
1027 ty = typ(y);
1028 if (tx == ty)
1029 switch(tx)
1030 {
1031 case t_INT:
1032 return equalii(x,y);
1033
1034 case t_REAL:
1035 return equalrr(x,y);
1036
1037 case t_FRAC: case t_INTMOD:
1038 return equalii(gel(x,2), gel(y,2)) && equalii(gel(x,1), gel(y,1));
1039
1040 case t_COMPLEX:
1041 return gequal(gel(x,2),gel(y,2)) && gequal(gel(x,1),gel(y,1));
1042 case t_PADIC:
1043 if (!equalii(gel(x,2),gel(y,2))) return 0;
1044 av = avma; i = gequal0(gsub(x,y)); set_avma(av);
1045 return i;
1046 case t_POLMOD:
1047 if (varn(gel(x,1)) != varn(gel(y,1))) break;
1048 return gequal(gel(x,2),gel(y,2)) && RgX_equal_var(gel(x,1),gel(y,1));
1049 case t_POL:
1050 if (varn(x) != varn(y)) break;
1051 return polequal(x,y);
1052 case t_SER:
1053 if (varn(x) != varn(y)) break;
1054 return serequal(x,y);
1055
1056 case t_FFELT:
1057 return FF_equal(x,y);
1058
1059 case t_QFR:
1060 case t_QFI:
1061 return equalii(gel(x,1),gel(y,1))
1062 && equalii(gel(x,2),gel(y,2))
1063 && equalii(gel(x,3),gel(y,3));
1064
1065 case t_QUAD:
1066 return ZX_equal(gel(x,1),gel(y,1))
1067 && gequal(gel(x,2),gel(y,2))
1068 && gequal(gel(x,3),gel(y,3));
1069
1070 case t_RFRAC:
1071 {
1072 GEN a = gel(x,1), b = gel(x,2), c = gel(y,1), d = gel(y,2);
1073 if (gequal(b,d)) return gequal(a,c); /* simple case */
1074 av = avma;
1075 a = simplify_shallow(gmul(a,d));
1076 b = simplify_shallow(gmul(b,c));
1077 return gc_bool(av, gequal(a,b));
1078 }
1079
1080 case t_STR:
1081 return !strcmp(GSTR(x),GSTR(y));
1082 case t_VEC: case t_COL: case t_MAT:
1083 return vecequal(x,y);
1084 case t_VECSMALL:
1085 return zv_equal(x,y);
1086 case t_LIST:
1087 return list_cmp(x, y, gequal);
1088 case t_CLOSURE:
1089 return closure_identical(x,y);
1090 case t_INFINITY:
1091 return gequal(gel(x,1),gel(y,1));
1092 }
1093 if (is_noncalc_t(tx) || is_noncalc_t(ty)) return 0;
1094 if (tx == t_INT && !signe(x)) return gequal0(y);
1095 if (ty == t_INT && !signe(y)) return gequal0(x);
1096 (void)&av; av = avma; /* emulate volatile */
1097 return gc_bool(av, gequal_try(x, y));
1098 }
1099
1100 int
gequalsg(long s,GEN x)1101 gequalsg(long s, GEN x)
1102 { pari_sp av = avma; return gc_bool(av, gequal(stoi(s), x)); }
1103
1104 /* a and b are t_INT, t_FRAC, t_REAL or t_COMPLEX of those. Check whether
1105 * a-b is invertible */
1106 int
cx_approx_equal(GEN a,GEN b)1107 cx_approx_equal(GEN a, GEN b)
1108 {
1109 pari_sp av = avma;
1110 GEN d;
1111 if (a == b) return 1;
1112 d = gsub(a,b);
1113 return gc_bool(av, gequal0(d) || (typ(d)==t_COMPLEX && gequal0(cxnorm(d))));
1114 }
1115 /*******************************************************************/
1116 /* */
1117 /* VALUATION */
1118 /* p is either a t_INT or a t_POL. */
1119 /* returns the largest exponent of p dividing x when this makes */
1120 /* sense : error for types real, integermod and polymod if p does */
1121 /* not divide the modulus, q-adic if q!=p. */
1122 /* */
1123 /*******************************************************************/
1124
1125 static long
minval(GEN x,GEN p)1126 minval(GEN x, GEN p)
1127 {
1128 long i,k, val = LONG_MAX, lx = lg(x);
1129 for (i=lontyp[typ(x)]; i<lx; i++)
1130 {
1131 k = gvaluation(gel(x,i),p);
1132 if (k < val) val = k;
1133 }
1134 return val;
1135 }
1136
1137 static int
intdvd(GEN x,GEN y,GEN * z)1138 intdvd(GEN x, GEN y, GEN *z) { GEN r; *z = dvmdii(x,y,&r); return (r==gen_0); }
1139
1140 /* x t_FRAC, p t_INT, return v_p(x) */
1141 static long
frac_val(GEN x,GEN p)1142 frac_val(GEN x, GEN p) {
1143 long v = Z_pval(gel(x,2),p);
1144 if (v) return -v;
1145 return Z_pval(gel(x,1),p);
1146 }
1147 long
Q_pval(GEN x,GEN p)1148 Q_pval(GEN x, GEN p)
1149 {
1150 if (lgefint(p) == 3) return Q_lval(x, uel(p,2));
1151 return (typ(x)==t_INT)? Z_pval(x, p): frac_val(x, p);
1152 }
1153
1154 static long
frac_lval(GEN x,ulong p)1155 frac_lval(GEN x, ulong p) {
1156 long v = Z_lval(gel(x,2),p);
1157 if (v) return -v;
1158 return Z_lval(gel(x,1),p);
1159 }
1160 long
Q_lval(GEN x,ulong p)1161 Q_lval(GEN x, ulong p){return (typ(x)==t_INT)? Z_lval(x, p): frac_lval(x, p);}
1162
1163 long
Q_pvalrem(GEN x,GEN p,GEN * y)1164 Q_pvalrem(GEN x, GEN p, GEN *y)
1165 {
1166 GEN a, b;
1167 long v;
1168 if (lgefint(p) == 3) return Q_lvalrem(x, uel(p,2), y);
1169 if (typ(x) == t_INT) return Z_pvalrem(x, p, y);
1170 a = gel(x,1);
1171 b = gel(x,2);
1172 v = Z_pvalrem(b, p, &b);
1173 if (v) { *y = isint1(b)? a: mkfrac(a, b); return -v; }
1174 v = Z_pvalrem(a, p, &a);
1175 *y = mkfrac(a, b); return v;
1176 }
1177 long
Q_lvalrem(GEN x,ulong p,GEN * y)1178 Q_lvalrem(GEN x, ulong p, GEN *y)
1179 {
1180 GEN a, b;
1181 long v;
1182 if (typ(x) == t_INT) return Z_lvalrem(x, p, y);
1183 a = gel(x,1);
1184 b = gel(x,2);
1185 v = Z_lvalrem(b, p, &b);
1186 if (v) { *y = isint1(b)? a: mkfrac(a, b); return -v; }
1187 v = Z_lvalrem(a, p, &a);
1188 *y = mkfrac(a, b); return v;
1189 }
1190
1191 long
gvaluation(GEN x,GEN p)1192 gvaluation(GEN x, GEN p)
1193 {
1194 long tx = typ(x), tp = typ(p);
1195 pari_sp av;
1196
1197 switch(tp)
1198 {
1199 case t_INT:
1200 if (signe(p) && !is_pm1(p)) break;
1201 pari_err_DOMAIN("gvaluation", "p", "=", p, p);
1202 case t_POL:
1203 if (degpol(p) > 0) break;
1204 default:
1205 pari_err_DOMAIN("gvaluation", "p", "=", p, p);
1206 }
1207
1208 switch(tx)
1209 {
1210 case t_INT:
1211 if (!signe(x)) return LONG_MAX;
1212 if (tp == t_POL) return 0;
1213 return Z_pval(x,p);
1214
1215 case t_REAL:
1216 if (tp == t_POL) return 0;
1217 break;
1218
1219 case t_FFELT:
1220 if (tp == t_POL) return FF_equal0(x)? LONG_MAX: 0;
1221 break;
1222
1223 case t_INTMOD: {
1224 GEN a = gel(x,1), b = gel(x,2);
1225 long val;
1226 if (tp == t_POL) return signe(b)? 0: LONG_MAX;
1227 av = avma;
1228 if (!intdvd(a, p, &a)) break;
1229 if (!intdvd(b, p, &b)) return gc_long(av,0);
1230 val = 1; while (intdvd(a,p,&a) && intdvd(b,p,&b)) val++;
1231 return gc_long(av,val);
1232 }
1233
1234 case t_FRAC:
1235 if (tp == t_POL) return 0;
1236 return frac_val(x, p);
1237
1238 case t_PADIC:
1239 if (tp == t_POL) return 0;
1240 if (!equalii(p,gel(x,2))) break;
1241 return valp(x);
1242
1243 case t_POLMOD: {
1244 GEN a = gel(x,1), b = gel(x,2);
1245 long v, val;
1246 if (tp == t_INT) return gvaluation(b,p);
1247 v = varn(p);
1248 if (varn(a) != v) return 0;
1249 av = avma;
1250 a = RgX_divrem(a, p, ONLY_DIVIDES);
1251 if (!a) break;
1252 if (typ(b) != t_POL || varn(b) != v ||
1253 !(b = RgX_divrem(b, p, ONLY_DIVIDES)) ) return gc_long(av,0);
1254 val = 1;
1255 while ((a = RgX_divrem(a, p, ONLY_DIVIDES)) &&
1256 (b = RgX_divrem(b, p, ONLY_DIVIDES)) ) val++;
1257 return gc_long(av,val);
1258 }
1259 case t_POL: {
1260 if (tp == t_POL) {
1261 long vp = varn(p), vx = varn(x);
1262 if (vp == vx)
1263 {
1264 long val;
1265 if (RgX_is_monomial(p)) return RgX_val(x) / degpol(p);
1266 av = avma;
1267 for (val=0; ; val++)
1268 {
1269 x = RgX_divrem(x,p,ONLY_DIVIDES);
1270 if (!x) return gc_long(av,val);
1271 if (gc_needed(av,1))
1272 {
1273 if(DEBUGMEM>1) pari_warn(warnmem,"gvaluation");
1274 x = gerepilecopy(av, x);
1275 }
1276 }
1277 }
1278 if (varncmp(vx, vp) > 0) return 0;
1279 }
1280 return minval(x,p);
1281 }
1282
1283 case t_SER: {
1284 if (tp == t_POL) {
1285 long vp = varn(p), vx = varn(x);
1286 if (vp == vx)
1287 {
1288 long val = RgX_val(p);
1289 if (!val) pari_err_DOMAIN("gvaluation", "p", "=", p, p);
1290 return (long)(valp(x) / val);
1291 }
1292 if (varncmp(vx, vp) > 0) return 0;
1293 }
1294 return minval(x,p);
1295 }
1296
1297 case t_RFRAC:
1298 return gvaluation(gel(x,1),p) - gvaluation(gel(x,2),p);
1299
1300 case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
1301 return minval(x,p);
1302 }
1303 pari_err_OP("valuation", x,p);
1304 return 0; /* LCOV_EXCL_LINE */
1305 }
1306 GEN
gpvaluation(GEN x,GEN p)1307 gpvaluation(GEN x, GEN p)
1308 {
1309 long v = gvaluation(x,p);
1310 return v == LONG_MAX? mkoo(): stoi(v);
1311 }
1312
1313 /* x is nonzero */
1314 long
u_lvalrem(ulong x,ulong p,ulong * py)1315 u_lvalrem(ulong x, ulong p, ulong *py)
1316 {
1317 ulong vx;
1318 if (p == 2) { vx = vals(x); *py = x >> vx; return vx; }
1319 for(vx = 0;;)
1320 {
1321 if (x % p) { *py = x; return vx; }
1322 x /= p; /* gcc is smart enough to make a single div */
1323 vx++;
1324 }
1325 }
1326 long
u_lval(ulong x,ulong p)1327 u_lval(ulong x, ulong p)
1328 {
1329 ulong vx;
1330 if (p == 2) return vals(x);
1331 for(vx = 0;;)
1332 {
1333 if (x % p) return vx;
1334 x /= p; /* gcc is smart enough to make a single div */
1335 vx++;
1336 }
1337 }
1338
1339 long
z_lval(long s,ulong p)1340 z_lval(long s, ulong p) { return u_lval(labs(s), p); }
1341 long
z_lvalrem(long s,ulong p,long * py)1342 z_lvalrem(long s, ulong p, long *py)
1343 {
1344 long v;
1345 if (s < 0)
1346 {
1347 ulong u = (ulong)-s;
1348 v = u_lvalrem(u, p, &u);
1349 *py = -(long)u;
1350 }
1351 else
1352 {
1353 ulong u = (ulong)s;
1354 v = u_lvalrem(u, p, &u);
1355 *py = (long)u;
1356 }
1357 return v;
1358 }
1359 /* assume |p| > 1 */
1360 long
z_pval(long s,GEN p)1361 z_pval(long s, GEN p)
1362 {
1363 if (lgefint(p) > 3) return 0;
1364 return z_lval(s, uel(p,2));
1365 }
1366 /* assume |p| > 1 */
1367 long
z_pvalrem(long s,GEN p,long * py)1368 z_pvalrem(long s, GEN p, long *py)
1369 {
1370 if (lgefint(p) > 3) { *py = s; return 0; }
1371 return z_lvalrem(s, uel(p,2), py);
1372 }
1373
1374 /* return v_q(x) and set *py = x / q^v_q(x), using divide & conquer */
1375 static long
Z_pvalrem_DC(GEN x,GEN q,GEN * py)1376 Z_pvalrem_DC(GEN x, GEN q, GEN *py)
1377 {
1378 GEN r, z = dvmdii(x, q, &r);
1379 long v;
1380 if (r != gen_0) { *py = x; return 0; }
1381 if (2 * lgefint(q) <= lgefint(z)+3) /* avoid squaring if pointless */
1382 v = Z_pvalrem_DC(z, sqri(q), py) << 1;
1383 else
1384 { v = 0; *py = z; }
1385 z = dvmdii(*py, q, &r);
1386 if (r != gen_0) return v + 1;
1387 *py = z; return v + 2;
1388 }
1389
1390 static const long VAL_DC_THRESHOLD = 16;
1391
1392 long
Z_lval(GEN x,ulong p)1393 Z_lval(GEN x, ulong p)
1394 {
1395 long vx;
1396 pari_sp av;
1397 if (p == 2) return vali(x);
1398 if (lgefint(x) == 3) return u_lval(uel(x,2), p);
1399 av = avma;
1400 for(vx = 0;;)
1401 {
1402 ulong r;
1403 GEN q = absdiviu_rem(x, p, &r);
1404 if (r) break;
1405 vx++; x = q;
1406 if (vx == VAL_DC_THRESHOLD) {
1407 if (p == 1) pari_err_DOMAIN("Z_lval", "p", "=", gen_1, gen_1);
1408 vx += Z_pvalrem_DC(x, sqru(p), &x) << 1;
1409 q = absdiviu_rem(x, p, &r); if (!r) vx++;
1410 break;
1411 }
1412 }
1413 return gc_long(av,vx);
1414 }
1415 long
Z_lvalrem(GEN x,ulong p,GEN * py)1416 Z_lvalrem(GEN x, ulong p, GEN *py)
1417 {
1418 long vx, sx;
1419 pari_sp av;
1420 if (p == 2) { vx = vali(x); *py = shifti(x, -vx); return vx; }
1421 if (lgefint(x) == 3) {
1422 ulong u;
1423 vx = u_lvalrem(uel(x,2), p, &u);
1424 *py = signe(x) < 0? utoineg(u): utoipos(u);
1425 return vx;
1426 }
1427 av = avma; (void)new_chunk(lgefint(x));
1428 sx = signe(x);
1429 for(vx = 0;;)
1430 {
1431 ulong r;
1432 GEN q = absdiviu_rem(x, p, &r);
1433 if (r) break;
1434 vx++; x = q;
1435 if (vx == VAL_DC_THRESHOLD) {
1436 if (p == 1) pari_err_DOMAIN("Z_lvalrem", "p", "=", gen_1, gen_1);
1437 vx += Z_pvalrem_DC(x, sqru(p), &x) << 1;
1438 q = absdiviu_rem(x, p, &r); if (!r) { vx++; x = q; }
1439 break;
1440 }
1441 }
1442 set_avma(av); *py = icopy(x); setsigne(*py, sx); return vx;
1443 }
1444
1445 /* Is |q| <= p ? */
1446 static int
isless_iu(GEN q,ulong p)1447 isless_iu(GEN q, ulong p) {
1448 long l = lgefint(q);
1449 return l==2 || (l == 3 && uel(q,2) <= p);
1450 }
1451
1452 long
u_lvalrem_stop(ulong * n,ulong p,int * stop)1453 u_lvalrem_stop(ulong *n, ulong p, int *stop)
1454 {
1455 ulong N = *n, q = N / p, r = N % p; /* gcc makes a single div */
1456 long v = 0;
1457 if (!r)
1458 {
1459 do { v++; N = q; q = N / p; r = N % p; } while (!r);
1460 *n = N;
1461 }
1462 *stop = q <= p; return v;
1463 }
1464 /* Assume n > 0. Return v_p(n), set *n := n/p^v_p(n). Set 'stop' if now
1465 * n < p^2 [implies n prime if no prime < p divides n] */
1466 long
Z_lvalrem_stop(GEN * n,ulong p,int * stop)1467 Z_lvalrem_stop(GEN *n, ulong p, int *stop)
1468 {
1469 pari_sp av;
1470 long v;
1471 ulong r;
1472 GEN N, q;
1473
1474 if (lgefint(*n) == 3)
1475 {
1476 r = (*n)[2];
1477 v = u_lvalrem_stop(&r, p, stop);
1478 if (v) *n = utoipos(r);
1479 return v;
1480 }
1481 av = avma; v = 0; q = absdiviu_rem(*n, p, &r);
1482 if (r) set_avma(av);
1483 else
1484 {
1485 do {
1486 v++; N = q;
1487 if (v == VAL_DC_THRESHOLD)
1488 {
1489 v += Z_pvalrem_DC(N,sqru(p),&N) << 1;
1490 q = absdiviu_rem(N, p, &r); if (!r) { v++; N = q; }
1491 break;
1492 }
1493 q = absdiviu_rem(N, p, &r);
1494 } while (!r);
1495 *n = N;
1496 }
1497 *stop = isless_iu(q,p); return v;
1498 }
1499
1500 /* x is a nonzero integer, |p| > 1 */
1501 long
Z_pvalrem(GEN x,GEN p,GEN * py)1502 Z_pvalrem(GEN x, GEN p, GEN *py)
1503 {
1504 long vx;
1505 pari_sp av;
1506
1507 if (lgefint(p) == 3) return Z_lvalrem(x, uel(p,2), py);
1508 if (lgefint(x) == 3) { *py = icopy(x); return 0; }
1509 av = avma; vx = 0; (void)new_chunk(lgefint(x));
1510 for(;;)
1511 {
1512 GEN r, q = dvmdii(x,p,&r);
1513 if (r != gen_0) { set_avma(av); *py = icopy(x); return vx; }
1514 vx++; x = q;
1515 }
1516 }
1517 long
u_pvalrem(ulong x,GEN p,ulong * py)1518 u_pvalrem(ulong x, GEN p, ulong *py)
1519 {
1520 if (lgefint(p) == 3) return u_lvalrem(x, uel(p,2), py);
1521 *py = x; return 0;
1522 }
1523 long
u_pval(ulong x,GEN p)1524 u_pval(ulong x, GEN p)
1525 {
1526 if (lgefint(p) == 3) return u_lval(x, uel(p,2));
1527 return 0;
1528 }
1529 long
Z_pval(GEN x,GEN p)1530 Z_pval(GEN x, GEN p) {
1531 long vx;
1532 pari_sp av;
1533
1534 if (lgefint(p) == 3) return Z_lval(x, uel(p,2));
1535 if (lgefint(x) == 3) return 0;
1536 av = avma; vx = 0;
1537 for(;;)
1538 {
1539 GEN r, q = dvmdii(x,p,&r);
1540 if (r != gen_0) return gc_long(av,vx);
1541 vx++; x = q;
1542 }
1543 }
1544
1545 /* return v_p(n!) = [n/p] + [n/p^2] + ... */
1546 long
factorial_lval(ulong n,ulong p)1547 factorial_lval(ulong n, ulong p)
1548 {
1549 ulong q, v;
1550 if (p == 2) return n - hammingl(n);
1551 q = p; v = 0;
1552 do { v += n/q; q *= p; } while (n >= q);
1553 return (long)v;
1554 }
1555
1556 /********** Same for "containers" ZX / ZV / ZC **********/
1557
1558 /* If the t_INT q divides the ZX/ZV x, return the quotient. Otherwise NULL.
1559 * Stack clean; assumes lg(x) > 1 */
1560 static GEN
gen_Z_divides(GEN x,GEN q,long imin)1561 gen_Z_divides(GEN x, GEN q, long imin)
1562 {
1563 long i, l;
1564 GEN y = cgetg_copy(x, &l);
1565
1566 y[1] = x[1]; /* Needed for ZX; no-op if ZV, overwritten in first iteration */
1567 for (i = imin; i < l; i++)
1568 {
1569 GEN r, xi = gel(x,i);
1570 if (!signe(xi)) { gel(y,i) = xi; continue; }
1571 gel(y,i) = dvmdii(xi, q, &r);
1572 if (r != gen_0) { set_avma((pari_sp)(y+l)); return NULL; }
1573 }
1574 return y;
1575 }
1576 /* If q divides the ZX/ZV x, return the quotient. Otherwise NULL.
1577 * Stack clean; assumes lg(x) > 1 */
1578 static GEN
gen_z_divides(GEN x,ulong q,long imin)1579 gen_z_divides(GEN x, ulong q, long imin)
1580 {
1581 long i, l;
1582 GEN y = cgetg_copy(x, &l);
1583
1584 y[1] = x[1]; /* Needed for ZX; no-op if ZV, overwritten in first iteration */
1585 for (i = imin; i < l; i++)
1586 {
1587 ulong r;
1588 GEN xi = gel(x,i);
1589 if (!signe(xi)) { gel(y,i) = xi; continue; }
1590 gel(y,i) = absdiviu_rem(xi, q, &r);
1591 if (r) { set_avma((pari_sp)(y+l)); return NULL; }
1592 affectsign_safe(xi, &gel(y,i));
1593 }
1594 return y;
1595 }
1596
1597 /* return v_q(x) and set *py = x / q^v_q(x), using divide & conquer */
1598 static long
gen_pvalrem_DC(GEN x,GEN q,GEN * py,long imin)1599 gen_pvalrem_DC(GEN x, GEN q, GEN *py, long imin)
1600 {
1601
1602 pari_sp av = avma;
1603 long v, i, l, lz = LONG_MAX;
1604 GEN y = cgetg_copy(x, &l);
1605
1606 y[1] = x[1];
1607 for (i = imin; i < l; i++)
1608 {
1609 GEN r, xi = gel(x,i);
1610 if (!signe(xi)) { gel(y,i) = xi; continue; }
1611 gel(y,i) = dvmdii(xi, q, &r);
1612 if (r != gen_0) { *py = x; return gc_long(av,0); }
1613 lz = minss(lz, lgefint(gel(y,i)));
1614 }
1615 if (2 * lgefint(q) <= lz+3) /* avoid squaring if pointless */
1616 v = gen_pvalrem_DC(y, sqri(q), py, imin) << 1;
1617 else
1618 { v = 0; *py = y; }
1619
1620 y = gen_Z_divides(*py, q, imin);
1621 if (!y) return v+1;
1622 *py = y; return v+2;
1623 }
1624
1625 static long
gen_2val(GEN x,long imin)1626 gen_2val(GEN x, long imin)
1627 {
1628 long i, lx = lg(x), v = LONG_MAX;
1629 for (i = imin; i < lx; i++)
1630 {
1631 GEN c = gel(x,i);
1632 long w;
1633 if (!signe(c)) continue;
1634 w = vali(c);
1635 if (w < v) { v = w; if (!v) break; }
1636 }
1637 return v;
1638 }
1639 static long
gen_lval(GEN x,ulong p,long imin)1640 gen_lval(GEN x, ulong p, long imin)
1641 {
1642 long i, lx, v;
1643 pari_sp av;
1644 GEN y;
1645 if (p == 2) return gen_2val(x, imin);
1646 av = avma;
1647 lx = lg(x); y = leafcopy(x);
1648 for(v = 0;; v++)
1649 for (i = imin; i < lx; i++)
1650 {
1651 ulong r;
1652 gel(y,i) = absdiviu_rem(gel(y,i), p, &r);
1653 if (r) return gc_long(av,v);
1654 }
1655 }
1656 long
ZX_lval(GEN x,ulong p)1657 ZX_lval(GEN x, ulong p) { return gen_lval(x, p, 2); }
1658 long
ZV_lval(GEN x,ulong p)1659 ZV_lval(GEN x, ulong p) { return gen_lval(x, p, 1); }
1660
1661 long
zx_lval(GEN f,long p)1662 zx_lval(GEN f, long p)
1663 {
1664 long i, l = lg(f), x = LONG_MAX;
1665 for(i=2; i<l; i++)
1666 {
1667 long y;
1668 if (f[i] == 0) continue;
1669 y = z_lval(f[i], p);
1670 if (y < x) { x = y; if (x == 0) return x; }
1671 }
1672 return x;
1673 }
1674
1675 static long
gen_pval(GEN x,GEN p,long imin)1676 gen_pval(GEN x, GEN p, long imin)
1677 {
1678 long i, lx, v;
1679 pari_sp av;
1680 GEN y;
1681 if (lgefint(p) == 3) return gen_lval(x, p[2], imin);
1682 av = avma;
1683 lx = lg(x); y = leafcopy(x);
1684 for(v = 0;; v++)
1685 {
1686 if (v == VAL_DC_THRESHOLD)
1687 {
1688 if (is_pm1(p)) pari_err_DOMAIN("gen_pval", "p", "=", p, p);
1689 v += gen_pvalrem_DC(y, p, &y, imin);
1690 return gc_long(av,v);
1691 }
1692
1693 for (i = imin; i < lx; i++)
1694 {
1695 GEN r; gel(y,i) = dvmdii(gel(y,i), p, &r);
1696 if (r != gen_0) return gc_long(av,v);
1697 }
1698 }
1699 }
1700 long
ZX_pval(GEN x,GEN p)1701 ZX_pval(GEN x, GEN p) { return gen_pval(x, p, 2); }
1702 long
ZV_pval(GEN x,GEN p)1703 ZV_pval(GEN x, GEN p) { return gen_pval(x, p, 1); }
1704 /* v = 0 (mod p) */
1705 int
ZV_Z_dvd(GEN v,GEN p)1706 ZV_Z_dvd(GEN v, GEN p)
1707 {
1708 pari_sp av = avma;
1709 long i, l = lg(v);
1710 for (i=1; i<l; i++)
1711 if (!dvdii(gel(v,i), p)) return gc_long(av,0);
1712 return gc_long(av,1);
1713 }
1714
1715 static long
gen_2valrem(GEN x,GEN * px,long imin)1716 gen_2valrem(GEN x, GEN *px, long imin)
1717 {
1718 long i, lx = lg(x), v = LONG_MAX;
1719 GEN z;
1720 for (i = imin; i < lx; i++)
1721 {
1722 GEN c = gel(x,i);
1723 long w;
1724 if (!signe(c)) continue;
1725 w = vali(c);
1726 if (w < v) {
1727 v = w;
1728 if (!v) { *px = x; return 0; } /* early abort */
1729 }
1730 }
1731 z = cgetg_copy(x, &lx); z[1] = x[1];
1732 for (i=imin; i<lx; i++) gel(z,i) = shifti(gel(x,i), -v);
1733 *px = z; return v;
1734 }
1735 static long
gen_lvalrem(GEN x,ulong p,GEN * px,long imin)1736 gen_lvalrem(GEN x, ulong p, GEN *px, long imin)
1737 {
1738 long i, lx, v;
1739 GEN y;
1740 if (p == 2) return gen_2valrem(x, px, imin);
1741 y = cgetg_copy(x, &lx);
1742 y[1] = x[1];
1743 x = leafcopy(x);
1744 for(v = 0;; v++)
1745 {
1746 if (v == VAL_DC_THRESHOLD)
1747 {
1748 if (p == 1) pari_err_DOMAIN("gen_lvalrem", "p", "=", gen_1, gen_1);
1749 v += gen_pvalrem_DC(x, sqru(p), px, imin) << 1;
1750 x = gen_z_divides(*px, p, imin);
1751 if (x) { *px = x; v++; }
1752 return v;
1753 }
1754
1755 for (i = imin; i < lx; i++)
1756 {
1757 ulong r; gel(y,i) = absdiviu_rem(gel(x,i), p, &r);
1758 if (r) { *px = x; return v; }
1759 affectsign_safe(gel(x,i), &gel(y,i));
1760 }
1761 swap(x, y);
1762 }
1763 }
1764 long
ZX_lvalrem(GEN x,ulong p,GEN * px)1765 ZX_lvalrem(GEN x, ulong p, GEN *px) { return gen_lvalrem(x,p,px, 2); }
1766 long
ZV_lvalrem(GEN x,ulong p,GEN * px)1767 ZV_lvalrem(GEN x, ulong p, GEN *px) { return gen_lvalrem(x,p,px, 1); }
1768
1769 static long
gen_pvalrem(GEN x,GEN p,GEN * px,long imin)1770 gen_pvalrem(GEN x, GEN p, GEN *px, long imin)
1771 {
1772 long i, lx, v;
1773 GEN y;
1774 if (lgefint(p) == 3) return gen_lvalrem(x, p[2], px, imin);
1775 y = cgetg_copy(x, &lx);
1776 y[1] = x[1];
1777 x = leafcopy(x);
1778 for(v = 0;; v++)
1779 {
1780 if (v == VAL_DC_THRESHOLD)
1781 {
1782 if (is_pm1(p)) pari_err_DOMAIN("gen_pvalrem", "p", "=", p, p);
1783 return v + gen_pvalrem_DC(x, p, px, imin);
1784 }
1785
1786 for (i = imin; i < lx; i++)
1787 {
1788 GEN r; gel(y,i) = dvmdii(gel(x,i), p, &r);
1789 if (r != gen_0) { *px = x; return v; }
1790 }
1791 swap(x, y);
1792 }
1793 }
1794 long
ZX_pvalrem(GEN x,GEN p,GEN * px)1795 ZX_pvalrem(GEN x, GEN p, GEN *px) { return gen_pvalrem(x,p,px, 2); }
1796 long
ZV_pvalrem(GEN x,GEN p,GEN * px)1797 ZV_pvalrem(GEN x, GEN p, GEN *px) { return gen_pvalrem(x,p,px, 1); }
1798
1799 /*******************************************************************/
1800 /* */
1801 /* NEGATION: Create -x */
1802 /* */
1803 /*******************************************************************/
1804
1805 GEN
gneg(GEN x)1806 gneg(GEN x)
1807 {
1808 long lx, i;
1809 GEN y;
1810
1811 switch(typ(x))
1812 {
1813 case t_INT:
1814 return signe(x)? negi(x): gen_0;
1815 case t_REAL:
1816 return mpneg(x);
1817
1818 case t_INTMOD: y=cgetg(3,t_INTMOD);
1819 gel(y,1) = icopy(gel(x,1));
1820 gel(y,2) = signe(gel(x,2))? subii(gel(y,1),gel(x,2)): gen_0;
1821 break;
1822
1823 case t_FRAC:
1824 y = cgetg(3, t_FRAC);
1825 gel(y,1) = negi(gel(x,1));
1826 gel(y,2) = icopy(gel(x,2)); break;
1827
1828 case t_COMPLEX:
1829 y=cgetg(3, t_COMPLEX);
1830 gel(y,1) = gneg(gel(x,1));
1831 gel(y,2) = gneg(gel(x,2));
1832 break;
1833
1834 case t_POLMOD:
1835 retmkpolmod(gneg(gel(x,2)), RgX_copy(gel(x,1)));
1836
1837 case t_RFRAC:
1838 y = cgetg(3, t_RFRAC);
1839 gel(y,1) = gneg(gel(x,1));
1840 gel(y,2) = RgX_copy(gel(x,2)); break;
1841
1842 case t_PADIC:
1843 if (!signe(gel(x,4))) return gcopy(x);
1844 y = cgetg(5, t_PADIC);
1845 y[1] = x[1];
1846 gel(y,2) = icopy(gel(x,2));
1847 gel(y,3) = icopy(gel(x,3));
1848 gel(y,4) = subii(gel(x,3),gel(x,4));
1849 break;
1850
1851 case t_QUAD:
1852 y=cgetg(4,t_QUAD);
1853 gel(y,1) = ZX_copy(gel(x,1));
1854 gel(y,2) = gneg(gel(x,2));
1855 gel(y,3) = gneg(gel(x,3)); break;
1856
1857 case t_FFELT: return FF_neg(x);
1858 case t_POL: return RgX_neg(x);
1859 case t_SER:
1860 y = cgetg_copy(x, &lx); y[1] = x[1];
1861 for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
1862 break;
1863 case t_VEC: return RgV_neg(x);
1864 case t_COL: return RgC_neg(x);
1865 case t_MAT: return RgM_neg(x);
1866 case t_INFINITY: return inf_get_sign(x) == 1? mkmoo(): mkoo();
1867 default:
1868 pari_err_TYPE("gneg",x);
1869 return NULL; /* LCOV_EXCL_LINE */
1870 }
1871 return y;
1872 }
1873
1874 GEN
gneg_i(GEN x)1875 gneg_i(GEN x)
1876 {
1877 long lx, i;
1878 GEN y;
1879
1880 switch(typ(x))
1881 {
1882 case t_INT:
1883 return signe(x)? negi(x): gen_0;
1884 case t_REAL:
1885 return mpneg(x);
1886
1887 case t_INTMOD: y=cgetg(3,t_INTMOD);
1888 gel(y,1) = gel(x,1);
1889 gel(y,2) = signe(gel(x,2))? subii(gel(y,1),gel(x,2)): gen_0;
1890 break;
1891
1892 case t_FRAC:
1893 y = cgetg(3, t_FRAC);
1894 gel(y,1) = negi(gel(x,1));
1895 gel(y,2) = gel(x,2); break;
1896
1897 case t_COMPLEX:
1898 y = cgetg(3, t_COMPLEX);
1899 gel(y,1) = gneg_i(gel(x,1));
1900 gel(y,2) = gneg_i(gel(x,2)); break;
1901
1902 case t_PADIC: y = cgetg(5,t_PADIC);
1903 y[1] = x[1];
1904 gel(y,2) = gel(x,2);
1905 gel(y,3) = gel(x,3);
1906 gel(y,4) = signe(gel(x,4))? subii(gel(x,3),gel(x,4)): gen_0; break;
1907
1908 case t_POLMOD:
1909 retmkpolmod(gneg_i(gel(x,2)), RgX_copy(gel(x,1)));
1910
1911 case t_FFELT: return FF_neg_i(x);
1912
1913 case t_QUAD: y=cgetg(4,t_QUAD);
1914 gel(y,1) = gel(x,1);
1915 gel(y,2) = gneg_i(gel(x,2));
1916 gel(y,3) = gneg_i(gel(x,3)); break;
1917
1918 case t_VEC: case t_COL: case t_MAT:
1919 y = cgetg_copy(x, &lx);
1920 for (i=1; i<lx; i++) gel(y,i) = gneg_i(gel(x,i));
1921 break;
1922
1923 case t_POL: case t_SER:
1924 y = cgetg_copy(x, &lx); y[1]=x[1];
1925 for (i=2; i<lx; i++) gel(y,i) = gneg_i(gel(x,i));
1926 break;
1927
1928 case t_RFRAC:
1929 y = cgetg(3, t_RFRAC);
1930 gel(y,1) = gneg_i(gel(x,1));
1931 gel(y,2) = gel(x,2); break;
1932
1933 default:
1934 pari_err_TYPE("gneg_i",x);
1935 return NULL; /* LCOV_EXCL_LINE */
1936 }
1937 return y;
1938 }
1939
1940 /******************************************************************/
1941 /* */
1942 /* ABSOLUTE VALUE */
1943 /* Create abs(x) if x is integer, real, fraction or complex. */
1944 /* Error otherwise. */
1945 /* */
1946 /******************************************************************/
1947 static int
is_negative(GEN x)1948 is_negative(GEN x) {
1949 switch(typ(x))
1950 {
1951 case t_INT: case t_REAL:
1952 return (signe(x) < 0);
1953 case t_FRAC:
1954 return (signe(gel(x,1)) < 0);
1955 }
1956 return 0;
1957 }
1958
1959 GEN
gabs(GEN x,long prec)1960 gabs(GEN x, long prec)
1961 {
1962 long lx, i;
1963 pari_sp av, tetpil;
1964 GEN y,p1;
1965
1966 switch(typ(x))
1967 {
1968 case t_INT: case t_REAL:
1969 return mpabs(x);
1970
1971 case t_FRAC:
1972 return absfrac(x);
1973
1974 case t_COMPLEX:
1975 av=avma; p1=cxnorm(x);
1976 switch(typ(p1))
1977 {
1978 case t_INT:
1979 if (!Z_issquareall(p1, &y)) break;
1980 return gerepileupto(av, y);
1981 case t_FRAC: {
1982 GEN a,b;
1983 if (!Z_issquareall(gel(p1,1), &a)) break;
1984 if (!Z_issquareall(gel(p1,2), &b)) break;
1985 return gerepileupto(av, gdiv(a,b));
1986 }
1987 }
1988 tetpil=avma;
1989 return gerepile(av,tetpil,gsqrt(p1,prec));
1990
1991 case t_QUAD:
1992 av = avma;
1993 return gerepileuptoleaf(av, gabs(quadtofp(x, prec), prec));
1994
1995 case t_POL:
1996 lx = lg(x); if (lx<=2) return RgX_copy(x);
1997 return is_negative(gel(x,lx-1))? gneg(x): RgX_copy(x);
1998
1999 case t_SER:
2000 if (!signe(x)) pari_err_DOMAIN("abs", "argument", "=", gen_0, x);
2001 if (valp(x)) pari_err_DOMAIN("abs", "series valuation", "!=", gen_0, x);
2002 return is_negative(gel(x,2))? gneg(x): gcopy(x);
2003
2004 case t_VEC: case t_COL: case t_MAT:
2005 y = cgetg_copy(x, &lx);
2006 for (i=1; i<lx; i++) gel(y,i) = gabs(gel(x,i),prec);
2007 return y;
2008 }
2009 pari_err_TYPE("gabs",x);
2010 return NULL; /* LCOV_EXCL_LINE */
2011 }
2012
2013 GEN
gmax(GEN x,GEN y)2014 gmax(GEN x, GEN y) { return gcopy(gmax_shallow(x,y)); }
2015 GEN
gmaxgs(GEN x,long s)2016 gmaxgs(GEN x, long s) { return (gcmpsg(s,x)>=0)? stoi(s): gcopy(x); }
2017
2018 GEN
gmin(GEN x,GEN y)2019 gmin(GEN x, GEN y) { return gcopy(gmin_shallow(x,y)); }
2020 GEN
gmings(GEN x,long s)2021 gmings(GEN x, long s) { return (gcmpsg(s,x)>0)? gcopy(x): stoi(s); }
2022
2023 long
vecindexmax(GEN x)2024 vecindexmax(GEN x)
2025 {
2026 long lx = lg(x), i0, i;
2027 GEN s;
2028
2029 if (lx==1) pari_err_DOMAIN("vecindexmax", "empty argument", "=", x,x);
2030 switch(typ(x))
2031 {
2032 case t_VEC: case t_COL:
2033 s = gel(x,i0=1);
2034 for (i=2; i<lx; i++)
2035 if (gcmp(gel(x,i),s) > 0) s = gel(x,i0=i);
2036 return i0;
2037 case t_VECSMALL:
2038 return vecsmall_indexmax(x);
2039 default: pari_err_TYPE("vecindexmax",x);
2040 }
2041 /* LCOV_EXCL_LINE */
2042 return 0;
2043 }
2044 long
vecindexmin(GEN x)2045 vecindexmin(GEN x)
2046 {
2047 long lx = lg(x), i0, i;
2048 GEN s;
2049
2050 if (lx==1) pari_err_DOMAIN("vecindexmin", "empty argument", "=", x,x);
2051 switch(typ(x))
2052 {
2053 case t_VEC: case t_COL:
2054 s = gel(x,i0=1);
2055 for (i=2; i<lx; i++)
2056 if (gcmp(gel(x,i),s) < 0) s = gel(x,i0=i);
2057 return i0;
2058 case t_VECSMALL:
2059 return vecsmall_indexmin(x);
2060 default: pari_err_TYPE("vecindexmin",x);
2061 }
2062 /* LCOV_EXCL_LINE */
2063 return 0;
2064 }
2065
2066 GEN
vecmax0(GEN x,GEN * pi)2067 vecmax0(GEN x, GEN *pi)
2068 {
2069 long i, lx = lg(x), tx = typ(x);
2070 if (!is_matvec_t(tx) && tx != t_VECSMALL) return gcopy(x);
2071 if (lx==1) pari_err_DOMAIN("vecmax", "empty argument", "=", x,x);
2072 switch(typ(x))
2073 {
2074 case t_VEC: case t_COL:
2075 i = vecindexmax(x); if (pi) *pi = utoipos(i);
2076 return gcopy(gel(x,i));
2077 case t_MAT: {
2078 long j, i0 = 1, j0 = 1, lx2 = lgcols(x);
2079 GEN s;
2080 if (lx2 == 1) pari_err_DOMAIN("vecmax", "empty argument", "=", x,x);
2081 s = gcoeff(x,i0,j0); i = 2;
2082 for (j=1; j<lx; j++,i=1)
2083 {
2084 GEN c = gel(x,j);
2085 for (; i<lx2; i++)
2086 if (gcmp(gel(c,i),s) > 0) { s = gel(c,i); j0=j; i0=i; }
2087 }
2088 if (pi) *pi = mkvec2(utoipos(i0), utoipos(j0));
2089 return gcopy(s);
2090 }
2091 case t_VECSMALL:
2092 i = vecsmall_indexmax(x); if (pi) *pi = utoipos(i);
2093 return stoi(x[i]);
2094 }
2095 return NULL;/*LCOV_EXCL_LINE*/
2096 }
2097 GEN
vecmin0(GEN x,GEN * pi)2098 vecmin0(GEN x, GEN *pi)
2099 {
2100 long i, lx = lg(x), tx = typ(x);
2101 if (!is_matvec_t(tx) && tx != t_VECSMALL) return gcopy(x);
2102 if (lx==1) pari_err_DOMAIN("vecmin", "empty argument", "=", x,x);
2103 switch(typ(x))
2104 {
2105 case t_VEC: case t_COL:
2106 i = vecindexmin(x); if (pi) *pi = utoipos(i);
2107 return gcopy(gel(x,i));
2108 case t_MAT: {
2109 long j, i0 = 1, j0 = 1, lx2 = lgcols(x);
2110 GEN s;
2111 if (lx2 == 1) pari_err_DOMAIN("vecmin", "empty argument", "=", x,x);
2112 s = gcoeff(x,i0,j0); i = 2;
2113 for (j=1; j<lx; j++,i=1)
2114 {
2115 GEN c = gel(x,j);
2116 for (; i<lx2; i++)
2117 if (gcmp(gel(c,i),s) < 0) { s = gel(c,i); j0=j; i0=i; }
2118 }
2119 if (pi) *pi = mkvec2(utoipos(i0), utoipos(j0));
2120 return gcopy(s);
2121 }
2122 case t_VECSMALL:
2123 i = vecsmall_indexmin(x); if (pi) *pi = utoipos(i);
2124 return stoi(x[i]);
2125 }
2126 return NULL;/*LCOV_EXCL_LINE*/
2127 }
2128
2129 GEN
vecmax(GEN x)2130 vecmax(GEN x) { return vecmax0(x, NULL); }
2131 GEN
vecmin(GEN x)2132 vecmin(GEN x) { return vecmin0(x, NULL); }
2133
2134 /*******************************************************************/
2135 /* */
2136 /* AFFECT long --> GEN */
2137 /* affect long s to GEN x. Useful for initialization. */
2138 /* */
2139 /*******************************************************************/
2140
2141 static void
padicaff0(GEN x)2142 padicaff0(GEN x)
2143 {
2144 if (signe(gel(x,4)))
2145 {
2146 x[1] = evalvalp(valp(x)+precp(x));
2147 affsi(0,gel(x,4));
2148 }
2149 }
2150
2151 void
gaffsg(long s,GEN x)2152 gaffsg(long s, GEN x)
2153 {
2154 switch(typ(x))
2155 {
2156 case t_INT: affsi(s,x); break;
2157 case t_REAL: affsr(s,x); break;
2158 case t_INTMOD: modsiz(s,gel(x,1),gel(x,2)); break;
2159 case t_FRAC: affsi(s,gel(x,1)); affsi(1,gel(x,2)); break;
2160 case t_COMPLEX: gaffsg(s,gel(x,1)); gaffsg(0,gel(x,2)); break;
2161 case t_PADIC: {
2162 long vx;
2163 GEN y;
2164 if (!s) { padicaff0(x); break; }
2165 vx = Z_pvalrem(stoi(s), gel(x,2), &y);
2166 setvalp(x,vx); modiiz(y,gel(x,3),gel(x,4));
2167 break;
2168 }
2169 case t_QUAD: gaffsg(s,gel(x,2)); gaffsg(0,gel(x,3)); break;
2170 default: pari_err_TYPE2("=",stoi(s),x);
2171 }
2172 }
2173
2174 /*******************************************************************/
2175 /* */
2176 /* GENERIC AFFECTATION */
2177 /* Affect the content of x to y, whenever possible */
2178 /* */
2179 /*******************************************************************/
2180 /* x PADIC, Y INT, return lift(x * Mod(1,Y)) */
2181 GEN
padic_to_Fp(GEN x,GEN Y)2182 padic_to_Fp(GEN x, GEN Y) {
2183 pari_sp av = avma;
2184 GEN p = gel(x,2), z;
2185 long vy, vx = valp(x);
2186 if (!signe(Y)) pari_err_INV("padic_to_Fp",Y);
2187 vy = Z_pvalrem(Y,p, &z);
2188 if (vx < 0 || !gequal1(z)) pari_err_OP("",x, mkintmod(gen_1,Y));
2189 if (vx >= vy) { set_avma(av); return gen_0; }
2190 z = gel(x,4);
2191 if (!signe(z) || vy > vx + precp(x)) pari_err_OP("",x, mkintmod(gen_1,Y));
2192 if (vx) z = mulii(z, powiu(p,vx));
2193 return gerepileuptoint(av, remii(z, Y));
2194 }
2195 ulong
padic_to_Fl(GEN x,ulong Y)2196 padic_to_Fl(GEN x, ulong Y) {
2197 GEN p = gel(x,2);
2198 ulong u, z;
2199 long vy, vx = valp(x);
2200 vy = u_pvalrem(Y,p, &u);
2201 if (vx < 0 || u != 1) pari_err_OP("",x, mkintmodu(1,Y));
2202 /* Y = p^vy */
2203 if (vx >= vy) return 0;
2204 z = umodiu(gel(x,4), Y);
2205 if (!z || vy > vx + precp(x)) pari_err_OP("",x, mkintmodu(1,Y));
2206 if (vx) {
2207 ulong pp = p[2];
2208 z = Fl_mul(z, upowuu(pp,vx), Y); /* p^vx < p^vy = Y */
2209 }
2210 return z;
2211 }
2212
2213 static void
croak(const char * s)2214 croak(const char *s) {
2215 char *t;
2216 t = stack_sprintf("gaffect [overwriting universal object: %s]",s);
2217 pari_err_BUG(t);
2218 }
2219
2220 void
gaffect(GEN x,GEN y)2221 gaffect(GEN x, GEN y)
2222 {
2223 long vx, i, lx, ly, tx = typ(x), ty = typ(y);
2224 pari_sp av;
2225 GEN p1, num, den;
2226
2227 if (tx == ty) switch(tx) {
2228 case t_INT:
2229 if (!is_universal_constant(y)) { affii(x,y); return; }
2230 /* y = gen_0, gnil, gen_1 or gen_2 */
2231 if (y==gen_0) croak("gen_0");
2232 if (y==gen_1) croak("gen_1");
2233 if (y==gen_m1) croak("gen_m1");
2234 if (y==gen_m2) croak("gen_m2");
2235 if (y==gen_2) croak("gen_2");
2236 croak("gnil)");
2237 case t_REAL: affrr(x,y); return;
2238 case t_INTMOD:
2239 if (!dvdii(gel(x,1),gel(y,1))) pari_err_OP("",x,y);
2240 modiiz(gel(x,2),gel(y,1),gel(y,2)); return;
2241 case t_FRAC:
2242 affii(gel(x,1),gel(y,1));
2243 affii(gel(x,2),gel(y,2)); return;
2244 case t_COMPLEX:
2245 gaffect(gel(x,1),gel(y,1));
2246 gaffect(gel(x,2),gel(y,2)); return;
2247 case t_PADIC:
2248 if (!equalii(gel(x,2),gel(y,2))) pari_err_OP("",x,y);
2249 modiiz(gel(x,4),gel(y,3),gel(y,4));
2250 setvalp(y,valp(x)); return;
2251 case t_QUAD:
2252 if (! ZX_equal(gel(x,1),gel(y,1))) pari_err_OP("",x,y);
2253 affii(gel(x,2),gel(y,2));
2254 affii(gel(x,3),gel(y,3)); return;
2255 case t_VEC: case t_COL: case t_MAT:
2256 lx = lg(x); if (lx != lg(y)) pari_err_DIM("gaffect");
2257 for (i=1; i<lx; i++) gaffect(gel(x,i),gel(y,i));
2258 return;
2259 }
2260
2261 /* Various conversions. Avoid them, use specialized routines ! */
2262
2263 if (!is_const_t(ty)) pari_err_TYPE2("=",x,y);
2264 switch(tx)
2265 {
2266 case t_INT:
2267 switch(ty)
2268 {
2269 case t_REAL:
2270 affir(x,y); break;
2271
2272 case t_INTMOD:
2273 modiiz(x,gel(y,1),gel(y,2)); break;
2274
2275 case t_COMPLEX:
2276 gaffect(x,gel(y,1)); gaffsg(0,gel(y,2)); break;
2277
2278 case t_PADIC:
2279 if (!signe(x)) { padicaff0(y); break; }
2280 av = avma;
2281 setvalp(y, Z_pvalrem(x,gel(y,2),&p1));
2282 affii(modii(p1,gel(y,3)), gel(y,4));
2283 set_avma(av); break;
2284
2285 case t_QUAD: gaffect(x,gel(y,2)); gaffsg(0,gel(y,3)); break;
2286 default: pari_err_TYPE2("=",x,y);
2287 }
2288 break;
2289
2290 case t_REAL:
2291 switch(ty)
2292 {
2293 case t_COMPLEX: gaffect(x,gel(y,1)); gaffsg(0,gel(y,2)); break;
2294 default: pari_err_TYPE2("=",x,y);
2295 }
2296 break;
2297
2298 case t_FRAC:
2299 switch(ty)
2300 {
2301 case t_REAL: rdiviiz(gel(x,1),gel(x,2), y); break;
2302 case t_INTMOD: av = avma;
2303 p1 = Fp_inv(gel(x,2),gel(y,1));
2304 affii(modii(mulii(gel(x,1),p1),gel(y,1)), gel(y,2));
2305 set_avma(av); break;
2306 case t_COMPLEX: gaffect(x,gel(y,1)); gaffsg(0,gel(y,2)); break;
2307 case t_PADIC:
2308 if (!signe(gel(x,1))) { padicaff0(y); break; }
2309 num = gel(x,1);
2310 den = gel(x,2);
2311 av = avma; vx = Z_pvalrem(num, gel(y,2), &num);
2312 if (!vx) vx = -Z_pvalrem(den,gel(y,2),&den);
2313 setvalp(y,vx);
2314 p1 = mulii(num,Fp_inv(den,gel(y,3)));
2315 affii(modii(p1,gel(y,3)), gel(y,4)); set_avma(av); break;
2316 case t_QUAD: gaffect(x,gel(y,2)); gaffsg(0,gel(y,3)); break;
2317 default: pari_err_TYPE2("=",x,y);
2318 }
2319 break;
2320
2321 case t_COMPLEX:
2322 if (!gequal0(gel(x,2))) pari_err_TYPE2("=",x,y);
2323 gaffect(gel(x,1), y);
2324 break;
2325
2326 case t_PADIC:
2327 switch(ty)
2328 {
2329 case t_INTMOD:
2330 av = avma; affii(padic_to_Fp(x, gel(y,1)), gel(y,2));
2331 set_avma(av); break;
2332 default: pari_err_TYPE2("=",x,y);
2333 }
2334 break;
2335
2336 case t_QUAD:
2337 switch(ty)
2338 {
2339 case t_INT: case t_INTMOD: case t_FRAC: case t_PADIC:
2340 pari_err_TYPE2("=",x,y);
2341
2342 case t_REAL:
2343 av = avma; affgr(quadtofp(x,realprec(y)), y); set_avma(av); break;
2344 case t_COMPLEX:
2345 ly = precision(y); if (!ly) pari_err_TYPE2("=",x,y);
2346 av = avma; gaffect(quadtofp(x,ly), y); set_avma(av); break;
2347 default: pari_err_TYPE2("=",x,y);
2348 }
2349 default: pari_err_TYPE2("=",x,y);
2350 }
2351 }
2352
2353 /*******************************************************************/
2354 /* */
2355 /* CONVERSION QUAD --> REAL, COMPLEX OR P-ADIC */
2356 /* */
2357 /*******************************************************************/
2358 GEN
quadtofp(GEN x,long prec)2359 quadtofp(GEN x, long prec)
2360 {
2361 GEN b, mb2, c, z, Q, u = gel(x,2), v = gel(x,3);
2362 pari_sp av;
2363 if (prec < LOWDEFAULTPREC) prec = LOWDEFAULTPREC;
2364 if (isintzero(v)) return cxcompotor(u, prec);
2365 av = avma; Q = gel(x,1); c = gel(Q,2); b = gel(Q,3);
2366 z = sqrtr_abs(itor(quad_disc(x), prec));
2367 shiftr_inplace(z, -1);
2368 mb2 = signe(b) ? real2n(-1, prec): NULL;
2369 if (signe(c) > 0) /* c = -D/4 or (1-D)/4 */
2370 z = mkcomplex(mb2? mb2: real_0(prec), z);
2371 else if (mb2)
2372 z = addrr(z, mb2);
2373 /* z = (-b + sqrt(D)) / 2 */
2374 return gerepileupto(av, gadd(u, gmul(v,z)));
2375 }
2376
2377 static GEN
qtop(GEN x,GEN p,long d)2378 qtop(GEN x, GEN p, long d)
2379 {
2380 GEN z, D, P, b, u = gel(x,2), v = gel(x,3);
2381 pari_sp av;
2382 if (gequal0(v)) return cvtop(u, p, d);
2383 P = gel(x,1);
2384 b = gel(P,3);
2385 av = avma; D = quad_disc(x);
2386 if (absequaliu(p,2)) d += 2;
2387 z = Qp_sqrt(cvtop(D,p,d));
2388 if (!z) pari_err_SQRTN("Qp_sqrt",D);
2389 z = gmul2n(gsub(z, b), -1);
2390
2391 z = gadd(u, gmul(v, z));
2392 if (typ(z) != t_PADIC) /* t_INTMOD for t_QUAD of t_INTMODs... */
2393 z = cvtop(z, p, d);
2394 return gerepileupto(av, z);
2395 }
2396 static GEN
ctop(GEN x,GEN p,long d)2397 ctop(GEN x, GEN p, long d)
2398 {
2399 pari_sp av = avma;
2400 GEN z, u = gel(x,1), v = gel(x,2);
2401 if (isrationalzero(v)) return cvtop(u, p, d);
2402 z = Qp_sqrt(cvtop(gen_m1, p, d - gvaluation(v, p))); /* = I */
2403 if (!z) pari_err_SQRTN("Qp_sqrt",gen_m1);
2404
2405 z = gadd(u, gmul(v, z));
2406 if (typ(z) != t_PADIC) /* t_INTMOD for t_COMPLEX of t_INTMODs... */
2407 z = cvtop(z, p, d);
2408 return gerepileupto(av, z);
2409 }
2410
2411 /* cvtop2(stoi(s), y) */
2412 GEN
cvstop2(long s,GEN y)2413 cvstop2(long s, GEN y)
2414 {
2415 GEN z, p = gel(y,2);
2416 long v, d = signe(gel(y,4))? precp(y): 0;
2417 if (!s) return zeropadic(p, d);
2418 v = z_pvalrem(s, p, &s);
2419 if (d <= 0) return zeropadic(p, v);
2420 z = cgetg(5, t_PADIC);
2421 z[1] = evalprecp(d) | evalvalp(v);
2422 gel(z,2) = p;
2423 gel(z,3) = gel(y,3);
2424 gel(z,4) = modsi(s, gel(y,3)); return z;
2425 }
2426
2427 /* cvtop(x, gel(y,2), precp(y)), shallow */
2428 GEN
cvtop2(GEN x,GEN y)2429 cvtop2(GEN x, GEN y)
2430 {
2431 GEN z, p = gel(y,2);
2432 long v, d = signe(gel(y,4))? precp(y): 0;
2433 switch(typ(x))
2434 {
2435 case t_INT:
2436 if (!signe(x)) return zeropadic(p, d);
2437 if (d <= 0) return zeropadic(p, Z_pval(x,p));
2438 v = Z_pvalrem(x, p, &x);
2439 z = cgetg(5, t_PADIC);
2440 z[1] = evalprecp(d) | evalvalp(v);
2441 gel(z,2) = p;
2442 gel(z,3) = gel(y,3);
2443 gel(z,4) = modii(x, gel(y,3)); return z;
2444
2445 case t_INTMOD:
2446 v = Z_pval(gel(x,1),p); if (v > d) v = d;
2447 return cvtop(gel(x,2), p, v);
2448
2449 case t_FRAC:
2450 {
2451 GEN num, den;
2452 if (d <= 0) return zeropadic(p, Q_pval(x,p));
2453 num = gel(x,1); v = Z_pvalrem(num, p, &num);
2454 den = gel(x,2); if (!v) v = -Z_pvalrem(den, p, &den);
2455 z = cgetg(5, t_PADIC);
2456 z[1] = evalprecp(d) | evalvalp(v);
2457 gel(z,2) = p;
2458 gel(z,3) = gel(y,3);
2459 if (!is_pm1(den)) num = mulii(num, Fp_inv(den, gel(y,3)));
2460 gel(z,4) = modii(num, gel(y,3)); return z;
2461 }
2462 case t_COMPLEX: return ctop(x, p, d);
2463 case t_QUAD: return qtop(x, p, d);
2464 }
2465 pari_err_TYPE("cvtop2",x);
2466 return NULL; /* LCOV_EXCL_LINE */
2467 }
2468
2469 /* assume is_const_t(tx) */
2470 GEN
cvtop(GEN x,GEN p,long d)2471 cvtop(GEN x, GEN p, long d)
2472 {
2473 GEN z;
2474 long v;
2475
2476 if (typ(p) != t_INT) pari_err_TYPE("cvtop",p);
2477 switch(typ(x))
2478 {
2479 case t_INT:
2480 if (!signe(x)) return zeropadic(p, d);
2481 if (d <= 0) return zeropadic(p, Z_pval(x,p));
2482 v = Z_pvalrem(x, p, &x);
2483 z = cgetg(5, t_PADIC);
2484 z[1] = evalprecp(d) | evalvalp(v);
2485 gel(z,2) = icopy(p);
2486 gel(z,3) = powiu(p, d);
2487 gel(z,4) = modii(x, gel(z,3)); return z; /* not memory-clean */
2488
2489 case t_INTMOD:
2490 v = Z_pval(gel(x,1),p); if (v > d) v = d;
2491 return cvtop(gel(x,2), p, v);
2492
2493 case t_FRAC:
2494 {
2495 GEN num, den;
2496 if (d <= 0) return zeropadic(p, Q_pval(x,p));
2497 num = gel(x,1); v = Z_pvalrem(num, p, &num);
2498 den = gel(x,2); if (!v) v = -Z_pvalrem(den, p, &den);
2499 z = cgetg(5, t_PADIC);
2500 z[1] = evalprecp(d) | evalvalp(v);
2501 gel(z,2) = icopy(p);
2502 gel(z,3) = powiu(p, d);
2503 if (!is_pm1(den)) num = mulii(num, Fp_inv(den, gel(z,3)));
2504 gel(z,4) = modii(num, gel(z,3)); return z; /* not memory-clean */
2505 }
2506 case t_COMPLEX: return ctop(x, p, d);
2507 case t_PADIC:
2508 p = gel(x,2); /* override */
2509 if (!signe(gel(x,4))) return zeropadic(p, d);
2510 z = cgetg(5,t_PADIC);
2511 z[1] = x[1]; setprecp(z,d);
2512 gel(z,2) = icopy(p);
2513 gel(z,3) = powiu(p, d);
2514 gel(z,4) = modii(gel(x,4), gel(z,3)); return z;
2515
2516 case t_QUAD: return qtop(x, p, d);
2517 }
2518 pari_err_TYPE("cvtop",x);
2519 return NULL; /* LCOV_EXCL_LINE */
2520 }
2521
2522 GEN
gcvtop(GEN x,GEN p,long r)2523 gcvtop(GEN x, GEN p, long r)
2524 {
2525 long i, lx;
2526 GEN y;
2527
2528 switch(typ(x))
2529 {
2530 case t_POL: case t_SER:
2531 y = cgetg_copy(x, &lx); y[1] = x[1];
2532 for (i=2; i<lx; i++) gel(y,i) = gcvtop(gel(x,i),p,r);
2533 return y;
2534 case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
2535 y = cgetg_copy(x, &lx);
2536 for (i=1; i<lx; i++) gel(y,i) = gcvtop(gel(x,i),p,r);
2537 return y;
2538 }
2539 return cvtop(x,p,r);
2540 }
2541
2542 long
gexpo_safe(GEN x)2543 gexpo_safe(GEN x)
2544 {
2545 long tx = typ(x), lx, e, f, i;
2546
2547 switch(tx)
2548 {
2549 case t_INT:
2550 return expi(x);
2551
2552 case t_FRAC:
2553 return expi(gel(x,1)) - expi(gel(x,2));
2554
2555 case t_REAL:
2556 return expo(x);
2557
2558 case t_COMPLEX:
2559 e = gexpo(gel(x,1));
2560 f = gexpo(gel(x,2)); return maxss(e, f);
2561
2562 case t_QUAD: {
2563 GEN p = gel(x,1); /* mod = X^2 + {0,1}* X - {D/4, (1-D)/4})*/
2564 long d = 1 + expi(gel(p,2))/2; /* ~ expo(sqrt(D)) */
2565 e = gexpo(gel(x,2));
2566 f = gexpo(gel(x,3)) + d; return maxss(e, f);
2567 }
2568 case t_POL: case t_SER:
2569 lx = lg(x); f = -(long)HIGHEXPOBIT;
2570 for (i=2; i<lx; i++) { e=gexpo(gel(x,i)); if (e>f) f=e; }
2571 return f;
2572 case t_VEC: case t_COL: case t_MAT:
2573 lx = lg(x); f = -(long)HIGHEXPOBIT;
2574 for (i=1; i<lx; i++) { e=gexpo(gel(x,i)); if (e>f) f=e; }
2575 return f;
2576 }
2577 return -1-(long)HIGHEXPOBIT;
2578 }
2579 long
gexpo(GEN x)2580 gexpo(GEN x)
2581 {
2582 long e = gexpo_safe(x);
2583 if (e < -(long)HIGHEXPOBIT) pari_err_TYPE("gexpo",x);
2584 return e;
2585 }
2586 GEN
gpexponent(GEN x)2587 gpexponent(GEN x)
2588 {
2589 long e = gexpo(x);
2590 return e == -(long)HIGHEXPOBIT? mkmoo(): stoi(e);
2591 }
2592
2593 long
sizedigit(GEN x)2594 sizedigit(GEN x)
2595 {
2596 return gequal0(x)? 0: (long) ((gexpo(x)+1) * LOG10_2) + 1;
2597 }
2598
2599 /* normalize series. avma is not updated */
2600 GEN
normalize(GEN x)2601 normalize(GEN x)
2602 {
2603 long i, lx = lg(x), vx=varn(x), vp=valp(x);
2604 GEN y, z;
2605
2606 if (typ(x) != t_SER) pari_err_TYPE("normalize",x);
2607 if (lx == 2) { setsigne(x,0); return x; }
2608 if (lx == 3) {
2609 z = gel(x,2);
2610 if (!gequal0(z)) { setsigne(x,1); return x; }
2611 if (isrationalzero(z)) return zeroser(vx,vp+1);
2612 if (isexactzero(z)) {
2613 /* dangerous case: already normalized ? */
2614 if (!signe(x)) return x;
2615 setvalp(x,vp+1); /* no: normalize */
2616 }
2617 setsigne(x,0); return x;
2618 }
2619 for (i=2; i<lx; i++)
2620 if (! isrationalzero(gel(x,i))) break;
2621 if (i == lx) return zeroser(vx,lx-2+vp);
2622 z = gel(x,i);
2623 while (i<lx && isexactzero(gel(x,i))) i++;
2624 if (i == lx)
2625 {
2626 i -= 3; y = x + i;
2627 stackdummy((pari_sp)y, (pari_sp)x);
2628 gel(y,2) = z;
2629 y[1] = evalsigne(0) | evalvalp(lx-2+vp) | evalvarn(vx);
2630 y[0] = evaltyp(t_SER) | _evallg(3);
2631 return y;
2632 }
2633
2634 i -= 2; y = x + i; lx -= i;
2635 y[1] = evalsigne(1) | evalvalp(vp+i) | evalvarn(vx);
2636 y[0] = evaltyp(t_SER) | evallg(lx);
2637
2638 stackdummy((pari_sp)y, (pari_sp)x);
2639 for (i = 2; i < lx; i++)
2640 if (!gequal0(gel(y, i))) return y;
2641 setsigne(y, 0); return y;
2642 }
2643
2644 GEN
normalizepol_approx(GEN x,long lx)2645 normalizepol_approx(GEN x, long lx)
2646 {
2647 long i;
2648 for (i = lx-1; i>1; i--)
2649 if (! gequal0(gel(x,i))) break;
2650 stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
2651 setlg(x, i+1); setsigne(x, i!=1); return x;
2652 }
2653
2654 GEN
normalizepol_lg(GEN x,long lx)2655 normalizepol_lg(GEN x, long lx)
2656 {
2657 long i, LX = 0;
2658 GEN KEEP = NULL;
2659
2660 for (i = lx-1; i>1; i--)
2661 {
2662 GEN z = gel(x,i);
2663 if (! gequal0(z) ) {
2664 if (!LX) LX = i+1;
2665 stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + LX));
2666 x[0] = evaltyp(t_POL) | evallg(LX);
2667 setsigne(x,1); return x;
2668 } else if (!isexactzero(z)) {
2669 if (!LX) LX = i+1; /* to be kept as leading coeff */
2670 } else if (!isrationalzero(z))
2671 KEEP = z; /* to be kept iff all other coeffs are exact 0s */
2672 }
2673 if (!LX) {
2674 if (KEEP) { /* e.g. Pol(Mod(0,2)) */
2675 gel(x,2) = KEEP;
2676 LX = 3;
2677 } else
2678 LX = 2; /* Pol(0) */
2679 }
2680 stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + LX));
2681 x[0] = evaltyp(t_POL) | evallg(LX);
2682 setsigne(x,0); return x;
2683 }
2684
2685 /* normalize polynomial x in place */
2686 GEN
normalizepol(GEN x)2687 normalizepol(GEN x)
2688 {
2689 return normalizepol_lg(x, lg(x));
2690 }
2691
2692 int
gsigne(GEN x)2693 gsigne(GEN x)
2694 {
2695 switch(typ(x))
2696 {
2697 case t_INT: case t_REAL: return signe(x);
2698 case t_FRAC: return signe(gel(x,1));
2699 case t_QUAD:
2700 {
2701 pari_sp av = avma;
2702 GEN T = gel(x,1), a = gel(x,2), b = gel(x,3);
2703 long sa, sb;
2704 if (signe(gel(T,2)) > 0) break;
2705 a = gmul2n(a,1);
2706 if (signe(gel(T,3))) a = gadd(a,b);
2707 /* a + b sqrt(D) > 0 ? */
2708 sa = gsigne(a);
2709 sb = gsigne(b); if (sa == sb) return gc_int(av,sa);
2710 if (sa == 0) return gc_int(av,sb);
2711 if (sb == 0) return gc_int(av,sa);
2712 /* different signs, take conjugate expression */
2713 sb = gsigne(gsub(gsqr(a), gmul(quad_disc(x), gsqr(b))));
2714 return gc_int(av, sb*sa);
2715 }
2716 case t_INFINITY: return inf_get_sign(x);
2717 }
2718 pari_err_TYPE("gsigne",x);
2719 return 0; /* LCOV_EXCL_LINE */
2720 }
2721
2722 /*******************************************************************/
2723 /* */
2724 /* LISTS */
2725 /* */
2726 /*******************************************************************/
2727 /* make sure L can hold l elements, at least doubling the previous max number
2728 * of components. */
2729 static void
ensure_nb(GEN L,long l)2730 ensure_nb(GEN L, long l)
2731 {
2732 long nmax = list_nmax(L), i, lw;
2733 GEN v, w;
2734 if (l <= nmax) return;
2735 if (nmax)
2736 {
2737 nmax <<= 1;
2738 if (l > nmax) nmax = l;
2739 w = list_data(L); lw = lg(w);
2740 v = newblock(nmax+1);
2741 v[0] = w[0];
2742 for (i=1; i < lw; i++) gel(v,i) = gel(w, i);
2743 killblock(w);
2744 }
2745 else /* unallocated */
2746 {
2747 nmax = 32;
2748 if (list_data(L))
2749 pari_err(e_MISC, "store list in variable before appending elements");
2750 v = newblock(nmax+1);
2751 v[0] = evaltyp(t_VEC) | _evallg(1);
2752 }
2753 list_data(L) = v;
2754 L[1] = evaltyp(list_typ(L))|evallg(nmax);
2755 }
2756
2757 void
listkill(GEN L)2758 listkill(GEN L)
2759 {
2760
2761 if (typ(L) != t_LIST) pari_err_TYPE("listkill",L);
2762 if (list_nmax(L)) {
2763 GEN v = list_data(L);
2764 long i, l = lg(v);
2765 for (i=1; i<l; i++) gunclone_deep(gel(v,i));
2766 killblock(v);
2767 L[1] = evaltyp(list_typ(L));
2768 list_data(L) = NULL;
2769 }
2770 }
2771
2772 GEN
mklist_typ(long t)2773 mklist_typ(long t)
2774 {
2775 GEN L = cgetg(3,t_LIST);
2776 L[1] = evaltyp(t);
2777 list_data(L) = NULL; return L;
2778 }
2779
2780 GEN
mklist(void)2781 mklist(void)
2782 {
2783 return mklist_typ(t_LIST_RAW);
2784 }
2785
2786 GEN
mkmap(void)2787 mkmap(void)
2788 {
2789 return mklist_typ(t_LIST_MAP);
2790 }
2791
2792 /* return a list with single element x, allocated on stack */
2793 GEN
mklistcopy(GEN x)2794 mklistcopy(GEN x)
2795 {
2796 GEN y = mklist();
2797 list_data(y) = mkveccopy(x);
2798 return y;
2799 }
2800
2801 GEN
listcreate_gp(long n)2802 listcreate_gp(long n)
2803 {
2804 (void) n; return mklist();
2805 }
2806
2807 GEN
listput(GEN L,GEN x,long index)2808 listput(GEN L, GEN x, long index)
2809 {
2810 long l;
2811 GEN z;
2812
2813 if (index < 0) pari_err_COMPONENT("listput", "<", gen_0, stoi(index));
2814 z = list_data(L);
2815 l = z? lg(z): 1;
2816
2817 x = gclone(x);
2818 if (!index || index >= l)
2819 {
2820 ensure_nb(L, l);
2821 z = list_data(L); /* it may change ! */
2822 index = l;
2823 l++;
2824 } else
2825 gunclone_deep( gel(z, index) );
2826 gel(z,index) = x;
2827 z[0] = evaltyp(t_VEC) | evallg(l); /*must be after gel(z,index) is set*/
2828 return gel(z,index);
2829 }
2830
2831 GEN
listput0(GEN L,GEN x,long index)2832 listput0(GEN L, GEN x, long index)
2833 {
2834 if (typ(L) != t_LIST || list_typ(L) != t_LIST_RAW)
2835 pari_err_TYPE("listput",L);
2836 return listput(L, x, index);
2837 }
2838
2839 GEN
listinsert(GEN L,GEN x,long index)2840 listinsert(GEN L, GEN x, long index)
2841 {
2842 long l, i;
2843 GEN z;
2844
2845 if (typ(L) != t_LIST || list_typ(L) != t_LIST_RAW)
2846 pari_err_TYPE("listinsert",L);
2847 z = list_data(L); l = z? lg(z): 1;
2848 if (index <= 0) pari_err_COMPONENT("listinsert", "<=", gen_0, stoi(index));
2849 if (index > l) index = l;
2850 ensure_nb(L, l);
2851 BLOCK_SIGINT_START
2852 z = list_data(L);
2853 for (i=l; i > index; i--) gel(z,i) = gel(z,i-1);
2854 z[0] = evaltyp(t_VEC) | evallg(l+1);
2855 gel(z,index) = gclone(x);
2856 BLOCK_SIGINT_END
2857 return gel(z,index);
2858 }
2859
2860 void
listpop(GEN L,long index)2861 listpop(GEN L, long index)
2862 {
2863 long l, i;
2864 GEN z;
2865
2866 if (typ(L) != t_LIST) pari_err_TYPE("listinsert",L);
2867 if (index < 0) pari_err_COMPONENT("listpop", "<", gen_0, stoi(index));
2868 z = list_data(L);
2869 if (!z || (l = lg(z)-1) == 0) return;
2870
2871 if (!index || index > l) index = l;
2872 BLOCK_SIGINT_START
2873 gunclone_deep( gel(z, index) );
2874 z[0] = evaltyp(t_VEC) | evallg(l);
2875 for (i=index; i < l; i++) z[i] = z[i+1];
2876 BLOCK_SIGINT_END
2877 }
2878
2879 void
listpop0(GEN L,long index)2880 listpop0(GEN L, long index)
2881 {
2882 if (typ(L) != t_LIST || list_typ(L) != t_LIST_RAW)
2883 pari_err_TYPE("listpop",L);
2884 listpop(L, index);
2885 }
2886
2887 /* return a copy fully allocated on stack. gclone from changevalue is
2888 * supposed to malloc() it */
2889 GEN
gtolist(GEN x)2890 gtolist(GEN x)
2891 {
2892 GEN y;
2893
2894 if (!x) return mklist();
2895 switch(typ(x))
2896 {
2897 case t_VEC: case t_COL:
2898 y = mklist();
2899 if (lg(x) == 1) return y;
2900 list_data(y) = gcopy(x);
2901 settyp(list_data(y), t_VEC);
2902 return y;
2903 case t_LIST:
2904 y = mklist();
2905 list_data(y) = list_data(x)? gcopy(list_data(x)): NULL;
2906 return y;
2907 default:
2908 return mklistcopy(x);
2909 }
2910 }
2911
2912 void
listsort(GEN L,long flag)2913 listsort(GEN L, long flag)
2914 {
2915 long i, l;
2916 pari_sp av = avma;
2917 GEN perm, v, vnew;
2918
2919 if (typ(L) != t_LIST) pari_err_TYPE("listsort",L);
2920 v = list_data(L); l = v? lg(v): 1;
2921 if (l < 3) return;
2922 if (flag)
2923 {
2924 long lnew;
2925 perm = gen_indexsort_uniq(L, (void*)&cmp_universal, cmp_nodata);
2926 lnew = lg(perm); /* may have changed since 'uniq' */
2927 vnew = cgetg(lnew,t_VEC);
2928 for (i=1; i<lnew; i++) {
2929 long c = perm[i];
2930 gel(vnew,i) = gel(v,c);
2931 gel(v,c) = NULL;
2932 }
2933 if (l != lnew) { /* was shortened */
2934 for (i=1; i<l; i++)
2935 if (gel(v,i)) gunclone_deep(gel(v,i));
2936 l = lnew;
2937 }
2938 }
2939 else
2940 {
2941 perm = gen_indexsort(L, (void*)&cmp_universal, cmp_nodata);
2942 vnew = cgetg(l,t_VEC);
2943 for (i=1; i<l; i++) gel(vnew,i) = gel(v,perm[i]);
2944 }
2945 for (i=1; i<l; i++) gel(v,i) = gel(vnew,i);
2946 v[0] = vnew[0]; set_avma(av);
2947 }
2948