Line data Source code
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 : /** (third part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : /********************************************************************/
25 : /** **/
26 : /** PRINCIPAL VARIABLE NUMBER **/
27 : /** **/
28 : /********************************************************************/
29 : static void
30 75691 : recvar(hashtable *h, GEN x)
31 : {
32 75691 : long i = 1, lx = lg(x);
33 : void *v;
34 75691 : switch(typ(x))
35 : {
36 12649 : case t_POL: case t_SER:
37 12649 : v = (void*)varn(x);
38 12649 : if (!hash_search(h, v)) hash_insert(h, v, NULL);
39 12649 : i = 2; break;
40 1309 : case t_POLMOD: case t_RFRAC:
41 1309 : case t_VEC: case t_COL: case t_MAT: break;
42 14 : case t_LIST:
43 14 : x = list_data(x);
44 14 : lx = x? lg(x): 1; break;
45 61719 : default:
46 61719 : return;
47 : }
48 88319 : for (; i < lx; i++) recvar(h, gel(x,i));
49 : }
50 :
51 : GEN
52 1344 : variables_vecsmall(GEN x)
53 : {
54 1344 : hashtable *h = hash_create_ulong(100, 1);
55 1344 : recvar(h, x);
56 1344 : return vars_sort_inplace(hash_keys(h));
57 : }
58 :
59 : GEN
60 98 : variables_vec(GEN x)
61 98 : { return x? vars_to_RgXV(variables_vecsmall(x)): gpolvar(NULL); }
62 :
63 : long
64 133413089 : gvar(GEN x)
65 : {
66 : long i, v, w, lx;
67 133413089 : switch(typ(x))
68 : {
69 50944266 : case t_POL: case t_SER: return varn(x);
70 109721 : case t_POLMOD: return varn(gel(x,1));
71 14706583 : case t_RFRAC: return varn(gel(x,2));
72 3953023 : case t_VEC: case t_COL: case t_MAT:
73 3953023 : lx = lg(x); break;
74 14 : case t_LIST:
75 14 : x = list_data(x);
76 14 : lx = x? lg(x): 1; break;
77 63699482 : default:
78 63699482 : return NO_VARIABLE;
79 : }
80 3953037 : v = NO_VARIABLE;
81 33800738 : for (i=1; i < lx; i++) { w = gvar(gel(x,i)); if (varncmp(w,v) < 0) v = w; }
82 3953037 : return v;
83 : }
84 : /* T main polynomial in R[X], A auxiliary in R[X] (possibly degree 0).
85 : * Guess and return the main variable of R */
86 : static long
87 10465 : var2_aux(GEN T, GEN A)
88 : {
89 10465 : long a = gvar2(T);
90 10465 : long b = (typ(A) == t_POL && varn(A) == varn(T))? gvar2(A): gvar(A);
91 10465 : if (varncmp(a, b) > 0) a = b;
92 10465 : return a;
93 : }
94 : static long
95 6895 : var2_rfrac(GEN x) { return var2_aux(gel(x,2), gel(x,1)); }
96 : static long
97 3570 : var2_polmod(GEN x) { return var2_aux(gel(x,1), gel(x,2)); }
98 :
99 : /* main variable of x, with the convention that the "natural" main
100 : * variable of a POLMOD is mute, so we want the next one. */
101 : static long
102 62846 : gvar9(GEN x)
103 62846 : { return (typ(x) == t_POLMOD)? var2_polmod(x): gvar(x); }
104 :
105 : /* main variable of the ring over wich x is defined */
106 : long
107 20672254 : gvar2(GEN x)
108 : {
109 : long i, v, w;
110 20672254 : switch(typ(x))
111 : {
112 7 : case t_POLMOD:
113 7 : return var2_polmod(x);
114 19229 : case t_POL: case t_SER:
115 19229 : v = NO_VARIABLE;
116 81011 : for (i=2; i < lg(x); i++) {
117 61782 : w = gvar9(gel(x,i));
118 61782 : if (varncmp(w,v) < 0) v=w;
119 : }
120 19229 : return v;
121 6895 : case t_RFRAC:
122 6895 : return var2_rfrac(x);
123 49 : case t_VEC: case t_COL: case t_MAT:
124 49 : v = NO_VARIABLE;
125 147 : for (i=1; i < lg(x); i++) {
126 98 : w = gvar2(gel(x,i));
127 98 : if (varncmp(w,v)<0) v=w;
128 : }
129 49 : return v;
130 : }
131 20646074 : return NO_VARIABLE;
132 : }
133 :
134 : /*******************************************************************/
135 : /* */
136 : /* PRECISION OF SCALAR OBJECTS */
137 : /* */
138 : /*******************************************************************/
139 : static long
140 9762998 : prec0(long e) { return (e < 0)? nbits2prec(-e): LOWDEFAULTPREC; }
141 : static long
142 985697394 : precREAL(GEN x) { return signe(x) ? realprec(x): prec0(expo(x)); }
143 : /* x t_REAL, y an exact noncomplex type. Return precision(|x| + |y|) */
144 : static long
145 2007084 : precrealexact(GEN x, GEN y)
146 : {
147 2007084 : long lx, ey = gexpo(y), ex, e;
148 2007084 : if (ey == -(long)HIGHEXPOBIT) return precREAL(x);
149 189297 : ex = expo(x);
150 189297 : e = ey - ex;
151 189297 : if (!signe(x)) return prec0((e >= 0)? -e: ex);
152 189248 : lx = realprec(x);
153 189248 : return (e > 0)? lx + nbits2extraprec(e): lx;
154 : }
155 : static long
156 23097458 : precCOMPLEX(GEN z)
157 : { /* ~ precision(|x| + |y|) */
158 23097458 : GEN x = gel(z,1), y = gel(z,2);
159 : long e, ex, ey, lz, lx, ly;
160 23097458 : if (typ(x) != t_REAL) {
161 2794205 : if (typ(y) != t_REAL) return 0;
162 1993795 : return precrealexact(y, x);
163 : }
164 20303253 : if (typ(y) != t_REAL) return precrealexact(x, y);
165 : /* x, y are t_REALs, cf addrr_sign */
166 20289964 : ex = expo(x);
167 20289964 : ey = expo(y);
168 20289964 : e = ey - ex;
169 20289964 : if (!signe(x)) {
170 596234 : if (!signe(y)) return prec0( minss(ex,ey) );
171 596136 : if (e <= 0) return prec0(ex);
172 596089 : lz = nbits2prec(e);
173 596089 : ly = realprec(y); if (lz > ly) lz = ly;
174 596089 : return lz;
175 : }
176 19693730 : if (!signe(y)) {
177 77236 : if (e >= 0) return prec0(ey);
178 77229 : lz = nbits2prec(-e);
179 77229 : lx = realprec(x); if (lz > lx) lz = lx;
180 77229 : return lz;
181 : }
182 19616494 : if (e < 0) { swap(x, y); e = -e; }
183 19616494 : lx = realprec(x);
184 19616494 : ly = realprec(y);
185 19616494 : if (e) {
186 16346046 : long d = nbits2extraprec(e), l = ly-d;
187 16346046 : return (l > lx)? lx + d: ly;
188 : }
189 3270448 : return minss(lx, ly);
190 : }
191 : long
192 997535321 : precision(GEN z)
193 : {
194 997535321 : switch(typ(z))
195 : {
196 979765841 : case t_REAL: return precREAL(z);
197 17347540 : case t_COMPLEX: return precCOMPLEX(z);
198 : }
199 421940 : return 0;
200 : }
201 :
202 : long
203 15067335 : gprecision(GEN x)
204 : {
205 : long i, k, l;
206 :
207 15067335 : switch(typ(x))
208 : {
209 4113766 : case t_REAL: return precREAL(x);
210 5749918 : case t_COMPLEX: return precCOMPLEX(x);
211 963123 : case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT:
212 : case t_PADIC: case t_QUAD: case t_POLMOD:
213 963123 : return 0;
214 :
215 602 : case t_POL: case t_SER:
216 602 : k = LONG_MAX;
217 2121 : for (i=lg(x)-1; i>1; i--)
218 : {
219 1519 : l = gprecision(gel(x,i));
220 1519 : if (l && l<k) k = l;
221 : }
222 602 : return (k==LONG_MAX)? 0: k;
223 4239856 : case t_VEC: case t_COL: case t_MAT:
224 4239856 : k = LONG_MAX;
225 15494215 : for (i=lg(x)-1; i>0; i--)
226 : {
227 11254359 : l = gprecision(gel(x,i));
228 11254359 : if (l && l<k) k = l;
229 : }
230 4239856 : return (k==LONG_MAX)? 0: k;
231 :
232 7 : case t_RFRAC:
233 : {
234 7 : k=gprecision(gel(x,1));
235 7 : l=gprecision(gel(x,2)); if (l && (!k || l<k)) k=l;
236 7 : return k;
237 : }
238 7 : case t_QFB:
239 7 : return gprecision(gel(x,4));
240 : }
241 56 : return 0;
242 : }
243 :
244 : static long
245 371 : vec_padicprec_relative(GEN x, long imin)
246 : {
247 : long s, t, i;
248 1197 : for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
249 : {
250 826 : t = padicprec_relative(gel(x,i)); if (t<s) s = t;
251 : }
252 371 : return s;
253 : }
254 : /* RELATIVE padic precision. Only accept decent types: don't try to make sense
255 : * of everything like padicprec */
256 : long
257 2275 : padicprec_relative(GEN x)
258 : {
259 2275 : switch(typ(x))
260 : {
261 413 : case t_INT: case t_FRAC:
262 413 : return LONG_MAX;
263 1491 : case t_PADIC:
264 1491 : return signe(padic_u(x))? precp(x): 0;
265 224 : case t_POLMOD: case t_VEC: case t_COL: case t_MAT:
266 224 : return vec_padicprec_relative(x, 1);
267 147 : case t_POL: case t_SER:
268 147 : return vec_padicprec_relative(x, 2);
269 : }
270 0 : pari_err_TYPE("padicprec_relative",x);
271 : return 0; /* LCOV_EXCL_LINE */
272 : }
273 :
274 : static long
275 826 : vec_padicprec(GEN x, GEN p, long imin)
276 : {
277 : long s, t, i;
278 4760 : for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
279 : {
280 3934 : t = padicprec(gel(x,i),p); if (t<s) s = t;
281 : }
282 826 : return s;
283 : }
284 : static long
285 14 : vec_serprec(GEN x, long v, long imin)
286 : {
287 : long s, t, i;
288 42 : for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
289 : {
290 28 : t = serprec(gel(x,i),v); if (t<s) s = t;
291 : }
292 14 : return s;
293 : }
294 :
295 : /* ABSOLUTE padic precision */
296 : long
297 4172 : padicprec(GEN x, GEN p)
298 : {
299 4172 : if (typ(p) != t_INT) pari_err_TYPE("padicprec",p);
300 4165 : switch(typ(x))
301 : {
302 42 : case t_INT: case t_FRAC:
303 42 : return LONG_MAX;
304 :
305 7 : case t_INTMOD:
306 7 : return Z_pval(gel(x,1),p);
307 :
308 3290 : case t_PADIC:
309 3290 : if (!equalii(padic_p(x),p)) pari_err_MODULUS("padicprec", padic_p(x), p);
310 3283 : return precp(x)+valp(x);
311 :
312 14 : case t_POL: case t_SER:
313 14 : return vec_padicprec(x, p, 2);
314 812 : case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_RFRAC:
315 : case t_VEC: case t_COL: case t_MAT:
316 812 : return vec_padicprec(x, p, 1);
317 : }
318 0 : pari_err_TYPE("padicprec",x);
319 : return 0; /* LCOV_EXCL_LINE */
320 : }
321 : GEN
322 105 : gppadicprec(GEN x, GEN p)
323 : {
324 105 : long v = padicprec(x,p);
325 91 : return v == LONG_MAX? mkoo(): stoi(v);
326 : }
327 :
328 : /* ABSOLUTE X-adic precision */
329 : long
330 70 : serprec(GEN x, long v)
331 : {
332 : long w;
333 70 : switch(typ(x))
334 : {
335 21 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
336 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_QFB:
337 21 : return LONG_MAX;
338 :
339 7 : case t_POL:
340 7 : w = varn(x);
341 7 : if (varncmp(v,w) <= 0) return LONG_MAX;
342 7 : return vec_serprec(x, v, 2);
343 42 : case t_SER:
344 42 : w = varn(x);
345 42 : if (w == v)
346 : {
347 35 : long l = lg(x); /* Mod(0,2) + O(x) */
348 35 : if (l == 3 && !signe(x) && !isinexact(gel(x,2))) l--;
349 35 : return l - 2 + valser(x);
350 : }
351 7 : if (varncmp(v,w) < 0) return LONG_MAX;
352 7 : return vec_serprec(x, v, 2);
353 0 : case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
354 0 : return vec_serprec(x, v, 1);
355 : }
356 0 : pari_err_TYPE("serprec",x);
357 : return 0; /* LCOV_EXCL_LINE */
358 : }
359 : GEN
360 42 : gpserprec(GEN x, long v)
361 : {
362 42 : long p = serprec(x,v);
363 42 : return p == LONG_MAX? mkoo(): stoi(p);
364 : }
365 :
366 : /* Degree of x (scalar, t_POL, t_RFRAC) wrt variable v if v >= 0,
367 : * wrt to main variable if v < 0. */
368 : long
369 180439 : poldegree(GEN x, long v)
370 : {
371 180439 : const long DEGREE0 = -LONG_MAX;
372 180439 : long tx = typ(x), lx,w,i,d;
373 :
374 180439 : if (is_scalar_t(tx)) return gequal0(x)? DEGREE0: 0;
375 174149 : switch(tx)
376 : {
377 174030 : case t_POL:
378 174030 : if (!signe(x)) return DEGREE0;
379 159687 : w = varn(x);
380 159687 : if (v < 0 || v == w) return degpol(x);
381 12891 : if (varncmp(v, w) < 0) return 0;
382 12835 : lx = lg(x); d = DEGREE0;
383 77796 : for (i=2; i<lx; i++)
384 : {
385 64961 : long e = poldegree(gel(x,i), v);
386 64961 : if (e > d) d = e;
387 : }
388 12835 : return d;
389 :
390 119 : case t_RFRAC:
391 : {
392 119 : GEN a = gel(x,1), b = gel(x,2);
393 119 : if (gequal0(a)) return DEGREE0;
394 112 : if (v < 0)
395 : {
396 112 : v = varn(b); d = -degpol(b);
397 112 : if (typ(a) == t_POL && varn(a) == v) d += degpol(a);
398 112 : return d;
399 : }
400 0 : return poldegree(a,v) - poldegree(b,v);
401 : }
402 : }
403 0 : pari_err_TYPE("degree",x);
404 : return 0; /* LCOV_EXCL_LINE */
405 : }
406 : GEN
407 31392 : gppoldegree(GEN x, long v)
408 : {
409 31392 : long d = poldegree(x,v);
410 31392 : return d == -LONG_MAX? mkmoo(): stoi(d);
411 : }
412 :
413 : /* assume v >= 0 and x is a POLYNOMIAL in v, return deg_v(x) */
414 : long
415 566892 : RgX_degree(GEN x, long v)
416 : {
417 566892 : long tx = typ(x), lx, w, i, d;
418 :
419 566892 : if (is_scalar_t(tx)) return gequal0(x)? -1: 0;
420 344351 : switch(tx)
421 : {
422 344351 : case t_POL:
423 344351 : if (!signe(x)) return -1;
424 344330 : w = varn(x);
425 344330 : if (v == w) return degpol(x);
426 128136 : if (varncmp(v, w) < 0) return 0;
427 128136 : lx = lg(x); d = -1;
428 548867 : for (i=2; i<lx; i++)
429 : {
430 420731 : long e = RgX_degree(gel(x,i), v);
431 420731 : if (e > d) d = e;
432 : }
433 128136 : return d;
434 :
435 0 : case t_RFRAC:
436 0 : w = varn(gel(x,2));
437 0 : if (varncmp(v, w) < 0) return 0;
438 0 : if (RgX_degree(gel(x,2),v)) pari_err_TYPE("RgX_degree", x);
439 0 : return RgX_degree(gel(x,1),v);
440 : }
441 0 : pari_err_TYPE("RgX_degree",x);
442 : return 0; /* LCOV_EXCL_LINE */
443 : }
444 :
445 : long
446 11314 : degree(GEN x)
447 : {
448 11314 : return poldegree(x,-1);
449 : }
450 :
451 : /* If v<0, leading coeff with respect to the main variable, otherwise wrt v. */
452 : GEN
453 1211 : pollead(GEN x, long v)
454 : {
455 1211 : long tx = typ(x), w;
456 : pari_sp av;
457 :
458 1211 : if (is_scalar_t(tx)) return gcopy(x);
459 1211 : w = varn(x);
460 1211 : switch(tx)
461 : {
462 1176 : case t_POL:
463 1176 : if (v < 0 || v == w)
464 : {
465 1141 : long l = lg(x);
466 1141 : return (l==2)? gen_0: gcopy(gel(x,l-1));
467 : }
468 35 : break;
469 :
470 35 : case t_SER:
471 35 : if (v < 0 || v == w) return signe(x)? gcopy(gel(x,2)): gen_0;
472 14 : if (varncmp(v, w) > 0) x = polcoef_i(x, valser(x), v);
473 14 : break;
474 :
475 0 : default:
476 0 : pari_err_TYPE("pollead",x);
477 : return NULL; /* LCOV_EXCL_LINE */
478 : }
479 49 : if (varncmp(v, w) < 0) return gcopy(x);
480 28 : av = avma; w = fetch_var_higher();
481 28 : x = gsubst(x, v, pol_x(w));
482 28 : x = pollead(x, w);
483 28 : delete_var(); return gc_upto(av, x);
484 : }
485 :
486 : /* returns 1 if there's a real component in the structure, 0 otherwise */
487 : int
488 14833 : isinexactreal(GEN x)
489 : {
490 : long i;
491 14833 : switch(typ(x))
492 : {
493 1246 : case t_REAL: return 1;
494 2723 : case t_COMPLEX: return (typ(gel(x,1))==t_REAL || typ(gel(x,2))==t_REAL);
495 :
496 10276 : case t_INT: case t_INTMOD: case t_FRAC:
497 : case t_FFELT: case t_PADIC: case t_QUAD:
498 10276 : case t_QFB: return 0;
499 :
500 0 : case t_RFRAC: case t_POLMOD:
501 0 : return isinexactreal(gel(x,1)) || isinexactreal(gel(x,2));
502 :
503 588 : case t_POL: case t_SER:
504 5411 : for (i=lg(x)-1; i>1; i--)
505 4872 : if (isinexactreal(gel(x,i))) return 1;
506 539 : return 0;
507 :
508 0 : case t_VEC: case t_COL: case t_MAT:
509 0 : for (i=lg(x)-1; i>0; i--)
510 0 : if (isinexactreal(gel(x,i))) return 1;
511 0 : return 0;
512 0 : default: return 0;
513 : }
514 : }
515 : /* Check if x is approximately real with precision e */
516 : int
517 1891216 : isrealappr(GEN x, long e)
518 : {
519 : long i;
520 1891216 : switch(typ(x))
521 : {
522 699723 : case t_INT: case t_REAL: case t_FRAC:
523 699723 : return 1;
524 1191493 : case t_COMPLEX:
525 1191493 : return (gexpo(gel(x,2)) < e);
526 :
527 0 : case t_POL: case t_SER:
528 0 : for (i=lg(x)-1; i>1; i--)
529 0 : if (! isrealappr(gel(x,i),e)) return 0;
530 0 : return 1;
531 :
532 0 : case t_RFRAC: case t_POLMOD:
533 0 : return isrealappr(gel(x,1),e) && isrealappr(gel(x,2),e);
534 :
535 0 : case t_VEC: case t_COL: case t_MAT:
536 0 : for (i=lg(x)-1; i>0; i--)
537 0 : if (! isrealappr(gel(x,i),e)) return 0;
538 0 : return 1;
539 0 : default: pari_err_TYPE("isrealappr",x); return 0;
540 : }
541 : }
542 :
543 : /* returns 1 if there's an inexact component in the structure, and
544 : * 0 otherwise. */
545 : int
546 129401112 : isinexact(GEN x)
547 : {
548 : long lx, i;
549 :
550 129401112 : switch(typ(x))
551 : {
552 642967 : case t_REAL: case t_PADIC: case t_SER:
553 642967 : return 1;
554 88079515 : case t_INT: case t_INTMOD: case t_FFELT: case t_FRAC:
555 : case t_QFB:
556 88079515 : return 0;
557 2445728 : case t_COMPLEX: case t_QUAD: case t_RFRAC: case t_POLMOD:
558 2445728 : return isinexact(gel(x,1)) || isinexact(gel(x,2));
559 38214061 : case t_POL:
560 123751633 : for (i=lg(x)-1; i>1; i--)
561 85707724 : if (isinexact(gel(x,i))) return 1;
562 38043909 : return 0;
563 18841 : case t_VEC: case t_COL: case t_MAT:
564 22999 : for (i=lg(x)-1; i>0; i--)
565 22488 : if (isinexact(gel(x,i))) return 1;
566 511 : return 0;
567 0 : case t_LIST:
568 0 : x = list_data(x); lx = x? lg(x): 1;
569 0 : for (i=1; i<lx; i++)
570 0 : if (isinexact(gel(x,i))) return 1;
571 0 : return 0;
572 : }
573 0 : return 0;
574 : }
575 :
576 : int
577 0 : isrationalzeroscalar(GEN g)
578 : {
579 0 : switch (typ(g))
580 : {
581 0 : case t_INT: return !signe(g);
582 0 : case t_COMPLEX: return isintzero(gel(g,1)) && isintzero(gel(g,2));
583 0 : case t_QUAD: return isintzero(gel(g,2)) && isintzero(gel(g,3));
584 : }
585 0 : return 0;
586 : }
587 :
588 : int
589 0 : iscomplex(GEN x)
590 : {
591 0 : switch(typ(x))
592 : {
593 0 : case t_INT: case t_REAL: case t_FRAC:
594 0 : return 0;
595 0 : case t_COMPLEX:
596 0 : return !gequal0(gel(x,2));
597 0 : case t_QUAD:
598 0 : return signe(gmael(x,1,2)) > 0;
599 : }
600 0 : pari_err_TYPE("iscomplex",x);
601 : return 0; /* LCOV_EXCL_LINE */
602 : }
603 :
604 : /*******************************************************************/
605 : /* */
606 : /* GENERIC REMAINDER */
607 : /* */
608 : /*******************************************************************/
609 : static int
610 1099 : is_realquad(GEN x) { GEN Q = gel(x,1); return signe(gel(Q,2)) < 0; }
611 : static int
612 188978 : is_realext(GEN x)
613 188978 : { long t = typ(x);
614 188978 : return (t == t_QUAD)? is_realquad(x): is_real_t(t);
615 : }
616 : /* euclidean quotient for scalars of admissible types */
617 : static GEN
618 875 : _quot(GEN x, GEN y)
619 : {
620 875 : GEN q = gdiv(x,y), f = gfloor(q);
621 637 : if (gsigne(y) < 0 && !gequal(f,q)) f = addiu(f, 1);
622 637 : return f;
623 : }
624 : /* y t_REAL, x \ y */
625 : static GEN
626 70 : _quotsr(long x, GEN y)
627 : {
628 : GEN q, f;
629 70 : if (!x) return gen_0;
630 70 : q = divsr(x,y); f = floorr(q);
631 70 : if (signe(y) < 0 && signe(subir(f,q))) f = addiu(f, 1);
632 70 : return f;
633 : }
634 : /* x t_REAL, x \ y */
635 : static GEN
636 280 : _quotrs(GEN x, long y)
637 : {
638 280 : GEN q = divrs(x,y), f = floorr(q);
639 280 : if (y < 0 && signe(subir(f,q))) f = addiu(f, 1);
640 280 : return f;
641 : }
642 : static GEN
643 7 : _quotri(GEN x, GEN y)
644 : {
645 7 : GEN q = divri(x,y), f = floorr(q);
646 7 : if (signe(y) < 0 && signe(subir(f,q))) f = addiu(f, 1);
647 7 : return f;
648 : }
649 : static GEN
650 70 : _quotsq(long x, GEN y)
651 : {
652 70 : GEN f = gfloor(gdivsg(x,y));
653 70 : if (gsigne(y) < 0) f = gaddgs(f, 1);
654 70 : return f;
655 : }
656 : static GEN
657 28 : _quotqs(GEN x, long y)
658 : {
659 28 : GEN f = gfloor(gdivgs(x,y));
660 28 : if (y < 0) f = addiu(f, 1);
661 28 : return f;
662 : }
663 :
664 : /* y t_FRAC, x \ y */
665 : static GEN
666 35 : _quotsf(long x, GEN y)
667 35 : { return truedivii(mulis(gel(y,2),x), gel(y,1)); }
668 : /* x t_FRAC, x \ y */
669 : static GEN
670 301 : _quotfs(GEN x, long y)
671 301 : { return truedivii(gel(x,1),mulis(gel(x,2),y)); }
672 : /* x t_FRAC, y t_INT, x \ y */
673 : static GEN
674 7 : _quotfi(GEN x, GEN y)
675 7 : { return truedivii(gel(x,1),mulii(gel(x,2),y)); }
676 :
677 : static GEN
678 777 : quot(GEN x, GEN y)
679 777 : { pari_sp av = avma; return gc_upto(av, _quot(x, y)); }
680 : static GEN
681 266 : quotrs(GEN x, long y)
682 266 : { pari_sp av = avma; return gc_leaf(av, _quotrs(x,y)); }
683 : static GEN
684 301 : quotfs(GEN x, long s)
685 301 : { pari_sp av = avma; return gc_leaf(av, _quotfs(x,s)); }
686 : static GEN
687 35 : quotsr(long x, GEN y)
688 35 : { pari_sp av = avma; return gc_leaf(av, _quotsr(x, y)); }
689 : static GEN
690 35 : quotsf(long x, GEN y)
691 35 : { pari_sp av = avma; return gc_leaf(av, _quotsf(x, y)); }
692 : static GEN
693 35 : quotsq(long x, GEN y)
694 35 : { pari_sp av = avma; return gc_leaf(av, _quotsq(x, y)); }
695 : static GEN
696 14 : quotqs(GEN x, long y)
697 14 : { pari_sp av = avma; return gc_leaf(av, _quotqs(x, y)); }
698 : static GEN
699 7 : quotfi(GEN x, GEN y)
700 7 : { pari_sp av = avma; return gc_leaf(av, _quotfi(x, y)); }
701 : static GEN
702 7 : quotri(GEN x, GEN y)
703 7 : { pari_sp av = avma; return gc_leaf(av, _quotri(x, y)); }
704 :
705 : static GEN
706 14 : modrs(GEN x, long y)
707 : {
708 14 : pari_sp av = avma;
709 14 : GEN q = _quotrs(x,y);
710 14 : if (!signe(q)) { set_avma(av); return rcopy(x); }
711 7 : return gc_leaf(av, subri(x, mulis(q,y)));
712 : }
713 : static GEN
714 35 : modsr(long x, GEN y)
715 : {
716 35 : pari_sp av = avma;
717 35 : GEN q = _quotsr(x,y);
718 35 : if (!signe(q)) return gc_stoi(av, x);
719 7 : return gc_leaf(av, subsr(x, mulir(q,y)));
720 : }
721 : static GEN
722 35 : modsf(long x, GEN y)
723 : {
724 35 : pari_sp av = avma;
725 35 : return gc_upto(av, Qdivii(modii(mulis(gel(y,2),x), gel(y,1)), gel(y,2)));
726 : }
727 :
728 : /* assume y a t_REAL, x a t_INT, t_FRAC or t_REAL.
729 : * Return x mod y or NULL if accuracy error */
730 : GEN
731 0 : modRr_safe(GEN x, GEN y)
732 : {
733 : GEN q, f;
734 : long e;
735 0 : if (isintzero(x)) return gen_0;
736 0 : q = gdiv(x,y); /* t_REAL */
737 :
738 0 : e = expo(q);
739 0 : if (e >= 0 && nbits2prec(e+1) > realprec(q)) return NULL;
740 0 : f = floorr(q);
741 0 : if (signe(y) < 0 && signe(subri(q,f))) f = addiu(f, 1);
742 0 : return signe(f)? gsub(x, mulir(f,y)): x;
743 : }
744 : GEN
745 5105155 : modRr_i(GEN x, GEN y, GEN iy)
746 : {
747 : GEN q, f;
748 : long e;
749 5105155 : if (isintzero(x)) return gen_0;
750 5105155 : q = gmul(x, iy); /* t_REAL */
751 :
752 5105155 : e = expo(q);
753 5105155 : if (e >= 0 && nbits2prec(e+1) > realprec(q)) return NULL;
754 5105155 : f = floorr(q);
755 5105155 : if (signe(y) < 0 && signe(subri(q,f))) f = addiu(f, 1);
756 5105155 : return signe(f)? gsub(x, mulir(f,y)): x;
757 : }
758 :
759 : GEN
760 46175890 : gmod(GEN x, GEN y)
761 : {
762 : pari_sp av;
763 : long ty, tx;
764 : GEN z;
765 :
766 46175890 : tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gmodsg(itos(x),y);
767 1161695 : ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gmodgs(x,itos(y));
768 1756064 : if (is_matvec_t(tx)) pari_APPLY_same(gmod(gel(x,i), y));
769 818744 : if (tx == t_POL || ty == t_POL) return grem(x,y);
770 511328 : if (!is_scalar_t(tx) || !is_scalar_t(ty)) pari_err_TYPE2("%",x,y);
771 511265 : switch(ty)
772 : {
773 510761 : case t_INT:
774 510761 : switch(tx)
775 : {
776 507661 : case t_INT: return modii(x,y);
777 7 : case t_INTMOD: z=cgetg(3, t_INTMOD);
778 7 : gel(z,1) = gcdii(gel(x,1),y);
779 7 : gel(z,2) = modii(gel(x,2),gel(z,1)); return z;
780 491 : case t_FRAC: return Fp_div(gel(x,1),gel(x,2),y);
781 2567 : case t_PADIC: return padic_to_Fp(x, y);
782 14 : case t_QUAD: if (!is_realquad(x)) break;
783 : case t_REAL:
784 14 : av = avma; /* NB: conflicting semantic with lift(x * Mod(1,y)). */
785 14 : return gc_upto(av, gsub(x, gmul(_quot(x,y),y)));
786 : }
787 21 : break;
788 126 : case t_QUAD:
789 126 : if (!is_realquad(y)) break;
790 : case t_REAL: case t_FRAC:
791 189 : if (!is_realext(x)) break;
792 84 : av = avma;
793 84 : return gc_upto(av, gsub(x, gmul(_quot(x,y),y)));
794 : }
795 441 : pari_err_TYPE2("%",x,y);
796 : return NULL; /* LCOV_EXCL_LINE */
797 : }
798 :
799 : GEN
800 22032998 : gmodgs(GEN x, long y)
801 : {
802 : ulong u;
803 22032998 : long i, tx = typ(x);
804 : GEN z;
805 43827220 : if (is_matvec_t(tx)) pari_APPLY_same(gmodgs(gel(x,i), y));
806 19654700 : if (!y) pari_err_INV("gmodgs",gen_0);
807 19654700 : switch(tx)
808 : {
809 19566496 : case t_INT: return modis(x,y);
810 14 : case t_REAL: return modrs(x,y);
811 :
812 21 : case t_INTMOD: z=cgetg(3, t_INTMOD);
813 21 : u = (ulong)labs(y);
814 21 : i = ugcdiu(gel(x,1), u);
815 21 : gel(z,1) = utoi(i);
816 21 : gel(z,2) = modis(gel(x,2), i); return z;
817 :
818 86767 : case t_FRAC:
819 86767 : u = (ulong)labs(y);
820 86767 : return utoi( Fl_div(umodiu(gel(x,1), u),
821 86767 : umodiu(gel(x,2), u), u) );
822 28 : case t_QUAD:
823 : {
824 28 : pari_sp av = avma;
825 28 : if (!is_realquad(x)) break;
826 14 : return gc_upto(av, gsub(x, gmulgs(_quotqs(x,y),y)));
827 : }
828 1318 : case t_PADIC: return padic_to_Fp(x, stoi(y));
829 14 : case t_POL: return scalarpol(Rg_get_0(x), varn(x));
830 14 : case t_POLMOD: return gmul(gen_0,x);
831 : }
832 42 : pari_err_TYPE2("%",x,stoi(y));
833 : return NULL; /* LCOV_EXCL_LINE */
834 : }
835 : GEN
836 45014195 : gmodsg(long x, GEN y)
837 : {
838 45014195 : switch(typ(y))
839 : {
840 45013838 : case t_INT: return modsi(x,y);
841 35 : case t_REAL: return modsr(x,y);
842 35 : case t_FRAC: return modsf(x,y);
843 63 : case t_QUAD:
844 : {
845 63 : pari_sp av = avma;
846 63 : if (!is_realquad(y)) break;
847 35 : return gc_upto(av, gsubsg(x, gmul(_quotsq(x,y),y)));
848 : }
849 112 : case t_POL:
850 112 : if (!signe(y)) pari_err_INV("gmodsg",y);
851 112 : return degpol(y)? gmulsg(x, Rg_get_1(y)): Rg_get_0(y);
852 : }
853 140 : pari_err_TYPE2("%",stoi(x),y);
854 : return NULL; /* LCOV_EXCL_LINE */
855 : }
856 : /* divisibility: return 1 if y | x, 0 otherwise */
857 : int
858 16055 : gdvd(GEN x, GEN y)
859 : {
860 16055 : pari_sp av = avma;
861 16055 : return gc_bool(av, gequal0( gmod(x,y) ));
862 : }
863 :
864 : GEN
865 1074361 : gmodulss(long x, long y)
866 : {
867 1074361 : if (!y) pari_err_INV("%",gen_0);
868 1074354 : y = labs(y);
869 1074354 : retmkintmod(utoi(umodsu(x, y)), utoipos(y));
870 : }
871 : GEN
872 1343357 : gmodulsg(long x, GEN y)
873 : {
874 1343357 : switch(typ(y))
875 : {
876 1127503 : case t_INT:
877 1127503 : if (!is_bigint(y)) return gmodulss(x,itos(y));
878 62380 : retmkintmod(modsi(x,y), absi(y));
879 215847 : case t_POL:
880 215847 : if (!signe(y)) pari_err_INV("%", y);
881 215840 : retmkpolmod(degpol(y)? stoi(x): gen_0,RgX_copy(y));
882 : }
883 7 : pari_err_TYPE2("%",stoi(x),y);
884 : return NULL; /* LCOV_EXCL_LINE */
885 : }
886 : GEN
887 1959902 : gmodulo(GEN x,GEN y)
888 : {
889 1959902 : long tx = typ(x), vx, vy;
890 1959902 : if (tx == t_INT && !is_bigint(x)) return gmodulsg(itos(x), y);
891 1461375 : if (is_matvec_t(tx)) pari_APPLY_same(gmodulo(gel(x,i), y));
892 579256 : switch(typ(y))
893 : {
894 98223 : case t_INT:
895 98223 : if (!is_const_t(tx)) return gmul(x, gmodulsg(1,y));
896 98174 : if (tx == t_INTMOD) return gmod(x,y);
897 98167 : retmkintmod(Rg_to_Fp(x,y), absi(y));
898 481033 : case t_POL:
899 481033 : vx = gvar(x); vy = varn(y);
900 481033 : if (varncmp(vy, vx) > 0) return gmul(x, gmodulsg(1,y));
901 480942 : if (vx == vy && tx == t_POLMOD) return grem(x,y);
902 468062 : retmkpolmod(grem(x,y), RgX_copy(y));
903 : }
904 0 : pari_err_TYPE2("%",x,y);
905 : return NULL; /* LCOV_EXCL_LINE */
906 : }
907 :
908 : /*******************************************************************/
909 : /* */
910 : /* GENERIC EUCLIDEAN DIVISION */
911 : /* */
912 : /*******************************************************************/
913 : GEN
914 6204961 : gdivent(GEN x, GEN y)
915 : {
916 : long tx, ty;
917 6204961 : tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gdiventsg(itos(x),y);
918 2188 : ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gdiventgs(x,itos(y));
919 1960 : if (is_matvec_t(tx)) pari_APPLY_same(gdivent(gel(x,i), y));
920 1610 : if (tx == t_POL || ty == t_POL) return gdeuc(x,y);
921 1148 : switch(ty)
922 : {
923 112 : case t_INT:
924 112 : switch(tx)
925 : {
926 7 : case t_INT: return truedivii(x,y);
927 7 : case t_REAL: return quotri(x,y);
928 7 : case t_FRAC: return quotfi(x,y);
929 21 : case t_QUAD:
930 21 : if (!is_realquad(x)) break;
931 7 : return quot(x,y);
932 : }
933 84 : break;
934 252 : case t_QUAD:
935 252 : if (!is_realext(x) || !is_realquad(y)) break;
936 : case t_REAL: case t_FRAC:
937 252 : return quot(x,y);
938 : }
939 868 : pari_err_TYPE2("\\",x,y);
940 : return NULL; /* LCOV_EXCL_LINE */
941 : }
942 :
943 : GEN
944 328076 : gdiventgs(GEN x, long y)
945 : {
946 328076 : switch(typ(x))
947 : {
948 268440 : case t_INT: return truedivis(x,y);
949 266 : case t_REAL: return quotrs(x,y);
950 301 : case t_FRAC: return quotfs(x,y);
951 42 : case t_QUAD: if (!is_realquad(x)) break;
952 14 : return quotqs(x,y);
953 28 : case t_POL: return gdivgs(x,y);
954 316613 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gdiventgs(gel(x,i),y));
955 : }
956 168 : pari_err_TYPE2("\\",x,stoi(y));
957 : return NULL; /* LCOV_EXCL_LINE */
958 : }
959 : GEN
960 6202773 : gdiventsg(long x, GEN y)
961 : {
962 6202773 : switch(typ(y))
963 : {
964 6202318 : case t_INT: return truedivsi(x,y);
965 35 : case t_REAL: return quotsr(x,y);
966 35 : case t_FRAC: return quotsf(x,y);
967 91 : case t_QUAD: if (!is_realquad(y)) break;
968 35 : return quotsq(x,y);
969 70 : case t_POL:
970 70 : if (!signe(y)) pari_err_INV("gdiventsg",y);
971 70 : return degpol(y)? Rg_get_0(y): gdivsg(x,gel(y,2));
972 : }
973 280 : pari_err_TYPE2("\\",stoi(x),y);
974 : return NULL; /* LCOV_EXCL_LINE */
975 : }
976 :
977 : /* with remainder */
978 : static GEN
979 518 : quotrem(GEN x, GEN y, GEN *r)
980 : {
981 518 : GEN q = quot(x,y);
982 448 : pari_sp av = avma;
983 448 : *r = gc_upto(av, gsub(x, gmul(q,y)));
984 448 : return q;
985 : }
986 :
987 : GEN
988 1064 : gdiventres(GEN x, GEN y)
989 : {
990 1064 : long tx = typ(x), ty = typ(y);
991 : GEN z;
992 :
993 1078 : if (is_matvec_t(tx)) pari_APPLY_same(gdiventres(gel(x,i), y));
994 1057 : z = cgetg(3,t_COL);
995 1057 : if (tx == t_POL || ty == t_POL)
996 : {
997 182 : gel(z,1) = poldivrem(x,y,(GEN*)(z+2));
998 168 : return z;
999 : }
1000 875 : switch(ty)
1001 : {
1002 252 : case t_INT:
1003 252 : switch(tx)
1004 : { /* equal to, but more efficient than next case */
1005 84 : case t_INT:
1006 84 : gel(z,1) = truedvmdii(x,y,(GEN*)(z+2));
1007 84 : return z;
1008 42 : case t_QUAD:
1009 42 : if (!is_realquad(x)) break;
1010 : case t_REAL: case t_FRAC:
1011 63 : gel(z,1) = quotrem(x,y,&gel(z,2));
1012 63 : return z;
1013 : }
1014 105 : break;
1015 154 : case t_QUAD:
1016 154 : if (!is_realext(x) || !is_realquad(y)) break;
1017 : case t_REAL: case t_FRAC:
1018 196 : gel(z,1) = quotrem(x,y,&gel(z,2));
1019 126 : return z;
1020 : }
1021 532 : pari_err_TYPE2("\\",x,y);
1022 : return NULL; /* LCOV_EXCL_LINE */
1023 : }
1024 :
1025 : GEN
1026 1057 : divrem(GEN x, GEN y, long v)
1027 : {
1028 1057 : pari_sp av = avma;
1029 : long vx, vy, vz;
1030 : GEN q, r;
1031 1057 : if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
1032 7 : vz = v; vx = varn(x); vy = varn(y);
1033 7 : if (varncmp(vx, vz) < 0 || varncmp(vy, vz) < 0) vz = fetch_var_higher();
1034 7 : if (vx != vz) x = swap_vars(x, v, vz);
1035 7 : if (vy != vz) y = swap_vars(y, v, vz);
1036 7 : q = RgX_divrem(x,y, &r);
1037 7 : if (vx != v || vy != v)
1038 : {
1039 7 : GEN X = pol_x(v);
1040 7 : q = gsubst(q, vz, X); /* poleval broken for t_RFRAC, subst is safe */
1041 7 : r = gsubst(r, vz, X);
1042 7 : if (vz != v) (void)delete_var();
1043 : }
1044 7 : return gc_GEN(av, mkcol2(q, r));
1045 : }
1046 :
1047 : GEN
1048 68081799 : diviiround(GEN x, GEN y)
1049 : {
1050 68081799 : pari_sp av1, av = avma;
1051 : GEN q,r;
1052 : int fl;
1053 :
1054 68081799 : q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
1055 68081792 : if (r==gen_0) return q;
1056 37196009 : av1 = avma;
1057 37196009 : fl = abscmpii(shifti(r,1),y);
1058 37196009 : set_avma(av1); cgiv(r);
1059 37196009 : if (fl >= 0) /* If 2*|r| >= |y| */
1060 : {
1061 20301996 : long sz = signe(x)*signe(y);
1062 20301996 : if (fl || sz > 0) q = gc_INT(av, addis(q,sz));
1063 : }
1064 37196009 : return q;
1065 : }
1066 :
1067 : static GEN
1068 518 : _abs(GEN x)
1069 : {
1070 518 : if (typ(x) == t_QUAD) return (gsigne(x) < 0)? gneg(x): x;
1071 364 : return R_abs_shallow(x);
1072 : }
1073 :
1074 : /* If x and y are not both scalars, same as gdivent.
1075 : * Otherwise, compute the quotient x/y, rounded to the nearest integer
1076 : * (towards +oo in case of tie). */
1077 : GEN
1078 1481201 : gdivround(GEN x, GEN y)
1079 : {
1080 : pari_sp av;
1081 1481201 : long tx = typ(x), ty = typ(y);
1082 : GEN q, r;
1083 :
1084 1481201 : if (tx == t_INT && ty == t_INT) return diviiround(x,y);
1085 187781 : av = avma;
1086 187781 : if (is_realext(x) && is_realext(y))
1087 : { /* same as diviiround, less efficient */
1088 : pari_sp av1;
1089 : int fl;
1090 259 : q = quotrem(x,y,&r); av1 = avma;
1091 259 : fl = gcmp(gmul2n(_abs(r),1), _abs(y));
1092 259 : set_avma(av1); cgiv(r);
1093 259 : if (fl >= 0) /* If 2*|r| >= |y| */
1094 : {
1095 84 : long sz = gsigne(y);
1096 84 : if (fl || sz > 0) q = gc_upto(av, gaddgs(q, sz));
1097 : }
1098 259 : return q;
1099 : }
1100 1622364 : if (is_matvec_t(tx)) pari_APPLY_same(gdivround(gel(x,i),y));
1101 931 : return gdivent(x,y);
1102 : }
1103 :
1104 : GEN
1105 0 : gdivmod(GEN x, GEN y, GEN *pr)
1106 : {
1107 0 : switch(typ(x))
1108 : {
1109 0 : case t_INT:
1110 0 : switch(typ(y))
1111 : {
1112 0 : case t_INT: return dvmdii(x,y,pr);
1113 0 : case t_POL: *pr=icopy(x); return gen_0;
1114 : }
1115 0 : break;
1116 0 : case t_POL: return poldivrem(x,y,pr);
1117 : }
1118 0 : pari_err_TYPE2("gdivmod",x,y);
1119 : return NULL;/*LCOV_EXCL_LINE*/
1120 : }
1121 :
1122 : /*******************************************************************/
1123 : /* */
1124 : /* SHIFT */
1125 : /* */
1126 : /*******************************************************************/
1127 :
1128 : /* Shift tronque si n<0 (multiplication tronquee par 2^n) */
1129 :
1130 : GEN
1131 47409545 : gshift(GEN x, long n)
1132 : {
1133 47409545 : switch(typ(x))
1134 : {
1135 38982355 : case t_INT: return shifti(x,n);
1136 7698080 : case t_REAL:return shiftr(x,n);
1137 2248816 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gshift(gel(x,i),n));
1138 : }
1139 156030 : return gmul2n(x,n);
1140 : }
1141 :
1142 : /*******************************************************************/
1143 : /* */
1144 : /* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
1145 : /* */
1146 : /*******************************************************************/
1147 :
1148 : /* Convert t_SER --> t_POL, ignoring valser. INTERNAL ! */
1149 : GEN
1150 11869399 : ser2pol_i(GEN x, long lx)
1151 : {
1152 11869399 : long i = lx-1;
1153 : GEN y;
1154 16452224 : while (i > 1 && isrationalzero(gel(x,i))) i--;
1155 11869399 : if (!signe(x))
1156 : { /* danger */
1157 77 : if (i == 1) return zeropol(varn(x));
1158 77 : y = cgetg(3,t_POL); y[1] = x[1] & ~VALSERBITS;
1159 77 : gel(y,2) = gel(x,2); return y;
1160 : }
1161 11869322 : y = cgetg(i+1, t_POL); y[1] = x[1] & ~VALSERBITS;
1162 50672304 : for ( ; i > 1; i--) gel(y,i) = gel(x,i);
1163 11869322 : return y;
1164 : }
1165 :
1166 : GEN
1167 899117 : ser2pol_i_normalize(GEN x, long l, long *v)
1168 : {
1169 899117 : long i = 2, j = l-1, k;
1170 : GEN y;
1171 899152 : while (i < l && gequal0(gel(x,i))) i++;
1172 899117 : *v = i - 2; if (i == l) return zeropol(varn(x));
1173 1148603 : while (j > i && gequal0(gel(x,j))) j--;
1174 899103 : l = j - *v + 1;
1175 899103 : y = cgetg(l, t_POL); y[1] = x[1] & ~VALSERBITS;
1176 4586855 : k = l; while (k > 2) gel(y, --k) = gel(x,j--);
1177 899103 : return y;
1178 : }
1179 :
1180 : GEN
1181 50015 : ser_inv(GEN b)
1182 : {
1183 50015 : pari_sp av = avma;
1184 50015 : long e, l = lg(b);
1185 : GEN x, y;
1186 50015 : y = ser2pol_i_normalize(b, l, &e);
1187 50015 : if (e)
1188 : {
1189 0 : pari_warn(warner,"normalizing a series with 0 leading term");
1190 0 : l -= e; if (l <= 2) pari_err_INV("inv_ser", b);
1191 : }
1192 50015 : y = RgXn_inv_i(y, l-2);
1193 50008 : x = RgX_to_ser(y, l); setvalser(x, - valser(b) - e);
1194 50008 : return gc_GEN(av, x);
1195 : }
1196 :
1197 : /* T t_POL in var v, mod out by T components of x which are t_POL/t_RFRAC in v.
1198 : * Recursively. Make sure that resulting polynomials of degree 0 in v are
1199 : * simplified (map K[X]_0 to K) */
1200 : static GEN
1201 210 : mod_r(GEN x, long v, GEN T)
1202 : {
1203 210 : long w, tx = typ(x);
1204 : GEN y;
1205 :
1206 210 : if (is_const_t(tx)) return x;
1207 189 : switch(tx)
1208 : {
1209 7 : case t_POLMOD:
1210 7 : w = varn(gel(x,1));
1211 7 : if (w == v) pari_err_PRIORITY("subst", gel(x,1), "=", v);
1212 7 : if (varncmp(v, w) < 0) return x;
1213 7 : return gmodulo(mod_r(gel(x,2),v,T), mod_r(gel(x,1),v,T));
1214 7 : case t_SER:
1215 7 : w = varn(x);
1216 7 : if (w == v) break; /* fail */
1217 7 : if (varncmp(v, w) < 0 || ser_isexactzero(x)) return x;
1218 21 : pari_APPLY_ser(mod_r(gel(x,i),v,T));
1219 147 : case t_POL:
1220 147 : w = varn(x);
1221 147 : if (w == v)
1222 : {
1223 119 : x = RgX_rem(x, T);
1224 119 : if (!degpol(x)) x = gel(x,2);
1225 119 : return x;
1226 : }
1227 28 : if (varncmp(v, w) < 0) return x;
1228 98 : pari_APPLY_pol(mod_r(gel(x,i),v,T));
1229 14 : case t_RFRAC:
1230 14 : x = gdiv(mod_r(gel(x,1),v,T), mod_r(gel(x,2),v,T));
1231 14 : if (typ(x) == t_POL && varn(x) == v && lg(x) == 3) x = gel(x,2);
1232 14 : return x;
1233 7 : case t_VEC: case t_COL: case t_MAT:
1234 21 : pari_APPLY_same(mod_r(gel(x,i),v,T));
1235 7 : case t_LIST:
1236 7 : y = mklist();
1237 7 : list_data(y) = list_data(x)? mod_r(list_data(x),v,T): NULL;
1238 7 : return y;
1239 : }
1240 0 : pari_err_TYPE("substpol",x);
1241 : return NULL;/*LCOV_EXCL_LINE*/
1242 : }
1243 : GEN
1244 8302 : gsubstpol(GEN x, GEN T, GEN y)
1245 : {
1246 8302 : pari_sp av = avma;
1247 : long v;
1248 : GEN z;
1249 8302 : if (typ(T) == t_POL && RgX_is_monomial(T) && gequal1(leading_coeff(T)))
1250 : { /* T = t^d */
1251 8267 : long d = degpol(T);
1252 8267 : v = varn(T); z = (d==1)? x: gdeflate(x, v, d);
1253 8253 : if (z) return gc_upto(av, gsubst(z, v, y));
1254 : }
1255 63 : v = fetch_var(); T = simplify_shallow(T);
1256 63 : if (typ(T) == t_RFRAC)
1257 21 : z = gsub(gel(T,1), gmul(pol_x(v), gel(T,2)));
1258 : else
1259 42 : z = gsub(T, pol_x(v));
1260 63 : z = mod_r(x, gvar(T), z);
1261 63 : z = gsubst(z, v, y); (void)delete_var();
1262 63 : return gc_upto(av, z);
1263 : }
1264 :
1265 : long
1266 1256813 : RgX_deflate_order(GEN x)
1267 : {
1268 1256813 : ulong d = 0, i, lx = (ulong)lg(x);
1269 2756400 : for (i=3; i<lx; i++)
1270 2377485 : if (!gequal0(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1271 378915 : return d? (long)d: 1;
1272 : }
1273 : long
1274 530252 : ZX_deflate_order(GEN x)
1275 : {
1276 530252 : ulong d = 0, i, lx = (ulong)lg(x);
1277 1654617 : for (i=3; i<lx; i++)
1278 1454935 : if (signe(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1279 199682 : return d? (long)d: 1;
1280 : }
1281 :
1282 : /* deflate (non-leaf) x recursively */
1283 : static GEN
1284 63 : vdeflate(GEN x, long v, long d)
1285 : {
1286 63 : long i, lx, tx = typ(x);
1287 63 : GEN y = cgetg_copy(x, &lx);
1288 98 : for(i = 1; i < lontyp[tx]; i++) y[i] = x[i];
1289 154 : for (; i < lx; i++)
1290 : {
1291 133 : gel(y,i) = gdeflate(gel(x,i),v,d);
1292 133 : if (!y[i]) return NULL;
1293 : }
1294 21 : return y;
1295 : }
1296 :
1297 : /* don't return NULL if substitution fails (fallback won't be able to handle
1298 : * t_SER anyway), fail with a meaningful message */
1299 : static GEN
1300 6818 : serdeflate(GEN x, long v, long d)
1301 : {
1302 6818 : long V, dy, lx, vx = varn(x);
1303 : pari_sp av;
1304 : GEN y;
1305 6818 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1306 6811 : if (varncmp(vx, v) > 0) return gcopy(x);
1307 6811 : av = avma;
1308 6811 : V = valser(x);
1309 6811 : lx = lg(x);
1310 6811 : if (lx == 2) return zeroser(v, V / d);
1311 6811 : y = ser2pol_i(x, lx);
1312 6811 : dy = degpol(y);
1313 6811 : if (V % d != 0 || (dy > 0 && RgX_deflate_order(y) % d != 0))
1314 : {
1315 14 : const char *s = stack_sprintf("valuation(x) %% %ld", d);
1316 14 : pari_err_DOMAIN("gdeflate", s, "!=", gen_0,x);
1317 : }
1318 6797 : if (dy > 0) y = RgX_deflate(y, d);
1319 6797 : y = RgX_to_ser(y, 3 + (lx-3)/d);
1320 6797 : setvalser(y, V/d); return gc_GEN(av, y);
1321 : }
1322 : static GEN
1323 8302 : poldeflate(GEN x, long v, long d)
1324 : {
1325 8302 : long vx = varn(x);
1326 : pari_sp av;
1327 8302 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1328 8274 : if (varncmp(vx, v) > 0 || degpol(x) <= 0) return gcopy(x);
1329 8239 : av = avma;
1330 : /* x nonconstant */
1331 8239 : if (RgX_deflate_order(x) % d != 0) return NULL;
1332 8211 : return gc_GEN(av, RgX_deflate(x,d));
1333 : }
1334 : static GEN
1335 21 : listdeflate(GEN x, long v, long d)
1336 : {
1337 21 : GEN y = NULL, z = mklist();
1338 21 : if (list_data(x))
1339 : {
1340 14 : y = vdeflate(list_data(x),v,d);
1341 14 : if (!y) return NULL;
1342 : }
1343 14 : list_data(z) = y; return z;
1344 : }
1345 : /* return NULL if substitution fails */
1346 : GEN
1347 15183 : gdeflate(GEN x, long v, long d)
1348 : {
1349 15183 : if (d <= 0) pari_err_DOMAIN("gdeflate", "degree", "<=", gen_0,stoi(d));
1350 15183 : switch(typ(x))
1351 : {
1352 28 : case t_INT:
1353 : case t_REAL:
1354 : case t_INTMOD:
1355 : case t_FRAC:
1356 : case t_FFELT:
1357 : case t_COMPLEX:
1358 : case t_PADIC:
1359 28 : case t_QUAD: return gcopy(x);
1360 8302 : case t_POL: return poldeflate(x,v,d);
1361 6818 : case t_SER: return serdeflate(x,v,d);
1362 7 : case t_POLMOD:
1363 7 : if (varncmp(varn(gel(x,1)), v) >= 0) return gcopy(x);
1364 : /* fall through */
1365 : case t_RFRAC:
1366 : case t_VEC:
1367 : case t_COL:
1368 14 : case t_MAT: return vdeflate(x,v,d);
1369 21 : case t_LIST: return listdeflate(x,v,d);
1370 : }
1371 0 : pari_err_TYPE("gdeflate",x);
1372 : return NULL; /* LCOV_EXCL_LINE */
1373 : }
1374 :
1375 : /* set *m to the largest d such that x0 = A(X^d); return A */
1376 : GEN
1377 730396 : RgX_deflate_max(GEN x, long *m)
1378 : {
1379 730396 : *m = RgX_deflate_order(x);
1380 730396 : return RgX_deflate(x, *m);
1381 : }
1382 : GEN
1383 330313 : ZX_deflate_max(GEN x, long *m)
1384 : {
1385 330313 : *m = ZX_deflate_order(x);
1386 330313 : return RgX_deflate(x, *m);
1387 : }
1388 :
1389 : static int
1390 25081 : serequalXk(GEN x)
1391 : {
1392 25081 : long i, l = lg(x);
1393 25081 : if (l == 2 || !isint1(gel(x,2))) return 0;
1394 11025 : for (i = 3; i < l; i++)
1395 8694 : if (!isintzero(gel(x,i))) return 0;
1396 2331 : return 1;
1397 : }
1398 :
1399 : static GEN
1400 84 : gsubst_v(GEN e, long v, GEN x)
1401 245 : { pari_APPLY_same(gsubst(e, v, gel(x,i))); }
1402 :
1403 : static GEN
1404 14 : constmat(GEN z, long n)
1405 : {
1406 14 : GEN y = cgetg(n+1, t_MAT), c = const_col(n, gcopy(z));
1407 : long i;
1408 35 : for (i = 1; i <= n; i++) gel(y, i) = c;
1409 14 : return y;
1410 : }
1411 : static GEN
1412 56 : scalarmat2(GEN o, GEN z, long n)
1413 : {
1414 : GEN y;
1415 : long i;
1416 56 : if (n == 0) return cgetg(1, t_MAT);
1417 56 : if (n == 1) retmkmat(mkcol(gcopy(o)));
1418 35 : y = cgetg(n+1, t_MAT); z = gcopy(z); o = gcopy(o);
1419 105 : for (i = 1; i <= n; i++) { gel(y, i) = const_col(n, z); gcoeff(y,i,i) = o; }
1420 35 : return y;
1421 : }
1422 : /* x * y^0, n = dim(y) if t_MAT, else -1 */
1423 : static GEN
1424 730968 : subst_higher(GEN x, GEN y, long n)
1425 : {
1426 730968 : GEN o = Rg_get_1(y);
1427 730968 : if (o == gen_1) return n < 0? gcopy(x): scalarmat(x,n);
1428 98 : x = gmul(x,o); return n < 0? x: scalarmat2(x, Rg_get_0(y), n);
1429 : }
1430 :
1431 : /* x t_POLMOD, v strictly lower priority than var(x) */
1432 : static GEN
1433 574 : subst_polmod(GEN x, long v, GEN y)
1434 : {
1435 : long l, i;
1436 574 : GEN a = gsubst(gel(x,2),v,y), b = gsubst(gel(x,1),v,y), z;
1437 :
1438 574 : if (typ(b) != t_POL) pari_err_TYPE2("substitution",x,y);
1439 574 : if (typ(a) != t_POL || varncmp(varn(a), varn(b)) >= 0) return gmodulo(a, b);
1440 539 : l = lg(a); z = cgetg(l,t_POL); z[1] = a[1];
1441 4151 : for (i = 2; i < l; i++) gel(z,i) = gmodulo(gel(a,i),b);
1442 539 : return normalizepol_lg(z, l);
1443 : }
1444 : /* Trunc to n terms; x + O(t^(n + v(x))). FIXME: export ? */
1445 : static GEN
1446 77 : sertrunc(GEN x, long n)
1447 : {
1448 77 : long i, l = n + 2;
1449 : GEN y;
1450 77 : if (l >= lg(x)) return x;
1451 14 : if (n <= 0) return zeroser(varn(x), n + valser(x));
1452 14 : y = cgetg(l, t_SER);
1453 28 : for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
1454 14 : y[1] = x[1]; return y;
1455 : }
1456 : /* FIXME: export ? */
1457 : static GEN
1458 2338 : sertrunc_copy(GEN x, long n)
1459 : {
1460 2338 : long i, l = minss(n + 2, lg(x));
1461 2338 : GEN y = cgetg(l, t_SER);
1462 15337 : for (i = 2; i < l; i++) gel(y,i) = gcopy(gel(x,i));
1463 2338 : y[1] = x[1]; return y;
1464 : }
1465 :
1466 : GEN
1467 2741615 : gsubst(GEN x, long v, GEN y)
1468 : {
1469 2741615 : long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
1470 : long l, vx, vy, ex, ey, i, j, k, jb, matn;
1471 : pari_sp av, av2;
1472 : GEN X, t, z;
1473 :
1474 2741615 : switch(ty)
1475 : {
1476 84 : case t_VEC: case t_COL:
1477 84 : return gsubst_v(x, v, y);
1478 175 : case t_MAT:
1479 175 : if (ly==1) return cgetg(1,t_MAT);
1480 168 : if (ly == lgcols(y)) { matn = ly - 1; break; }
1481 : /* fall through */
1482 : case t_QFB:
1483 7 : pari_err_TYPE2("substitution",x,y);
1484 2741356 : default: matn = -1;
1485 : }
1486 2741517 : if (is_scalar_t(tx))
1487 : {
1488 383296 : if (tx == t_POLMOD && varncmp(v, varn(gel(x,1))) > 0)
1489 : {
1490 574 : av = avma;
1491 574 : return gc_upto(av, subst_polmod(x, v, y));
1492 : }
1493 382722 : return subst_higher(x, y, matn);
1494 : }
1495 :
1496 2358221 : switch(tx)
1497 : {
1498 2112432 : case t_POL:
1499 2112432 : vx = varn(x);
1500 2112432 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1501 2111109 : if (varncmp(vx, v) < 0)
1502 : {
1503 167770 : if (lx == 2) return zeropol(vx);
1504 167553 : av = avma;
1505 167553 : if (lx == 3)
1506 : {
1507 28581 : z = gsubst(gel(x,2), v, y);
1508 28581 : if (typ(z) == t_POL && varn(z) == vx) return z;
1509 28581 : return gc_upto(av, gmul(pol_1(vx), z));
1510 : }
1511 138972 : z = cgetg(lx-1, t_VEC);
1512 742776 : for (i = 2; i < lx; i++) gel(z,i-1) = gsubst(gel(x,i),v,y);
1513 138972 : return gc_upto(av, poleval(z, pol_x(vx)));
1514 : }
1515 : /* v = vx */
1516 1943339 : if (lx == 2)
1517 : {
1518 29729 : z = Rg_get_0(y);
1519 29729 : return matn >= 0? constmat(z, matn): z;
1520 : }
1521 1913610 : if (lx == 3)
1522 : {
1523 346916 : x = subst_higher(gel(x,2), y, matn);
1524 346916 : if (matn >= 0) return x;
1525 346902 : vy = gvar(y);
1526 346902 : return (vy == NO_VARIABLE)? x: gmul(x, pol_1(vy));
1527 : }
1528 1566694 : return matn >= 0? RgX_RgM_eval(x, y): poleval(x,y);
1529 :
1530 31164 : case t_SER:
1531 31164 : vx = varn(x);
1532 31164 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1533 31157 : ex = valser(x);
1534 31157 : if (varncmp(vx, v) < 0)
1535 : {
1536 56 : if (lx == 2) return matn >= 0? scalarmat(x, matn): gcopy(x);
1537 56 : av = avma; X = pol_x(vx);
1538 56 : av2 = avma;
1539 56 : z = gadd(gsubst(gel(x,lx-1),v,y), zeroser(vx,1));
1540 224 : for (i = lx-2; i>=2; i--)
1541 : {
1542 168 : z = gadd(gmul(z,X), gsubst(gel(x,i),v,y));
1543 168 : if (gc_needed(av2,1))
1544 : {
1545 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1546 0 : z = gc_upto(av2, z);
1547 : }
1548 : }
1549 56 : if (ex) z = gmul(z, pol_xnall(ex,vx));
1550 56 : return gc_upto(av, z);
1551 : }
1552 31101 : switch(ty) /* here vx == v */
1553 : {
1554 25214 : case t_SER:
1555 25214 : vy = varn(y); ey = valser(y);
1556 25214 : if (ey < 1 || lx == 2) return zeroser(vy, ey*(ex+lx-2));
1557 25207 : if (ey == 1 && serequalXk(y)
1558 2331 : && (varncmp(vx,vy) >= 0 || varncmp(gvar2(x), vy) >= 0))
1559 : { /* y = t + O(t^N) */
1560 2331 : if (lx > ly)
1561 : { /* correct number of significant terms */
1562 1995 : l = ly;
1563 1995 : if (!ex)
1564 1974 : for (i = 3; i < lx; i++)
1565 1974 : if (++l >= lx || !gequal0(gel(x,i))) break;
1566 1995 : lx = l;
1567 : }
1568 2331 : z = sertrunc_copy(x, lx - 2); if (vx != vy) setvarn(z,vy);
1569 2331 : return z;
1570 : }
1571 22876 : if (vy != vx)
1572 : {
1573 28 : long nx = lx - 2, n = minss(ey * nx, ly - 2);
1574 28 : av = avma; z = gel(x, nx+1);
1575 91 : for (i = nx; i > 1; i--)
1576 : {
1577 63 : z = gadd(gmul(y,z), gel(x,i));
1578 63 : if (gc_needed(av,1))
1579 : {
1580 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1581 0 : z = gc_upto(av, z);
1582 : }
1583 : }
1584 28 : if (ex)
1585 : {
1586 21 : if (ex < 0) { y = ginv(y); ex = -ex; }
1587 21 : z = gmul(z, gpowgs(sertrunc(y, n), ex));
1588 : }
1589 28 : if (lg(z)-2 > n) z = sertrunc_copy(z, n);
1590 28 : return gc_upto(av,z);
1591 : }
1592 22848 : l = (lx-2)*ey + 2;
1593 22848 : if (ex) { if (l>ly) l = ly; }
1594 22792 : else if (lx != 3)
1595 : {
1596 22820 : for (i = 3; i < lx; i++)
1597 22806 : if (!isexactzero(gel(x,i))) break;
1598 22792 : l = minss(l, (i-2)*ey + (gequal0(y)? 2 : ly));
1599 : }
1600 22848 : av = avma; t = leafcopy(y);
1601 22848 : if (l < ly) setlg(t, l);
1602 22848 : z = scalarser(gen_1, varn(y), l-2);
1603 22848 : gel(z,2) = gel(x,2); /* ensure lg(z) = l even if x[2] = 0 */
1604 90790 : for (i = 3, jb = ey; jb <= l-2; i++,jb += ey)
1605 : {
1606 67949 : if (i < lx) {
1607 151676 : for (j = jb+2; j < minss(l, jb+ly); j++)
1608 83818 : gel(z,j) = gadd(gel(z,j), gmul(gel(x,i),gel(t,j-jb)));
1609 : }
1610 106232 : for (j = minss(ly-1, l-1-jb-ey); j > 1; j--)
1611 : {
1612 38290 : GEN a = gmul(gel(t,2), gel(y,j));
1613 87479 : for (k=2; k<j; k++) a = gadd(a, gmul(gel(t,j-k+2), gel(y,k)));
1614 38290 : gel(t,j) = a;
1615 : }
1616 67942 : if (gc_needed(av,1))
1617 : {
1618 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst");
1619 0 : (void)gc_all(av,2, &z,&t);
1620 : }
1621 : }
1622 22841 : if (!ex) return gc_GEN(av,z);
1623 56 : if (ex < 0) { ex = -ex; y = ginv(y); }
1624 56 : return gc_upto(av, gmul(z, gpowgs(sertrunc(y, l-2), ex)));
1625 :
1626 5845 : case t_POL: case t_RFRAC:
1627 : {
1628 5845 : long N, n = lx-2;
1629 5845 : vy = gvar(y); ey = gval(y,vy);
1630 5845 : if (ey == LONG_MAX)
1631 : { /* y = 0 */
1632 49 : if (ex < 0) pari_err_INV("gsubst",y);
1633 35 : if (!n) return gcopy(x);
1634 28 : if (ex > 0) return Rg_get_0(ty == t_RFRAC? gel(y,2): y);
1635 14 : y = Rg_get_1(ty == t_RFRAC? gel(y,2): y);
1636 14 : return gmul(y, gel(x,2));
1637 : }
1638 5796 : if (ey < 1 || n == 0) return zeroser(vy, ey*(ex+n));
1639 5789 : av = avma;
1640 5789 : n *= ey;
1641 5789 : N = ex? n: maxss(n-ey,1);
1642 5789 : y = (ty == t_RFRAC)? rfrac_to_ser_i(y, N+2): RgX_to_ser(y, N+2);
1643 5789 : if (lg(y)-2 > n) setlg(y, n+2);
1644 5789 : x = ser2pol_i(x, lx);
1645 5789 : if (varncmp(vy,vx) > 0)
1646 49 : z = gadd(poleval(x, y), zeroser(vy,n));
1647 : else
1648 : {
1649 5740 : z = RgXn_eval(x, ser2rfrac_i(y), n);
1650 5740 : if (varn(z) == vy) z = RgX_to_ser(z, n+2);
1651 0 : else z = scalarser(z, vy, n);
1652 : }
1653 5789 : if (!ex) return gc_GEN(av, z);
1654 5684 : return gc_upto(av, gmul(z, gpowgs(y,ex)));
1655 : }
1656 :
1657 42 : default:
1658 42 : if (isexactzero(y))
1659 : {
1660 35 : if (ex < 0) pari_err_INV("gsubst",y);
1661 14 : if (ex > 0) return gcopy(y);
1662 7 : if (lx > 2) return gadd(gel(x,2), y); /*add maps to correct ring*/
1663 : }
1664 7 : pari_err_TYPE2("substitution",x,y);
1665 : }
1666 0 : break;
1667 :
1668 1904 : case t_RFRAC:
1669 : {
1670 1904 : GEN a = gel(x,1), b = gel(x,2);
1671 1904 : av = avma;
1672 1904 : a = gsubst(a, v, y);
1673 1904 : b = gsubst(b, v, y); return gc_upto(av, gdiv(a, b));
1674 : }
1675 :
1676 666842 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gsubst(gel(x,i),v,y));
1677 56 : case t_LIST:
1678 56 : z = mklist();
1679 56 : list_data(z) = list_data(x)? gsubst(list_data(x),v,y): NULL;
1680 56 : return z;
1681 : }
1682 0 : return gcopy(x);
1683 : }
1684 :
1685 : /* Return P(x * h), not memory clean */
1686 : GEN
1687 4193 : ser_unscale(GEN P, GEN h)
1688 : {
1689 4193 : long l = lg(P);
1690 4193 : GEN Q = cgetg(l,t_SER);
1691 4193 : Q[1] = P[1];
1692 4193 : if (l != 2)
1693 : {
1694 4193 : long i = 2;
1695 4193 : GEN hi = gpowgs(h, valser(P));
1696 4193 : gel(Q,i) = gmul(gel(P,i), hi);
1697 200508 : for (i++; i<l; i++) { hi = gmul(hi,h); gel(Q,i) = gmul(gel(P,i), hi); }
1698 : }
1699 4193 : return Q;
1700 : }
1701 :
1702 : static int
1703 1358 : safe_polmod(GEN r)
1704 : {
1705 1358 : GEN a = gel(r,1), b = gel(r,2);
1706 1358 : long t = typ(a);
1707 1358 : return 0;
1708 : if (gvar2(b) != NO_VARIABLE) return 0;
1709 : if (is_scalar_t(t)) return 1;
1710 : return (t == t_POL && varn(a) == varn(b) && gvar2(a) == NO_VARIABLE);
1711 : }
1712 : GEN
1713 980 : gsubstvec(GEN e, GEN v, GEN r)
1714 : {
1715 980 : pari_sp av = avma;
1716 980 : long i, j, k, l = lg(v);
1717 : GEN w, z, R;
1718 980 : if ( !is_vec_t(typ(v)) ) pari_err_TYPE("substvec",v);
1719 980 : if ( !is_vec_t(typ(r)) ) pari_err_TYPE("substvec",r);
1720 980 : if (lg(r)!=l) pari_err_DIM("substvec");
1721 980 : w = cgetg(l, t_VECSMALL);
1722 980 : z = cgetg(l, t_VECSMALL);
1723 980 : R = cgetg(l, t_VEC); k = 0;
1724 4298 : for(i = j = 1; i < l; i++)
1725 : {
1726 3318 : GEN T = gel(v,i), ri = gel(r,i);
1727 3318 : if (!gequalX(T)) pari_err_TYPE("substvec [not a variable]", T);
1728 3318 : if (gvar(ri) == NO_VARIABLE || (typ(ri) == t_POLMOD && safe_polmod(ri)))
1729 : { /* no need to take precautions */
1730 1855 : e = gsubst(e, varn(T), ri);
1731 1855 : if (is_vec_t(typ(ri)) && k++) e = shallowconcat1(e);
1732 : }
1733 : else
1734 : {
1735 1463 : w[j] = varn(T);
1736 1463 : z[j] = fetch_var_higher();
1737 1463 : gel(R,j) = ri; j++;
1738 : }
1739 : }
1740 2443 : for(i = 1; i < j; i++) e = gsubst(e,w[i],pol_x(z[i]));
1741 2443 : for(i = 1; i < j; i++)
1742 : {
1743 1463 : e = gsubst(e,z[i],gel(R,i));
1744 1463 : if (is_vec_t(typ(gel(R,i))) && k++) e = shallowconcat1(e);
1745 : }
1746 2443 : for(i = 1; i < j; i++) (void)delete_var();
1747 980 : return k > 1? gc_GEN(av, e): gc_upto(av, e);
1748 : }
1749 :
1750 : /*******************************************************************/
1751 : /* */
1752 : /* SERIE RECIPROQUE D'UNE SERIE */
1753 : /* */
1754 : /*******************************************************************/
1755 :
1756 : GEN
1757 98 : serreverse(GEN x)
1758 : {
1759 98 : long v=varn(x), lx = lg(x), i, mi;
1760 98 : pari_sp av0 = avma, av;
1761 : GEN a, y, u;
1762 :
1763 98 : if (typ(x)!=t_SER) pari_err_TYPE("serreverse",x);
1764 98 : if (valser(x)!=1) pari_err_DOMAIN("serreverse", "valuation", "!=", gen_1,x);
1765 91 : if (lx < 3) pari_err_DOMAIN("serreverse", "x", "=", gen_0,x);
1766 91 : y = ser_normalize(x);
1767 91 : if (y == x) a = NULL; else { a = gel(x,2); x = y; }
1768 91 : av = avma;
1769 252 : mi = lx-1; while (mi>=3 && gequal0(gel(x,mi))) mi--;
1770 91 : u = cgetg(lx,t_SER);
1771 91 : y = cgetg(lx,t_SER);
1772 91 : u[1] = y[1] = evalsigne(1) | _evalvalser(1) | evalvarn(v);
1773 91 : gel(u,2) = gel(y,2) = gen_1;
1774 91 : if (lx > 3)
1775 : {
1776 84 : gel(u,3) = gmulsg(-2,gel(x,3));
1777 84 : gel(y,3) = gneg(gel(x,3));
1778 : }
1779 1113 : for (i=3; i<lx-1; )
1780 : {
1781 : pari_sp av2;
1782 : GEN p1;
1783 1022 : long j, k, K = minss(i,mi);
1784 8456 : for (j=3; j<i+1; j++)
1785 : {
1786 7434 : av2 = avma; p1 = gel(x,j);
1787 39291 : for (k = maxss(3,j+2-mi); k < j; k++)
1788 31857 : p1 = gadd(p1, gmul(gel(u,k),gel(x,j-k+2)));
1789 7434 : p1 = gneg(p1);
1790 7434 : gel(u,j) = gc_upto(av2, gadd(gel(u,j), p1));
1791 : }
1792 1022 : av2 = avma;
1793 1022 : p1 = gmulsg(i,gel(x,i+1));
1794 8309 : for (k = 2; k < K; k++)
1795 : {
1796 7287 : GEN p2 = gmul(gel(x,k+1),gel(u,i-k+2));
1797 7287 : p1 = gadd(p1, gmulsg(k,p2));
1798 : }
1799 1022 : i++;
1800 1022 : gel(u,i) = gc_upto(av2, gneg(p1));
1801 1022 : gel(y,i) = gdivgu(gel(u,i), i-1);
1802 1022 : if (gc_needed(av,2))
1803 : {
1804 0 : GEN dummy = cgetg(1,t_VEC);
1805 0 : if(DEBUGMEM>1) pari_warn(warnmem,"serreverse");
1806 0 : for(k=i+1; k<lx; k++) gel(u,k) = gel(y,k) = dummy;
1807 0 : (void)gc_all(av,2, &u,&y);
1808 : }
1809 : }
1810 91 : if (a) y = ser_unscale(y, ginv(a));
1811 91 : return gc_GEN(av0,y);
1812 : }
1813 :
1814 : /*******************************************************************/
1815 : /* */
1816 : /* DERIVATION ET INTEGRATION */
1817 : /* */
1818 : /*******************************************************************/
1819 : GEN
1820 29344 : derivser(GEN x)
1821 : {
1822 29344 : long i, vx = varn(x), e = valser(x), lx = lg(x);
1823 : GEN y;
1824 29344 : if (ser_isexactzero(x))
1825 : {
1826 7 : x = gcopy(x);
1827 7 : if (e) setvalser(x,e-1);
1828 7 : return x;
1829 : }
1830 29337 : if (e)
1831 : {
1832 602 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|evalvalser(e-1) | evalvarn(vx);
1833 22960 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i+e-2,gel(x,i));
1834 : } else {
1835 28735 : if (lx == 3) return zeroser(vx, 0);
1836 24815 : lx--;
1837 24815 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|_evalvalser(0) | evalvarn(vx);
1838 78652 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1839 : }
1840 25417 : return normalizeser(y);
1841 : }
1842 :
1843 : static GEN
1844 56 : rfrac_deriv(GEN x, long v)
1845 : {
1846 56 : pari_sp av = avma;
1847 56 : GEN y = cgetg(3,t_RFRAC), a = gel(x,1), b = gel(x,2), bp, b0, t, T;
1848 56 : long vx = varn(b);
1849 :
1850 56 : bp = deriv(b, v);
1851 56 : t = simplify_shallow(RgX_gcd(bp, b));
1852 56 : if (typ(t) != t_POL || varn(t) != vx)
1853 : {
1854 35 : if (gequal1(t)) b0 = b;
1855 : else
1856 : {
1857 0 : b0 = RgX_Rg_div(b, t);
1858 0 : bp = RgX_Rg_div(bp, t);
1859 : }
1860 35 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1861 35 : if (isexactzero(a)) return gc_upto(av, a);
1862 35 : if (b0 == b)
1863 : {
1864 35 : gel(y,1) = gc_upto((pari_sp)y, a);
1865 35 : gel(y,2) = RgX_sqr(b);
1866 : }
1867 : else
1868 : {
1869 0 : gel(y,1) = a;
1870 0 : gel(y,2) = RgX_Rg_mul(RgX_sqr(b0), t);
1871 0 : y = gc_GEN(av, y);
1872 : }
1873 35 : return y;
1874 : }
1875 21 : b0 = gdivexact(b, t);
1876 21 : bp = gdivexact(bp,t);
1877 21 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1878 21 : if (isexactzero(a)) return gc_upto(av, a);
1879 14 : T = RgX_gcd(a, t);
1880 14 : if (typ(T) != t_POL || varn(T) != vx)
1881 : {
1882 0 : a = gdiv(a, T);
1883 0 : t = gdiv(t, T);
1884 : }
1885 14 : else if (!gequal1(T))
1886 : {
1887 0 : a = gdivexact(a, T);
1888 0 : t = gdivexact(t, T);
1889 : }
1890 14 : gel(y,1) = a;
1891 14 : gel(y,2) = gmul(RgX_sqr(b0), t);
1892 14 : return gc_GEN(av, y);
1893 : }
1894 :
1895 : GEN
1896 112910 : deriv(GEN x, long v)
1897 : {
1898 112910 : long tx = typ(x);
1899 112910 : if (is_const_t(tx))
1900 39424 : switch(tx)
1901 : {
1902 14 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
1903 14 : case t_FFELT: return FF_zero(x);
1904 39396 : default: return gen_0;
1905 : }
1906 73486 : if (v < 0)
1907 : {
1908 49 : if (tx == t_CLOSURE) return closure_deriv(x);
1909 49 : v = gvar9(x);
1910 : }
1911 73486 : switch(tx)
1912 : {
1913 14 : case t_POLMOD:
1914 : {
1915 14 : GEN a = gel(x,2), b = gel(x,1);
1916 14 : if (v == varn(b)) return Rg_get_0(b);
1917 7 : retmkpolmod(deriv(a,v), RgX_copy(b));
1918 : }
1919 73227 : case t_POL:
1920 73227 : switch(varncmp(varn(x), v))
1921 : {
1922 0 : case 1: return Rg_get_0(x);
1923 65702 : case 0: return RgX_deriv(x);
1924 : }
1925 112287 : pari_APPLY_pol(deriv(gel(x,i),v));
1926 :
1927 147 : case t_SER:
1928 147 : switch(varncmp(varn(x), v))
1929 : {
1930 0 : case 1: return Rg_get_0(x);
1931 133 : case 0: return derivser(x);
1932 : }
1933 14 : if (ser_isexactzero(x)) return gcopy(x);
1934 28 : pari_APPLY_ser(deriv(gel(x,i),v));
1935 :
1936 56 : case t_RFRAC:
1937 56 : return rfrac_deriv(x,v);
1938 :
1939 42 : case t_VEC: case t_COL: case t_MAT:
1940 84 : pari_APPLY_same(deriv(gel(x,i),v));
1941 : }
1942 0 : pari_err_TYPE("deriv",x);
1943 : return NULL; /* LCOV_EXCL_LINE */
1944 : }
1945 :
1946 : /* n-th derivative of t_SER x, n > 0 */
1947 : static GEN
1948 210 : derivnser(GEN x, long n)
1949 : {
1950 210 : long i, vx = varn(x), e = valser(x), lx = lg(x);
1951 : GEN y;
1952 210 : if (ser_isexactzero(x))
1953 : {
1954 7 : x = gcopy(x);
1955 7 : if (e) setvalser(x,e-n);
1956 7 : return x;
1957 : }
1958 203 : if (e < 0 || e >= n)
1959 : {
1960 161 : y = cgetg(lx,t_SER);
1961 161 : y[1] = evalsigne(1)| evalvalser(e-n) | evalvarn(vx);
1962 756 : for (i=0; i<lx-2; i++)
1963 595 : gel(y,i+2) = gmul(muls_interval(i+e-n+1,i+e), gel(x,i+2));
1964 : } else {
1965 42 : if (lx <= n+2) return zeroser(vx, 0);
1966 42 : lx -= n;
1967 42 : y = cgetg(lx,t_SER);
1968 42 : y[1] = evalsigne(1)|_evalvalser(0) | evalvarn(vx);
1969 245 : for (i=0; i<lx-2; i++)
1970 203 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n-e));
1971 : }
1972 203 : return normalizeser(y);
1973 : }
1974 :
1975 : static GEN
1976 42 : rfrac_derivn(GEN x, long n, long vs)
1977 : {
1978 42 : pari_sp av = avma;
1979 42 : GEN u = gel(x,1), v = gel(x,2);
1980 42 : GEN dv = deriv(v, vs);
1981 : long i;
1982 112 : for (i=1; i<=n; i++)
1983 : {
1984 70 : GEN du = deriv(u, vs);
1985 70 : u = gadd(gmul(du, v), gmulsg (-i, gmul(dv, u)));
1986 : }
1987 42 : v = gpowgs(v, n+1);
1988 42 : return gc_upto(av, gdiv(u, v));
1989 : }
1990 :
1991 : GEN
1992 1379 : derivn(GEN x, long n, long v)
1993 : {
1994 : long tx;
1995 1379 : if (n < 0) pari_err_DOMAIN("derivn","n","<", gen_0, stoi(n));
1996 1372 : if (n == 0) return gcopy(x);
1997 1372 : tx = typ(x);
1998 1372 : if (is_const_t(tx))
1999 49 : switch(tx)
2000 : {
2001 21 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
2002 21 : case t_FFELT: return FF_zero(x);
2003 7 : default: return gen_0;
2004 : }
2005 1323 : if (v < 0)
2006 : {
2007 1085 : if (tx == t_CLOSURE) return closure_derivn(x, n);
2008 973 : v = gvar9(x);
2009 : }
2010 1211 : switch(tx)
2011 : {
2012 21 : case t_POLMOD:
2013 : {
2014 21 : GEN a = gel(x,2), b = gel(x,1);
2015 21 : if (v == varn(b)) return Rg_get_0(b);
2016 14 : retmkpolmod(derivn(a,n,v), RgX_copy(b));
2017 : }
2018 861 : case t_POL:
2019 861 : switch(varncmp(varn(x), v))
2020 : {
2021 0 : case 1: return Rg_get_0(x);
2022 833 : case 0: return RgX_derivn(x,n);
2023 : }
2024 84 : pari_APPLY_pol(derivn(gel(x,i),n,v));
2025 :
2026 217 : case t_SER:
2027 217 : switch(varncmp(varn(x), v))
2028 : {
2029 0 : case 1: return Rg_get_0(x);
2030 210 : case 0: return derivnser(x, n);
2031 : }
2032 7 : if (ser_isexactzero(x)) return gcopy(x);
2033 28 : pari_APPLY_ser(derivn(gel(x,i),n,v));
2034 :
2035 42 : case t_RFRAC:
2036 42 : return rfrac_derivn(x, n, v);
2037 :
2038 63 : case t_VEC: case t_COL: case t_MAT:
2039 126 : pari_APPLY_same(derivn(gel(x,i),n,v));
2040 : }
2041 7 : pari_err_TYPE("derivn",x);
2042 : return NULL; /* LCOV_EXCL_LINE */
2043 : }
2044 :
2045 : static long
2046 833 : lookup(GEN v, long vx)
2047 : {
2048 833 : long i, l = lg(v);
2049 1491 : for(i=1; i<l; i++)
2050 1253 : if (varn(gel(v,i)) == vx) return i;
2051 238 : return 0;
2052 : }
2053 :
2054 : static GEN
2055 119 : serdiffop(GEN x, GEN v, GEN dv) { pari_APPLY_ser(diffop(gel(x,i),v,dv)); }
2056 : GEN
2057 3535 : diffop(GEN x, GEN v, GEN dv)
2058 : {
2059 : pari_sp av;
2060 3535 : long i, idx, lx, tx = typ(x), vx;
2061 : GEN y;
2062 3535 : if (!is_vec_t(typ(v))) pari_err_TYPE("diffop",v);
2063 3535 : if (!is_vec_t(typ(dv))) pari_err_TYPE("diffop",dv);
2064 3535 : if (lg(v)!=lg(dv)) pari_err_DIM("diffop");
2065 3535 : if (is_const_t(tx)) return gen_0;
2066 1148 : switch(tx)
2067 : {
2068 84 : case t_POLMOD:
2069 84 : av = avma;
2070 84 : vx = varn(gel(x,1)); idx = lookup(v,vx);
2071 84 : if (idx) /*Assume the users know what they are doing */
2072 0 : y = gmodulo(diffop(gel(x,2),v,dv), gel(x,1));
2073 : else
2074 : {
2075 84 : GEN m = gel(x,1), pol=gel(x,2);
2076 84 : GEN u = gneg(gdiv(diffop(m,v,dv),RgX_deriv(m)));
2077 84 : y = diffop(pol,v,dv);
2078 84 : if (typ(pol)==t_POL && varn(pol)==varn(m))
2079 70 : y = gadd(y, gmul(u,RgX_deriv(pol)));
2080 84 : y = gmodulo(y, gel(x,1));
2081 : }
2082 84 : return gc_upto(av, y);
2083 952 : case t_POL:
2084 952 : if (signe(x)==0) return gen_0;
2085 742 : vx = varn(x); idx = lookup(v,vx);
2086 742 : av = avma; lx = lg(x);
2087 742 : y = diffop(gel(x,lx-1),v,dv);
2088 2842 : for (i=lx-2; i>=2; i--) y = gadd(gmul(y,pol_x(vx)),diffop(gel(x,i),v,dv));
2089 742 : if (idx) y = gadd(y, gmul(gel(dv,idx),RgX_deriv(x)));
2090 742 : return gc_upto(av, y);
2091 :
2092 7 : case t_SER:
2093 7 : if (signe(x)==0) return gen_0;
2094 7 : vx = varn(x); idx = lookup(v,vx);
2095 7 : if (!idx) return gen_0;
2096 7 : av = avma;
2097 7 : if (ser_isexactzero(x)) y = x;
2098 : else
2099 : {
2100 7 : y = serdiffop(x, v, dv); /* y is probably invalid */
2101 7 : y = gsubst(y, vx, pol_x(vx)); /* Fix that */
2102 : }
2103 7 : y = gadd(y, gmul(gel(dv,idx), derivser(x)));
2104 7 : return gc_upto(av, y);
2105 :
2106 105 : case t_RFRAC: {
2107 105 : GEN a = gel(x,1), b = gel(x,2), ap, bp;
2108 105 : av = avma;
2109 105 : ap = diffop(a, v, dv); bp = diffop(b, v, dv);
2110 105 : y = gsub(gdiv(ap,b),gdiv(gmul(a,bp),gsqr(b)));
2111 105 : return gc_upto(av, y);
2112 : }
2113 :
2114 0 : case t_VEC: case t_COL: case t_MAT:
2115 0 : pari_APPLY_same(diffop(gel(x,i),v,dv));
2116 : }
2117 0 : pari_err_TYPE("diffop",x);
2118 : return NULL; /* LCOV_EXCL_LINE */
2119 : }
2120 :
2121 : GEN
2122 42 : diffop0(GEN x, GEN v, GEN dv, long n)
2123 : {
2124 42 : pari_sp av=avma;
2125 : long i;
2126 245 : for(i=1; i<=n; i++)
2127 203 : x = gc_upto(av, diffop(x,v,dv));
2128 42 : return x;
2129 : }
2130 :
2131 : /********************************************************************/
2132 : /** **/
2133 : /** TAYLOR SERIES **/
2134 : /** **/
2135 : /********************************************************************/
2136 : /* swap vars (vx,v) in x (assume vx < v, vx main variable in x), then call
2137 : * act(data, v, x). FIXME: use in other places */
2138 : static GEN
2139 28 : swapvar_act(GEN x, long vx, long v, GEN (*act)(void*, long, GEN), void *data)
2140 : {
2141 28 : long v0 = fetch_var();
2142 28 : GEN y = act(data, v, gsubst(x,vx,pol_x(v0)));
2143 21 : y = gsubst(y,v0,pol_x(vx));
2144 21 : (void)delete_var(); return y;
2145 : }
2146 : /* x + O(v^data) */
2147 : static GEN
2148 14 : tayl_act(void *data, long v, GEN x) { return gadd(zeroser(v, (long)data), x); }
2149 : static GEN
2150 14 : integ_act(void *data, long v, GEN x) { (void)data; return integ(x,v); }
2151 :
2152 : GEN
2153 21 : tayl(GEN x, long v, long precS)
2154 : {
2155 21 : long vx = gvar9(x);
2156 : pari_sp av;
2157 :
2158 21 : if (varncmp(v, vx) <= 0) return gadd(zeroser(v,precS), x);
2159 14 : av = avma;
2160 14 : return gc_upto(av, swapvar_act(x, vx, v, tayl_act, (void*)precS));
2161 : }
2162 :
2163 : GEN
2164 7210 : ggrando(GEN x, long n)
2165 : {
2166 : long m, v;
2167 :
2168 7210 : switch(typ(x))
2169 : {
2170 4130 : case t_INT:/* bug 3 + O(1) */
2171 4130 : if (signe(x) <= 0) pari_err_DOMAIN("O", "x", "<=", gen_0, x);
2172 4130 : if (!is_pm1(x)) return zeropadic(x,n);
2173 : /* +/-1 = x^0 */
2174 91 : v = m = 0; break;
2175 3073 : case t_POL:
2176 3073 : if (!signe(x)) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2177 3073 : v = varn(x);
2178 3073 : m = n * RgX_val(x); break;
2179 7 : case t_RFRAC:
2180 7 : if (gequal0(gel(x,1))) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2181 7 : v = gvar(x);
2182 7 : m = n * gval(x,v); break;
2183 0 : default: pari_err_TYPE("O", x);
2184 : v = m = 0; /* LCOV_EXCL_LINE */
2185 : }
2186 3171 : return zeroser(v,m);
2187 : }
2188 :
2189 : /*******************************************************************/
2190 : /* */
2191 : /* FORMAL INTEGRATION */
2192 : /* */
2193 : /*******************************************************************/
2194 : GEN
2195 105 : RgX_integ(GEN x)
2196 : {
2197 105 : long i, lx = lg(x);
2198 : GEN y;
2199 105 : if (lx == 2) return RgX_copy(x);
2200 91 : y = cgetg(lx+1, t_POL); y[1] = x[1]; gel(y,2) = gen_0;
2201 273 : for (i=3; i<=lx; i++) gel(y,i) = gdivgu(gel(x,i-1),i-2);
2202 91 : return y;
2203 : }
2204 :
2205 : static void
2206 35 : err_intformal(GEN x)
2207 35 : { pari_err_DOMAIN("intformal", "residue(series, pole)", "!=", gen_0, x); }
2208 :
2209 : GEN
2210 29988 : integser(GEN x)
2211 : {
2212 29988 : long i, lx = lg(x), vx = varn(x), e = valser(x);
2213 : GEN y;
2214 29988 : if (lx == 2) return zeroser(vx, e+1);
2215 26047 : y = cgetg(lx, t_SER);
2216 110698 : for (i=2; i<lx; i++)
2217 : {
2218 84658 : long j = i+e-1;
2219 84658 : GEN c = gel(x,i);
2220 84658 : if (j)
2221 84343 : c = gdivgs(c, j);
2222 : else
2223 : { /* should be isexactzero, but try to avoid error */
2224 315 : if (!gequal0(c)) err_intformal(x);
2225 308 : c = gen_0;
2226 : }
2227 84651 : gel(y,i) = c;
2228 : }
2229 26040 : y[1] = evalsigne(1) | evalvarn(vx) | evalvalser(e+1); return y;
2230 : }
2231 :
2232 : GEN
2233 350 : integ(GEN x, long v)
2234 : {
2235 350 : long tx = typ(x), vx, n;
2236 350 : pari_sp av = avma;
2237 : GEN y, p1;
2238 :
2239 350 : if (v < 0) { v = gvar9(x); if (v == NO_VARIABLE) v = 0; }
2240 350 : if (is_scalar_t(tx))
2241 : {
2242 63 : if (tx == t_POLMOD)
2243 : {
2244 14 : GEN a = gel(x,2), b = gel(x,1);
2245 14 : vx = varn(b);
2246 14 : if (varncmp(v, vx) > 0) retmkpolmod(integ(a,v), RgX_copy(b));
2247 7 : if (v == vx) pari_err_PRIORITY("intformal",x,"=",v);
2248 : }
2249 49 : return deg1pol(x, gen_0, v);
2250 : }
2251 :
2252 287 : switch(tx)
2253 : {
2254 112 : case t_POL:
2255 112 : vx = varn(x);
2256 112 : if (v == vx) return RgX_integ(x);
2257 42 : if (lg(x) == 2) {
2258 14 : if (varncmp(vx, v) < 0) v = vx;
2259 14 : return zeropol(v);
2260 : }
2261 28 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2262 84 : pari_APPLY_pol(integ(gel(x,i),v));
2263 :
2264 77 : case t_SER:
2265 77 : vx = varn(x);
2266 77 : if (v == vx) return integser(x);
2267 21 : if (lg(x) == 2) {
2268 14 : if (varncmp(vx, v) < 0) v = vx;
2269 14 : return zeroser(v, valser(x));
2270 : }
2271 7 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2272 28 : pari_APPLY_ser(integ(gel(x,i),v));
2273 :
2274 56 : case t_RFRAC:
2275 : {
2276 56 : GEN a = gel(x,1), b = gel(x,2), c, d, s;
2277 56 : vx = varn(b);
2278 56 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2279 49 : if (varncmp(vx, v) < 0)
2280 14 : return gc_upto(av, swapvar_act(x, vx, v, integ_act, NULL));
2281 :
2282 35 : n = degpol(b);
2283 35 : if (typ(a) == t_POL && varn(a) == vx) n += degpol(a);
2284 35 : y = integ(gadd(x, zeroser(v,n + 2)), v);
2285 35 : y = gdiv(gtrunc(gmul(b, y)), b);
2286 35 : if (typ(y) != t_RFRAC) pari_err_BUG("intformal(t_RFRAC)");
2287 35 : c = gel(y,1); d = gel(y,2);
2288 35 : s = gsub(gmul(deriv(c,v),d), gmul(c,deriv(d,v)));
2289 : /* (c'd-cd')/d^2 = y' = x = a/b ? */
2290 35 : if (!gequal(gmul(s,b), gmul(a,gsqr(d)))) err_intformal(x);
2291 7 : if (typ(y)==t_RFRAC && lg(gel(y,1)) == lg(gel(y,2)))
2292 : {
2293 7 : GEN p2 = leading_coeff(gel(y,2));
2294 7 : p1 = gel(y,1);
2295 7 : if (typ(p1) == t_POL && varn(p1) == vx) p1 = leading_coeff(p1);
2296 7 : y = gsub(y, gdiv(p1,p2));
2297 : }
2298 7 : return gc_upto(av,y);
2299 : }
2300 :
2301 42 : case t_VEC: case t_COL: case t_MAT:
2302 84 : pari_APPLY_same(integ(gel(x,i),v));
2303 : }
2304 0 : pari_err_TYPE("integ",x);
2305 : return NULL; /* LCOV_EXCL_LINE */
2306 : }
2307 :
2308 : /*******************************************************************/
2309 : /* */
2310 : /* FLOOR */
2311 : /* */
2312 : /*******************************************************************/
2313 : static GEN
2314 742 : quad_floor(GEN x)
2315 : {
2316 742 : GEN Q = gel(x,1), D = quad_disc(x), u, v, b, d, z;
2317 742 : if (signe(D) < 0) return NULL;
2318 714 : x = Q_remove_denom(x, &d);
2319 714 : u = gel(x,2);
2320 714 : v = gel(x,3); b = gel(Q,3);
2321 714 : if (typ(u) != t_INT || typ(v) != t_INT) return NULL;
2322 : /* x0 = (2u + v*(-b + sqrt(D))) / (2d) */
2323 707 : z = sqrti(mulii(D, sqri(v)));
2324 707 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2325 : /* z = floor(v * sqrt(D)) */
2326 707 : z = addii(subii(shifti(u,1), mulii(v,b)), z);
2327 707 : return truedivii(z, d? shifti(d,1): gen_2);
2328 : }
2329 : GEN
2330 5539697 : gfloor(GEN x)
2331 : {
2332 5539697 : switch(typ(x))
2333 : {
2334 5306248 : case t_INT: return icopy(x);
2335 35 : case t_POL: return RgX_copy(x);
2336 214087 : case t_REAL: return floorr(x);
2337 18277 : case t_FRAC: return truedivii(gel(x,1),gel(x,2));
2338 735 : case t_QUAD:
2339 : {
2340 735 : pari_sp av = avma;
2341 : GEN y;
2342 735 : if (!(y = quad_floor(x))) break;
2343 700 : return gc_INT(av, y);
2344 : }
2345 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2346 98 : case t_VEC: case t_COL: case t_MAT:
2347 1533 : pari_APPLY_same(gfloor(gel(x,i)));
2348 : }
2349 238 : pari_err_TYPE("gfloor",x);
2350 : return NULL; /* LCOV_EXCL_LINE */
2351 : }
2352 :
2353 : GEN
2354 426531 : gfrac(GEN x)
2355 : {
2356 : pari_sp av;
2357 : GEN y;
2358 426531 : switch(typ(x))
2359 : {
2360 23868 : case t_INT: return gen_0;
2361 7 : case t_POL: return pol_0(varn(x));
2362 164416 : case t_REAL: av = avma; return gc_leaf(av, subri(x, floorr(x)));
2363 234533 : case t_FRAC: retmkfrac(modii(gel(x,1),gel(x,2)), icopy(gel(x,2)));
2364 7 : case t_QUAD:
2365 7 : av = avma; if (!(y = quad_floor(x))) break;
2366 7 : return gc_upto(av, gsub(x, y));
2367 7 : case t_RFRAC: retmkrfrac(grem(gel(x,1),gel(x,2)), gcopy(gel(x,2)));
2368 3665 : case t_VEC: case t_COL: case t_MAT:
2369 34174 : pari_APPLY_same(gfrac(gel(x,i)));
2370 : }
2371 28 : pari_err_TYPE("gfrac",x);
2372 : return NULL; /* LCOV_EXCL_LINE */
2373 : }
2374 :
2375 : /* assume x t_REAL */
2376 : GEN
2377 2721 : ceilr(GEN x)
2378 : {
2379 2721 : pari_sp av = avma;
2380 2721 : GEN y = floorr(x);
2381 2721 : if (cmpri(x, y)) return gc_INT(av, addui(1,y));
2382 29 : return y;
2383 : }
2384 :
2385 : GEN
2386 235307 : gceil(GEN x)
2387 : {
2388 : pari_sp av;
2389 : GEN y;
2390 235307 : switch(typ(x))
2391 : {
2392 18347 : case t_INT: return icopy(x);
2393 21 : case t_POL: return RgX_copy(x);
2394 2642 : case t_REAL: return ceilr(x);
2395 214171 : case t_FRAC:
2396 214171 : av = avma; y = divii(gel(x,1),gel(x,2));
2397 214171 : if (signe(gel(x,1)) > 0) y = gc_INT(av, addui(1,y));
2398 214171 : return y;
2399 49 : case t_QUAD:
2400 49 : if (!is_realquad(x)) break;
2401 42 : if (gequal0(gel(x,3))) return gceil(gel(x,2));
2402 35 : av = avma; return gc_upto(av, addiu(gfloor(x), 1));
2403 7 : case t_RFRAC:
2404 7 : return gdeuc(gel(x,1),gel(x,2));
2405 35 : case t_VEC: case t_COL: case t_MAT:
2406 105 : pari_APPLY_same(gceil(gel(x,i)));
2407 : }
2408 42 : pari_err_TYPE("gceil",x);
2409 : return NULL; /* LCOV_EXCL_LINE */
2410 : }
2411 :
2412 : GEN
2413 6237 : round0(GEN x, GEN *pte)
2414 : {
2415 6237 : if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
2416 6230 : return ground(x);
2417 : }
2418 :
2419 : /* x t_REAL, return q=floor(x+1/2), set e = expo(x-q) */
2420 : static GEN
2421 54676539 : round_i(GEN x, long *pe)
2422 : {
2423 : long e;
2424 54676539 : GEN B, q,r, m = mantissa_real(x, &e); /* x = m/2^e */
2425 54676539 : if (e <= 0)
2426 : {
2427 13773517 : if (e) m = shifti(m,-e);
2428 13773517 : if (pe) *pe = -e;
2429 13773517 : return m;
2430 : }
2431 40903022 : B = int2n(e-1);
2432 40903022 : m = addii(m, B);
2433 40903022 : q = shifti(m, -e);
2434 40903022 : r = remi2n(m, e);
2435 : /* 2^e (x+1/2) = m = 2^e q + r, sgn(r)=sgn(m), |r|<2^e */
2436 40903022 : if (!signe(r))
2437 32184 : { if (pe) *pe = -1; }
2438 : else
2439 : {
2440 40870838 : if (signe(m) < 0)
2441 : {
2442 16562759 : q = subiu(q,1);
2443 16562759 : r = addii(r, B);
2444 : }
2445 : else
2446 24308079 : r = subii(r, B);
2447 : /* |x - q| = |r| / 2^e */
2448 40870838 : if (pe) *pe = signe(r)? expi(r) - e: -e;
2449 40870838 : cgiv(r);
2450 : }
2451 40903022 : return q;
2452 : }
2453 : /* assume x a t_REAL */
2454 : GEN
2455 3635924 : roundr(GEN x)
2456 : {
2457 3635924 : long ex, s = signe(x);
2458 : pari_sp av;
2459 3635924 : if (!s || (ex=expo(x)) < -1) return gen_0;
2460 2761909 : if (ex == -1) return s>0? gen_1:
2461 297618 : absrnz_equal2n(x)? gen_0: gen_m1;
2462 1851600 : av = avma; x = round_i(x, &ex);
2463 1851600 : if (ex >= 0) pari_err_PREC( "roundr (precision loss in truncation)");
2464 1851600 : return gc_INT(av, x);
2465 : }
2466 : GEN
2467 314432 : roundr_safe(GEN x)
2468 : {
2469 314432 : long ex, s = signe(x);
2470 : pari_sp av;
2471 :
2472 314432 : if (!s || (ex = expo(x)) < -1) return gen_0;
2473 314388 : if (ex == -1) return s>0? gen_1:
2474 0 : absrnz_equal2n(x)? gen_0: gen_m1;
2475 314361 : av = avma; x = round_i(x, NULL);
2476 314361 : return gc_INT(av, x);
2477 : }
2478 :
2479 : GEN
2480 3424936 : ground(GEN x)
2481 : {
2482 : pari_sp av;
2483 : GEN y;
2484 :
2485 3424936 : switch(typ(x))
2486 : {
2487 836488 : case t_INT: return icopy(x);
2488 14 : case t_INTMOD: return gcopy(x);
2489 2032854 : case t_REAL: return roundr(x);
2490 50690 : case t_FRAC: return diviiround(gel(x,1), gel(x,2));
2491 49 : case t_QUAD:
2492 : {
2493 49 : GEN Q = gel(x,1), u, v, b, d, z;
2494 49 : av = avma;
2495 49 : if (is_realquad(x)) return gc_upto(av, gfloor(gadd(x, ghalf)));
2496 7 : u = gel(x,2);
2497 7 : v = gel(x,3); b = gel(Q,3);
2498 7 : u = ground(gsub(u, gmul2n(gmul(v,b),-1)));
2499 7 : v = Q_remove_denom(v, &d);
2500 7 : if (!d) d = gen_1;
2501 : /* Im x = v sqrt(|D|) / (2d),
2502 : * Im(round(x)) = floor((d + v sqrt(|D|)) / (2d))
2503 : * = floor(floor(d + v sqrt(|D|)) / (2d)) */
2504 7 : z = sqrti(mulii(sqri(v), quad_disc(x)));
2505 7 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2506 : /* z = floor(v * sqrt(|D|)) */
2507 7 : v = truedivii(addii(z, d), shifti(d,1));
2508 7 : return gc_GEN(av, signe(v)? mkcomplex(u,v): u);
2509 : }
2510 14 : case t_POLMOD:
2511 14 : retmkpolmod(ground(gel(x,2)), RgX_copy(gel(x,1)));
2512 :
2513 4096 : case t_COMPLEX:
2514 4096 : av = avma; y = cgetg(3, t_COMPLEX);
2515 4096 : gel(y,2) = ground(gel(x,2));
2516 4096 : if (!signe(gel(y,2))) { set_avma(av); return ground(gel(x,1)); }
2517 217 : gel(y,1) = ground(gel(x,1)); return y;
2518 :
2519 98 : case t_POL:
2520 4032 : pari_APPLY_pol(ground(gel(x,i)));
2521 182 : case t_SER:
2522 182 : if (ser_isexactzero(x)) return gcopy(x);
2523 42 : pari_APPLY_ser(ground(gel(x,i)));
2524 21 : case t_RFRAC:
2525 21 : av = avma;
2526 21 : return gc_upto(av, gdiv(ground(gel(x,1)), ground(gel(x,2))));
2527 500423 : case t_VEC: case t_COL: case t_MAT:
2528 2501659 : pari_APPLY_same(ground(gel(x,i)));
2529 : }
2530 7 : pari_err_TYPE("ground",x);
2531 : return NULL; /* LCOV_EXCL_LINE */
2532 : }
2533 :
2534 : /* e = number of error bits on integral part */
2535 : GEN
2536 92451965 : grndtoi(GEN x, long *e)
2537 : {
2538 : GEN y;
2539 : long i, lx, e1;
2540 : pari_sp av;
2541 :
2542 92451965 : if (e) *e = -(long)HIGHEXPOBIT;
2543 92451965 : switch(typ(x))
2544 : {
2545 10436191 : case t_INT: return icopy(x);
2546 57743649 : case t_REAL: {
2547 57743649 : long ex = expo(x);
2548 57743649 : if (!signe(x) || ex < -1)
2549 : {
2550 5233071 : if (e) *e = ex;
2551 5233071 : return gen_0;
2552 : }
2553 52510578 : av = avma; x = round_i(x, e);
2554 52510578 : return gc_INT(av, x);
2555 : }
2556 27863 : case t_FRAC:
2557 27863 : y = diviiround(gel(x,1), gel(x,2)); if (!e) return y;
2558 26505 : av = avma; *e = gexpo(gsub(x, y)); set_avma(av);
2559 26505 : return y;
2560 7 : case t_INTMOD: return gcopy(x);
2561 7 : case t_QUAD:
2562 7 : y = ground(x); av = avma; if (!e) return y;
2563 7 : *e = gexpo(gsub(x,y)); set_avma(avma); return y;
2564 1634890 : case t_COMPLEX:
2565 1634890 : av = avma; y = cgetg(3, t_COMPLEX);
2566 1634890 : gel(y,2) = grndtoi(gel(x,2), e);
2567 1634890 : if (!signe(gel(y,2))) {
2568 244855 : set_avma(av);
2569 244855 : y = grndtoi(gel(x,1), e? &e1: NULL);
2570 : }
2571 : else
2572 1390035 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL);
2573 1634890 : if (e && e1 > *e) *e = e1;
2574 1634890 : return y;
2575 :
2576 7 : case t_POLMOD:
2577 7 : retmkpolmod(grndtoi(gel(x,2), e), RgX_copy(gel(x,1)));
2578 :
2579 317699 : case t_POL:
2580 317699 : y = cgetg_copy(x, &lx); y[1] = x[1];
2581 3232957 : for (i=2; i<lx; i++)
2582 : {
2583 2915258 : gel(y,i) = grndtoi(gel(x,i), &e1);
2584 2915258 : if (e && e1 > *e) *e = e1;
2585 : }
2586 317699 : return normalizepol_lg(y, lx);
2587 168 : case t_SER:
2588 168 : if (ser_isexactzero(x)) return gcopy(x);
2589 154 : y = cgetg_copy(x, &lx); y[1] = x[1];
2590 441 : for (i=2; i<lx; i++)
2591 : {
2592 287 : gel(y,i) = grndtoi(gel(x,i), &e1);
2593 287 : if (e && e1 > *e) *e = e1;
2594 : }
2595 154 : return normalizeser(y);
2596 7 : case t_RFRAC:
2597 7 : y = cgetg(3,t_RFRAC);
2598 7 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2599 7 : gel(y,2) = grndtoi(gel(x,2), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2600 7 : return y;
2601 22291470 : case t_VEC: case t_COL: case t_MAT:
2602 22291470 : y = cgetg_copy(x, &lx);
2603 79692034 : for (i=1; i<lx; i++)
2604 : {
2605 57400564 : gel(y,i) = grndtoi(gel(x,i), e? &e1: NULL);
2606 57400564 : if (e && e1 > *e) *e = e1;
2607 : }
2608 22291470 : return y;
2609 : }
2610 7 : pari_err_TYPE("grndtoi",x);
2611 : return NULL; /* LCOV_EXCL_LINE */
2612 : }
2613 :
2614 : /* trunc(x * 2^s). lindep() sanity checks rely on this function to return a
2615 : * t_INT or fail when fed a non-t_COMPLEX input; so do not make this one
2616 : * recursive [ or change the lindep call ] */
2617 : GEN
2618 117762989 : gtrunc2n(GEN x, long s)
2619 : {
2620 : GEN z;
2621 117762989 : switch(typ(x))
2622 : {
2623 37620417 : case t_INT: return shifti(x, s);
2624 62858477 : case t_REAL: return trunc2nr(x, s);
2625 497 : case t_FRAC: {
2626 : pari_sp av;
2627 497 : GEN a = gel(x,1), b = gel(x,2), q;
2628 497 : if (s == 0) return divii(a, b);
2629 497 : av = avma;
2630 497 : if (s < 0) q = divii(shifti(a, s), b);
2631 : else {
2632 : GEN r;
2633 497 : q = dvmdii(a, b, &r);
2634 497 : q = addii(shifti(q,s), divii(shifti(r,s), b));
2635 : }
2636 497 : return gc_INT(av, q);
2637 : }
2638 17283598 : case t_COMPLEX:
2639 17283598 : z = cgetg(3, t_COMPLEX);
2640 17283598 : gel(z,2) = gtrunc2n(gel(x,2), s);
2641 17283598 : if (!signe(gel(z,2))) {
2642 1617489 : set_avma((pari_sp)(z + 3));
2643 1617489 : return gtrunc2n(gel(x,1), s);
2644 : }
2645 15666109 : gel(z,1) = gtrunc2n(gel(x,1), s);
2646 15666109 : return z;
2647 0 : default: pari_err_TYPE("gtrunc2n",x);
2648 : return NULL; /* LCOV_EXCL_LINE */
2649 : }
2650 : }
2651 :
2652 : /* e = number of error bits on integral part */
2653 : GEN
2654 1142511 : gcvtoi(GEN x, long *e)
2655 : {
2656 1142511 : long tx = typ(x), lx, e1;
2657 : GEN y;
2658 :
2659 1142511 : if (tx == t_REAL)
2660 : {
2661 1142378 : long ex = expo(x); if (ex < 0) { *e = ex; return gen_0; }
2662 1142287 : e1 = ex - bit_prec(x) + 1;
2663 1142287 : y = mantissa2nr(x, e1);
2664 1142287 : if (e1 <= 0) { pari_sp av = avma; e1 = expo(subri(x,y)); set_avma(av); }
2665 1142287 : *e = e1; return y;
2666 : }
2667 133 : *e = -(long)HIGHEXPOBIT;
2668 133 : if (is_matvec_t(tx))
2669 : {
2670 : long i;
2671 28 : y = cgetg_copy(x, &lx);
2672 84 : for (i=1; i<lx; i++)
2673 : {
2674 56 : gel(y,i) = gcvtoi(gel(x,i),&e1);
2675 56 : if (e1 > *e) *e = e1;
2676 : }
2677 28 : return y;
2678 : }
2679 105 : return gtrunc(x);
2680 : }
2681 :
2682 : int
2683 673894 : isint(GEN n, GEN *ptk)
2684 : {
2685 673894 : switch(typ(n))
2686 : {
2687 400817 : case t_INT: *ptk = n; return 1;
2688 221865 : case t_REAL: {
2689 221865 : pari_sp av0 = avma;
2690 221865 : GEN z = floorr(n);
2691 221865 : pari_sp av = avma;
2692 221865 : long s = signe(subri(n, z));
2693 221865 : if (s) return gc_bool(av0,0);
2694 189063 : *ptk = z; return gc_bool(av,1);
2695 : }
2696 33502 : case t_FRAC: return 0;
2697 17521 : case t_COMPLEX: return gequal0(gel(n,2)) && isint(gel(n,1),ptk);
2698 0 : case t_QUAD: return gequal0(gel(n,3)) && isint(gel(n,2),ptk);
2699 189 : default: pari_err_TYPE("isint",n);
2700 : return 0; /* LCOV_EXCL_LINE */
2701 : }
2702 : }
2703 :
2704 : int
2705 326463 : issmall(GEN n, long *ptk)
2706 : {
2707 326463 : pari_sp av = avma;
2708 : GEN z;
2709 : long k;
2710 326463 : if (!isint(n, &z)) return 0;
2711 324881 : k = itos_or_0(z); set_avma(av);
2712 324881 : if (k || lgefint(z) == 2) { *ptk = k; return 1; }
2713 0 : return 0;
2714 : }
2715 :
2716 : /* smallest integer greater than any incarnations of the real x
2717 : * Avoid "precision loss in truncation" */
2718 : GEN
2719 1058914 : ceil_safe(GEN x)
2720 : {
2721 1058914 : pari_sp av = avma;
2722 1058914 : long e, tx = typ(x);
2723 : GEN y;
2724 :
2725 1058914 : if (is_rational_t(tx)) return gceil(x);
2726 1058616 : y = gcvtoi(x,&e);
2727 1058616 : if (gsigne(x) >= 0)
2728 : {
2729 1058114 : if (e < 0) e = 0;
2730 1058114 : y = addii(y, int2n(e));
2731 : }
2732 1058616 : return gc_INT(av, y);
2733 : }
2734 : /* largest integer smaller than any incarnations of the real x
2735 : * Avoid "precision loss in truncation" */
2736 : GEN
2737 72059 : floor_safe(GEN x)
2738 : {
2739 72059 : pari_sp av = avma;
2740 72059 : long e, tx = typ(x);
2741 : GEN y;
2742 :
2743 72059 : if (is_rational_t(tx)) return gfloor(x);
2744 72059 : y = gcvtoi(x,&e);
2745 72059 : if (gsigne(x) <= 0)
2746 : {
2747 21 : if (e < 0) e = 0;
2748 21 : y = subii(y, int2n(e));
2749 : }
2750 72059 : return gc_INT(av, y);
2751 : }
2752 :
2753 : GEN
2754 11641 : ser2rfrac_i(GEN x)
2755 : {
2756 11641 : long e = valser(x);
2757 11641 : GEN a = ser2pol_i(x, lg(x));
2758 11641 : if (e) {
2759 6972 : if (e > 0) a = RgX_shift_shallow(a, e);
2760 0 : else a = gred_rfrac_simple(a, pol_xn(-e, varn(a)));
2761 : }
2762 11641 : return a;
2763 : }
2764 :
2765 : static GEN
2766 707 : ser2rfrac(GEN x)
2767 : {
2768 707 : pari_sp av = avma;
2769 707 : return gc_GEN(av, ser2rfrac_i(x));
2770 : }
2771 :
2772 : /* x t_PADIC, truncate to rational (t_INT/t_FRAC) */
2773 : GEN
2774 689519 : padic_to_Q(GEN x)
2775 : {
2776 689519 : GEN u = padic_u(x), p;
2777 : long v;
2778 689519 : if (!signe(u)) return gen_0;
2779 683443 : v = valp(x);
2780 683443 : if (!v) return icopy(u);
2781 258545 : p = padic_p(x);
2782 258545 : if (v > 0)
2783 : {
2784 258426 : pari_sp av = avma;
2785 258426 : return gc_INT(av, mulii(u, powiu(p,v)));
2786 : }
2787 119 : retmkfrac(icopy(u), powiu(p,-v));
2788 : }
2789 : GEN
2790 14 : padic_to_Q_shallow(GEN x)
2791 : {
2792 14 : GEN u = padic_u(x), p, q, q2;
2793 : long v;
2794 14 : if (!signe(u)) return gen_0;
2795 14 : q = padic_pd(x); q2 = shifti(q,-1);
2796 14 : v = valp(x);
2797 14 : u = Fp_center_i(u, q, q2);
2798 14 : if (!v) return u;
2799 7 : p = padic_p(x);
2800 7 : if (v > 0) return mulii(u, powiu(p,v));
2801 7 : return mkfrac(u, powiu(p,-v));
2802 : }
2803 : static GEN
2804 735 : Qp_to_Q(GEN c)
2805 : {
2806 735 : switch(typ(c))
2807 : {
2808 721 : case t_INT:
2809 721 : case t_FRAC: break;
2810 14 : case t_PADIC: c = padic_to_Q_shallow(c); break;
2811 0 : default: pari_err_TYPE("padic_to_Q", c);
2812 : }
2813 735 : return c;
2814 : }
2815 : GEN
2816 931 : QpV_to_QV(GEN x) { pari_APPLY_same(Qp_to_Q(gel(x,i))); }
2817 :
2818 : GEN
2819 597763 : gtrunc(GEN x)
2820 : {
2821 597763 : switch(typ(x))
2822 : {
2823 84175 : case t_INT: return icopy(x);
2824 49098 : case t_REAL: return truncr(x);
2825 56 : case t_FRAC: return divii(gel(x,1),gel(x,2));
2826 432710 : case t_PADIC: return padic_to_Q(x);
2827 42 : case t_POL: return RgX_copy(x);
2828 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2829 679 : case t_SER: return ser2rfrac(x);
2830 187215 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gtrunc(gel(x,i)));
2831 : }
2832 56 : pari_err_TYPE("gtrunc",x);
2833 : return NULL; /* LCOV_EXCL_LINE */
2834 : }
2835 :
2836 : GEN
2837 217 : trunc0(GEN x, GEN *pte)
2838 : {
2839 217 : if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
2840 189 : return gtrunc(x);
2841 : }
2842 : /*******************************************************************/
2843 : /* */
2844 : /* CONVERSIONS --> INT, POL & SER */
2845 : /* */
2846 : /*******************************************************************/
2847 :
2848 : /* return a_(n-1) B^(n-1) + ... + a_0, where B = 2^32.
2849 : * The a_i are 32bits integers */
2850 : GEN
2851 24804 : mkintn(long n, ...)
2852 : {
2853 : va_list ap;
2854 : GEN x, y;
2855 : long i;
2856 : #ifdef LONG_IS_64BIT
2857 21294 : long e = (n&1);
2858 21294 : n = (n+1) >> 1;
2859 : #endif
2860 24804 : va_start(ap,n);
2861 24804 : x = cgetipos(n+2);
2862 24804 : y = int_MSW(x);
2863 87714 : for (i=0; i <n; i++)
2864 : {
2865 : #ifdef LONG_IS_64BIT
2866 49140 : ulong a = (e && !i)? 0: (ulong) va_arg(ap, unsigned int);
2867 49140 : ulong b = (ulong) va_arg(ap, unsigned int);
2868 49140 : *y = (a << 32) | b;
2869 : #else
2870 13770 : *y = (ulong) va_arg(ap, unsigned int);
2871 : #endif
2872 62910 : y = int_precW(y);
2873 : }
2874 24804 : va_end(ap);
2875 24804 : return int_normalize(x, 0);
2876 : }
2877 :
2878 : /* 2^32 a + b */
2879 : GEN
2880 328944 : uu32toi(ulong a, ulong b)
2881 : {
2882 : #ifdef LONG_IS_64BIT
2883 269077 : return utoi((a<<32) | b);
2884 : #else
2885 59867 : return uutoi(a, b);
2886 : #endif
2887 : }
2888 : /* - (2^32 a + b), assume a or b != 0 */
2889 : GEN
2890 114643 : uu32toineg(ulong a, ulong b)
2891 : {
2892 : #ifdef LONG_IS_64BIT
2893 98108 : return utoineg((a<<32) | b);
2894 : #else
2895 16535 : return uutoineg(a, b);
2896 : #endif
2897 : }
2898 :
2899 : /* return a_(n-1) x^(n-1) + ... + a_0 */
2900 : GEN
2901 5629980 : mkpoln(long n, ...)
2902 : {
2903 : va_list ap;
2904 : GEN x, y;
2905 5629980 : long i, l = n+2;
2906 5629980 : va_start(ap,n);
2907 5629980 : x = cgetg(l, t_POL); y = x + 2;
2908 5629980 : x[1] = evalvarn(0);
2909 24684543 : for (i=n-1; i >= 0; i--) gel(y,i) = va_arg(ap, GEN);
2910 5629980 : va_end(ap); return normalizepol_lg(x, l);
2911 : }
2912 :
2913 : /* return [a_1, ..., a_n] */
2914 : GEN
2915 1221225 : mkvecn(long n, ...)
2916 : {
2917 : va_list ap;
2918 : GEN x;
2919 : long i;
2920 1221225 : va_start(ap,n);
2921 1221225 : x = cgetg(n+1, t_VEC);
2922 7852382 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2923 1221225 : va_end(ap); return x;
2924 : }
2925 :
2926 : GEN
2927 1379 : mkcoln(long n, ...)
2928 : {
2929 : va_list ap;
2930 : GEN x;
2931 : long i;
2932 1379 : va_start(ap,n);
2933 1379 : x = cgetg(n+1, t_COL);
2934 11032 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2935 1379 : va_end(ap); return x;
2936 : }
2937 :
2938 : GEN
2939 153621 : mkvecsmalln(long n, ...)
2940 : {
2941 : va_list ap;
2942 : GEN x;
2943 : long i;
2944 153621 : va_start(ap,n);
2945 153621 : x = cgetg(n+1, t_VECSMALL);
2946 1086230 : for (i=1; i <= n; i++) x[i] = va_arg(ap, long);
2947 153621 : va_end(ap); return x;
2948 : }
2949 :
2950 : GEN
2951 45306395 : scalarpol(GEN x, long v)
2952 : {
2953 : GEN y;
2954 45306395 : if (isrationalzero(x)) return zeropol(v);
2955 5856567 : y = cgetg(3,t_POL);
2956 5856567 : y[1] = gequal0(x)? evalvarn(v)
2957 5856567 : : evalvarn(v) | evalsigne(1);
2958 5856567 : gel(y,2) = gcopy(x); return y;
2959 : }
2960 : GEN
2961 3975112 : scalarpol_shallow(GEN x, long v)
2962 : {
2963 : GEN y;
2964 3975112 : if (isrationalzero(x)) return zeropol(v);
2965 3863509 : y = cgetg(3,t_POL);
2966 3863509 : y[1] = gequal0(x)? evalvarn(v)
2967 3863509 : : evalvarn(v) | evalsigne(1);
2968 3863509 : gel(y,2) = x; return y;
2969 : }
2970 :
2971 : /* x0 + x1*T, do not assume x1 != 0 */
2972 : GEN
2973 334445 : deg1pol(GEN x1, GEN x0,long v)
2974 : {
2975 334445 : GEN x = cgetg(4,t_POL);
2976 334445 : x[1] = evalsigne(1) | evalvarn(v);
2977 334445 : gel(x,2) = x0 == gen_0? x0: gcopy(x0); /* gen_0 frequent */
2978 334445 : gel(x,3) = gcopy(x1); return normalizepol_lg(x,4);
2979 : }
2980 :
2981 : /* same, no copy */
2982 : GEN
2983 9726617 : deg1pol_shallow(GEN x1, GEN x0,long v)
2984 : {
2985 9726617 : GEN x = cgetg(4,t_POL);
2986 9726617 : x[1] = evalsigne(1) | evalvarn(v);
2987 9726617 : gel(x,2) = x0;
2988 9726617 : gel(x,3) = x1; return normalizepol_lg(x,4);
2989 : }
2990 :
2991 : GEN
2992 324129 : deg2pol_shallow(GEN x2, GEN x1, GEN x0, long v)
2993 : {
2994 324129 : GEN x = cgetg(5,t_POL);
2995 324129 : x[1] = evalsigne(1) | evalvarn(v);
2996 324129 : gel(x,2) = x0;
2997 324129 : gel(x,3) = x1;
2998 324129 : gel(x,4) = x2;
2999 324129 : return normalizepol_lg(x,5);
3000 : }
3001 :
3002 : static GEN
3003 3479032 : _gtopoly(GEN x, long v, int reverse)
3004 : {
3005 3479032 : long tx = typ(x);
3006 : GEN y;
3007 :
3008 3479032 : if (v<0) v = 0;
3009 3479032 : switch(tx)
3010 : {
3011 21 : case t_POL:
3012 21 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3013 21 : y = RgX_copy(x); break;
3014 28 : case t_SER:
3015 28 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3016 28 : y = ser2rfrac(x);
3017 28 : if (typ(y) != t_POL)
3018 0 : pari_err_DOMAIN("gtopoly", "valuation", "<", gen_0, x);
3019 28 : break;
3020 42 : case t_RFRAC:
3021 : {
3022 42 : GEN a = gel(x,1), b = gel(x,2);
3023 42 : long vb = varn(b);
3024 42 : if (varncmp(vb, v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3025 42 : if (typ(a) != t_POL || varn(a) != vb) return zeropol(v);
3026 21 : y = RgX_div(a,b); break;
3027 : }
3028 337183 : case t_VECSMALL: x = zv_to_ZV(x); /* fall through */
3029 3478325 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
3030 : {
3031 3478325 : long j, k, lx = lg(x);
3032 : GEN z;
3033 3478325 : if (tx == t_QFB) lx--;
3034 3478325 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("gtopoly", x, "<=", v);
3035 3478325 : y = cgetg(lx+1, t_POL);
3036 3478325 : y[1] = evalvarn(v);
3037 3478325 : if (reverse) {
3038 2434005 : x--;
3039 26224373 : for (j=2; j<=lx; j++) gel(y,j) = gel(x,j);
3040 : } else {
3041 5369188 : for (j=2, k=lx; k>=2; j++) gel(y,j) = gel(x,--k);
3042 : }
3043 3478325 : z = RgX_copy(normalizepol_lg(y,lx+1));
3044 3478325 : settyp(y, t_VECSMALL);/* left on stack */
3045 3478325 : return z;
3046 : }
3047 616 : default:
3048 616 : if (is_scalar_t(tx)) return scalarpol(x,v);
3049 7 : pari_err_TYPE("gtopoly",x);
3050 : return NULL; /* LCOV_EXCL_LINE */
3051 : }
3052 70 : setvarn(y,v); return y;
3053 : }
3054 :
3055 : GEN
3056 2434040 : gtopolyrev(GEN x, long v) { return _gtopoly(x,v,1); }
3057 :
3058 : GEN
3059 1044992 : gtopoly(GEN x, long v) { return _gtopoly(x,v,0); }
3060 :
3061 : static GEN
3062 1092 : gtovecpost(GEN x, long n)
3063 : {
3064 1092 : long i, imax, lx, tx = typ(x);
3065 1092 : GEN y = zerovec(n);
3066 :
3067 1092 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,1) = gcopy(x); return y; }
3068 343 : switch(tx)
3069 : {
3070 56 : case t_POL:
3071 56 : lx = lg(x); imax = minss(lx-2, n);
3072 224 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,lx-i));
3073 56 : return y;
3074 28 : case t_SER:
3075 28 : lx = lg(x); imax = minss(lx-2, n); x++;
3076 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3077 28 : return y;
3078 28 : case t_QFB:
3079 28 : lx = lg(x)-1; /* remove discriminant */
3080 28 : imax = minss(lx-1, n);
3081 112 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3082 28 : return y;
3083 28 : case t_VEC: case t_COL:
3084 28 : lx = lg(x); imax = minss(lx-1, n);
3085 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3086 28 : return y;
3087 63 : case t_LIST:
3088 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3089 56 : x = list_data(x); lx = x? lg(x): 1;
3090 56 : imax = minss(lx-1, n);
3091 252 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3092 56 : return y;
3093 56 : case t_STR:
3094 : {
3095 56 : char *s = GSTR(x);
3096 56 : imax = minss(strlen(s), n); s--;
3097 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3098 56 : return y;
3099 : }
3100 28 : case t_VECSMALL:
3101 28 : lx = lg(x); imax = minss(lx-1, n);
3102 84 : for (i=1; i<=imax; i++) gel(y,i) = stoi(x[i]);
3103 28 : return y;
3104 56 : default: pari_err_TYPE("gtovec",x);
3105 : return NULL; /*LCOV_EXCL_LINE*/
3106 : }
3107 : }
3108 :
3109 : static GEN
3110 3773 : init_vectopre(long a, long n, GEN y, long *imax)
3111 : {
3112 3773 : if (n <= a) { *imax = n; return y; }
3113 2975 : *imax = a; return y + n - a;
3114 : }
3115 : /* assume n > 0 */
3116 : static GEN
3117 3829 : gtovecpre(GEN x, long n)
3118 : {
3119 3829 : long a, i, imax, lx, tx = typ(x);
3120 3829 : GEN y = zerovec(n), y0;
3121 :
3122 3829 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,n) = gcopy(x); return y; }
3123 3773 : switch(tx)
3124 : {
3125 56 : case t_POL:
3126 56 : lx = lg(x); a = lx-2;
3127 56 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x -= a-imax;
3128 224 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,lx-i));
3129 56 : return y;
3130 3458 : case t_SER:
3131 3458 : a = lg(x)-2; x++;
3132 3458 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3133 138796 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3134 3458 : return y;
3135 28 : case t_QFB:
3136 28 : a = lg(x)-2; /* remove discriminant */
3137 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3138 112 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3139 28 : return y;
3140 28 : case t_VEC: case t_COL:
3141 28 : a = lg(x)-1;
3142 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3143 84 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3144 28 : return y;
3145 63 : case t_LIST:
3146 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3147 56 : x = list_data(x); a = x? lg(x)-1: 0;
3148 56 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3149 252 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3150 56 : return y;
3151 56 : case t_STR:
3152 : {
3153 56 : char *s = GSTR(x);
3154 56 : a = strlen(s);
3155 56 : y0 = init_vectopre(a, n, y, &imax); s--; if (imax == n) s += a-imax;
3156 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3157 56 : return y;
3158 : }
3159 28 : case t_VECSMALL:
3160 28 : a = lg(x)-1;
3161 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3162 84 : for (i=1; i<=imax; i++) gel(y0,i) = stoi(x[i]);
3163 28 : return y;
3164 56 : default: pari_err_TYPE("gtovec",x);
3165 : return NULL; /*LCOV_EXCL_LINE*/
3166 : }
3167 : }
3168 : GEN
3169 123878 : gtovec0(GEN x, long n)
3170 : {
3171 123878 : if (!n) return gtovec(x);
3172 4921 : if (n > 0) return gtovecpost(x, n);
3173 3829 : return gtovecpre(x, -n);
3174 : }
3175 :
3176 : GEN
3177 119447 : gtovec(GEN x)
3178 : {
3179 119447 : long i, lx, tx = typ(x);
3180 : GEN y;
3181 :
3182 119447 : if (is_scalar_t(tx)) return mkveccopy(x);
3183 119314 : switch(tx)
3184 : {
3185 15477 : case t_POL:
3186 15477 : lx=lg(x); y=cgetg(lx-1,t_VEC);
3187 1501535 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,lx-i));
3188 15477 : return y;
3189 385 : case t_SER:
3190 385 : lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
3191 12264 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
3192 385 : return y;
3193 28 : case t_RFRAC: return mkveccopy(x);
3194 70049 : case t_QFB:
3195 70049 : retmkvec3(icopy(gel(x,1)), icopy(gel(x,2)), icopy(gel(x,3)));
3196 31283 : case t_VEC: case t_COL: case t_MAT:
3197 31283 : lx=lg(x); y=cgetg(lx,t_VEC);
3198 1665027 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3199 31283 : return y;
3200 426 : case t_LIST:
3201 426 : if (list_typ(x) == t_LIST_MAP) return mapdomain(x);
3202 412 : x = list_data(x); lx = x? lg(x): 1;
3203 412 : y = cgetg(lx, t_VEC);
3204 20373 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3205 412 : return y;
3206 105 : case t_STR:
3207 : {
3208 105 : char *s = GSTR(x);
3209 105 : lx = strlen(s)+1; y = cgetg(lx, t_VEC);
3210 105 : s--;
3211 340239 : for (i=1; i<lx; i++) gel(y,i) = chartoGENstr(s[i]);
3212 105 : return y;
3213 : }
3214 1498 : case t_VECSMALL:
3215 1498 : return vecsmall_to_vec(x);
3216 63 : case t_ERROR:
3217 63 : lx=lg(x); y = cgetg(lx,t_VEC);
3218 63 : gel(y,1) = errname(x);
3219 168 : for (i=2; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3220 63 : return y;
3221 0 : default: pari_err_TYPE("gtovec",x);
3222 : return NULL; /*LCOV_EXCL_LINE*/
3223 : }
3224 : }
3225 :
3226 : GEN
3227 308 : gtovecrev0(GEN x, long n)
3228 : {
3229 308 : GEN y = gtovec0(x, -n);
3230 280 : vecreverse_inplace(y);
3231 280 : return y;
3232 : }
3233 : GEN
3234 0 : gtovecrev(GEN x) { return gtovecrev0(x, 0); }
3235 :
3236 : GEN
3237 4004 : gtocol0(GEN x, long n)
3238 : {
3239 : GEN y;
3240 4004 : if (!n) return gtocol(x);
3241 3654 : y = gtovec0(x, n);
3242 3598 : settyp(y, t_COL); return y;
3243 : }
3244 : GEN
3245 350 : gtocol(GEN x)
3246 : {
3247 : long lx, tx, i, j, h;
3248 : GEN y;
3249 350 : tx = typ(x);
3250 350 : if (tx != t_MAT) { y = gtovec(x); settyp(y, t_COL); return y; }
3251 14 : lx = lg(x); if (lx == 1) return cgetg(1, t_COL);
3252 14 : h = lgcols(x); y = cgetg(h, t_COL);
3253 42 : for (i = 1 ; i < h; i++) {
3254 28 : gel(y,i) = cgetg(lx, t_VEC);
3255 112 : for (j = 1; j < lx; j++) gmael(y,i,j) = gcopy(gcoeff(x,i,j));
3256 : }
3257 14 : return y;
3258 : }
3259 :
3260 : GEN
3261 294 : gtocolrev0(GEN x, long n)
3262 : {
3263 294 : GEN y = gtocol0(x, -n);
3264 266 : long ly = lg(y), lim = ly>>1, i;
3265 763 : for (i = 1; i <= lim; i++) swap(gel(y,i), gel(y,ly-i));
3266 266 : return y;
3267 : }
3268 : GEN
3269 0 : gtocolrev(GEN x) { return gtocolrev0(x, 0); }
3270 :
3271 : static long
3272 19145 : Itos(GEN x)
3273 : {
3274 19145 : if (typ(x) != t_INT) pari_err_TYPE("vectosmall",x);
3275 19145 : return itos(x);
3276 : }
3277 :
3278 : static GEN
3279 98 : gtovecsmallpost(GEN x, long n)
3280 : {
3281 : long i, imax, lx;
3282 98 : GEN y = zero_Flv(n);
3283 :
3284 98 : switch(typ(x))
3285 : {
3286 7 : case t_INT:
3287 7 : y[1] = itos(x); return y;
3288 14 : case t_POL:
3289 14 : lx=lg(x); imax = minss(lx-2, n);
3290 56 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,lx-i));
3291 14 : return y;
3292 7 : case t_SER:
3293 7 : lx=lg(x); imax = minss(lx-2, n); x++;
3294 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3295 7 : return y;
3296 7 : case t_VEC: case t_COL:
3297 7 : lx=lg(x); imax = minss(lx-1, n);
3298 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3299 7 : return y;
3300 14 : case t_LIST:
3301 14 : x = list_data(x); lx = x? lg(x): 1;
3302 14 : imax = minss(lx-1, n);
3303 63 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3304 14 : return y;
3305 7 : case t_VECSMALL:
3306 7 : lx=lg(x);
3307 7 : imax = minss(lx-1, n);
3308 21 : for (i=1; i<=imax; i++) y[i] = x[i];
3309 7 : return y;
3310 14 : case t_STR:
3311 : {
3312 14 : unsigned char *s = (unsigned char*)GSTR(x);
3313 14 : imax = minss(strlen((const char *)s), n); s--;
3314 56 : for (i=1; i<=imax; i++) y[i] = (long)s[i];
3315 14 : return y;
3316 : }
3317 28 : default: pari_err_TYPE("gtovecsmall",x);
3318 : return NULL; /*LCOV_EXCL_LINE*/
3319 : }
3320 : }
3321 : static GEN
3322 98 : gtovecsmallpre(GEN x, long n)
3323 : {
3324 98 : GEN y = zero_Flv(n), y0;
3325 : long a, i, imax, lx;
3326 :
3327 98 : switch(typ(x))
3328 : {
3329 7 : case t_INT:
3330 7 : y[n] = itos(x); return y;
3331 14 : case t_POL:
3332 14 : lx = lg(x); a = lx-2;
3333 14 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x -= a-imax;
3334 56 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,lx-i));
3335 14 : return y;
3336 7 : case t_SER:
3337 7 : a = lg(x)-2; x++;
3338 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3339 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3340 7 : return y;
3341 7 : case t_VEC: case t_COL:
3342 7 : a = lg(x)-1;
3343 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3344 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3345 7 : return y;
3346 14 : case t_LIST:
3347 14 : x = list_data(x); a = x? lg(x)-1: 0;
3348 14 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3349 63 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3350 14 : return y;
3351 7 : case t_VECSMALL:
3352 7 : a = lg(x)-1;
3353 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3354 21 : for (i=1; i<=imax; i++) y0[i] = x[i];
3355 7 : return y;
3356 14 : case t_STR:
3357 : {
3358 14 : unsigned char *s = (unsigned char*)GSTR(x);
3359 14 : a = strlen((const char *)s);
3360 14 : y0 = init_vectopre(a, n, y, &imax); s--; if (imax == n) s += a-imax;
3361 56 : for (i=1; i<=imax; i++) y0[i] = (long)s[i];
3362 14 : return y;
3363 : }
3364 28 : default: pari_err_TYPE("gtovecsmall",x);
3365 : return NULL; /*LCOV_EXCL_LINE*/
3366 : }
3367 : }
3368 :
3369 : GEN
3370 7819 : gtovecsmall0(GEN x, long n)
3371 : {
3372 7819 : if (!n) return gtovecsmall(x);
3373 196 : if (n > 0) return gtovecsmallpost(x, n);
3374 98 : return gtovecsmallpre(x, -n);
3375 : }
3376 :
3377 : GEN
3378 19187 : gtovecsmall(GEN x)
3379 : {
3380 : GEN V;
3381 : long l, i;
3382 :
3383 19187 : switch(typ(x))
3384 : {
3385 126 : case t_INT: return mkvecsmall(itos(x));
3386 140 : case t_STR:
3387 : {
3388 140 : unsigned char *s = (unsigned char*)GSTR(x);
3389 140 : l = strlen((const char *)s);
3390 140 : V = cgetg(l+1, t_VECSMALL);
3391 140 : s--;
3392 2541 : for (i=1; i<=l; i++) V[i] = (long)s[i];
3393 140 : return V;
3394 : }
3395 13139 : case t_VECSMALL: return leafcopy(x);
3396 14 : case t_LIST:
3397 14 : x = list_data(x);
3398 14 : if (!x) return cgetg(1, t_VECSMALL);
3399 : /* fall through */
3400 : case t_VEC: case t_COL:
3401 5733 : l = lg(x); V = cgetg(l,t_VECSMALL);
3402 24577 : for(i=1; i<l; i++) V[i] = Itos(gel(x,i));
3403 5733 : return V;
3404 14 : case t_POL:
3405 14 : l = lg(x); V = cgetg(l-1,t_VECSMALL);
3406 63 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,l-i));
3407 14 : return V;
3408 7 : case t_SER:
3409 7 : l = lg(x); V = cgetg(l-1,t_VECSMALL); x++;
3410 21 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,i));
3411 7 : return V;
3412 28 : default:
3413 28 : pari_err_TYPE("vectosmall",x);
3414 : return NULL; /* LCOV_EXCL_LINE */
3415 : }
3416 : }
3417 :
3418 : GEN
3419 327 : compo(GEN x, long n)
3420 : {
3421 327 : long tx = typ(x);
3422 327 : ulong l, lx = (ulong)lg(x);
3423 :
3424 327 : if (!is_recursive_t(tx))
3425 : {
3426 63 : if (tx == t_VECSMALL)
3427 : {
3428 21 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3429 21 : if ((ulong)n >= lx) pari_err_COMPONENT("", ">", utoi(lx-1), stoi(n));
3430 7 : return stoi(x[n]);
3431 : }
3432 42 : pari_err_TYPE("component [leaf]", x);
3433 : }
3434 264 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3435 257 : if (tx == t_LIST) {
3436 28 : x = list_data(x); tx = t_VEC;
3437 28 : lx = (ulong)(x? lg(x): 1);
3438 : }
3439 257 : l = (ulong)lontyp[tx] + (ulong)n-1; /* beware overflow */
3440 257 : if (l >= lx) pari_err_COMPONENT("", ">", utoi(lx-lontyp[tx]), stoi(n));
3441 173 : return gcopy(gel(x,l));
3442 : }
3443 :
3444 : /* assume x a t_POL */
3445 : static GEN
3446 2601098 : _polcoef(GEN x, long n, long v)
3447 : {
3448 2601098 : long i, w, lx = lg(x), dx = lx-3;
3449 : GEN z;
3450 2601098 : if (dx < 0) return gen_0;
3451 2023122 : if (v < 0 || v == (w=varn(x)))
3452 1340412 : return (n < 0 || n > dx)? gen_0: gel(x,n+2);
3453 682710 : if (varncmp(w,v) > 0) return n? gen_0: x;
3454 : /* w < v */
3455 681870 : z = cgetg(lx, t_POL); z[1] = x[1];
3456 2730175 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3457 681870 : z = normalizepol_lg(z, lx);
3458 681870 : switch(lg(z))
3459 : {
3460 28295 : case 2: z = gen_0; break;
3461 412555 : case 3: z = gel(z,2); break;
3462 : }
3463 681870 : return z;
3464 : }
3465 :
3466 : /* assume x a t_SER */
3467 : static GEN
3468 111958 : _sercoef(GEN x, long n, long v)
3469 : {
3470 111958 : long i, w = varn(x), lx = lg(x), dx = lx-3, N;
3471 : GEN z;
3472 111958 : if (v < 0) v = w;
3473 111958 : N = v == w? n - valser(x): n;
3474 111958 : if (dx < 0)
3475 : {
3476 21 : if (N >= 0) pari_err_DOMAIN("polcoef", "t_SER", "=", x, x);
3477 14 : return gen_0;
3478 : }
3479 111937 : if (v == w)
3480 : {
3481 111895 : if (!dx && !signe(x) && !isinexact(gel(x,2))) dx = -1;
3482 111895 : if (N > dx)
3483 28 : pari_err_DOMAIN("polcoef", "degree", ">", stoi(dx+valser(x)), stoi(n));
3484 111867 : return (N < 0)? gen_0: gel(x,N+2);
3485 : }
3486 42 : if (varncmp(w,v) > 0) return N? gen_0: x;
3487 : /* w < v */
3488 28 : z = cgetg(lx, t_SER); z[1] = x[1];
3489 91 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3490 28 : return normalizeser(z);
3491 : }
3492 :
3493 : /* assume x a t_RFRAC(n) */
3494 : static GEN
3495 21 : _rfraccoef(GEN x, long n, long v)
3496 : {
3497 21 : GEN p = gel(x,1), q = gel(x,2), z;
3498 21 : long vp = gvar(p), vq = gvar(q), vz;
3499 21 : if (v < 0) v = varncmp(vp, vq) < 0? vp: vq;
3500 21 : vz = v;
3501 21 : if (varncmp(vp, v) < 0 || varncmp(vq, v) < 0) vz = fetch_var_higher();
3502 21 : if (vp != v) p = swap_vars(p, v, vz);
3503 21 : if (vq != v) q = swap_vars(q, v, vz);
3504 21 : n += degpol(q);
3505 21 : z = gdiv(_polcoef(p, n, vz), leading_coeff(q));
3506 21 : if (vz != v) (void)delete_var();
3507 21 : if (!RgX_is_monomial(q)) pari_err_TYPE("polcoef", x);
3508 21 : return z;
3509 : }
3510 :
3511 : GEN
3512 3613484 : polcoef_i(GEN x, long n, long v)
3513 : {
3514 3613484 : long tx = typ(x);
3515 3613484 : switch(tx)
3516 : {
3517 2601077 : case t_POL: return _polcoef(x,n,v);
3518 111958 : case t_SER: return _sercoef(x,n,v);
3519 21 : case t_RFRAC: return _rfraccoef(x,n,v);
3520 : }
3521 900428 : if (!is_scalar_t(tx)) pari_err_TYPE("polcoef", x);
3522 900232 : return n? gen_0: x;
3523 : }
3524 :
3525 : /* with respect to the main variable if v<0, with respect to the variable v
3526 : * otherwise. v ignored if x is not a polynomial/series. */
3527 : GEN
3528 727580 : polcoef(GEN x, long n, long v)
3529 : {
3530 727580 : pari_sp av = avma;
3531 727580 : x = polcoef_i(x,n,v);
3532 727349 : if (x == gen_0) return x;
3533 130795 : return (avma == av)? gcopy(x): gc_GEN(av, x);
3534 : }
3535 :
3536 : static GEN
3537 242315 : vecdenom(GEN v, long imin, long imax)
3538 : {
3539 242315 : long i = imin;
3540 : GEN s;
3541 242315 : if (imin > imax) return gen_1;
3542 242315 : s = denom_i(gel(v,i));
3543 2104290 : for (i++; i<=imax; i++)
3544 : {
3545 1861975 : GEN t = denom_i(gel(v,i));
3546 1861975 : if (t != gen_1) s = glcm(s,t);
3547 : }
3548 242315 : return s;
3549 : }
3550 : static GEN denompol(GEN x, long v);
3551 : static GEN
3552 14 : vecdenompol(GEN v, long imin, long imax, long vx)
3553 : {
3554 14 : long i = imin;
3555 : GEN s;
3556 14 : if (imin > imax) return gen_1;
3557 14 : s = denompol(gel(v,i), vx);
3558 14 : for (i++; i<=imax; i++)
3559 : {
3560 0 : GEN t = denompol(gel(v,i), vx);
3561 0 : if (t != gen_1) s = glcm(s,t);
3562 : }
3563 14 : return s;
3564 : }
3565 : GEN
3566 12319313 : denom_i(GEN x)
3567 : {
3568 12319313 : switch(typ(x))
3569 : {
3570 4730715 : case t_INT:
3571 : case t_REAL:
3572 : case t_INTMOD:
3573 : case t_FFELT:
3574 : case t_PADIC:
3575 : case t_SER:
3576 4730715 : case t_VECSMALL: return gen_1;
3577 81713 : case t_FRAC: return gel(x,2);
3578 294 : case t_COMPLEX: return vecdenom(x,1,2);
3579 69069 : case t_QUAD: return vecdenom(x,2,3);
3580 42 : case t_POLMOD: return denom_i(gel(x,2));
3581 7263569 : case t_RFRAC: return gel(x,2);
3582 945 : case t_POL: return pol_1(varn(x));
3583 172952 : case t_VEC: case t_COL: case t_MAT: return vecdenom(x, 1, lg(x)-1);
3584 : }
3585 14 : pari_err_TYPE("denom",x);
3586 : return NULL; /* LCOV_EXCL_LINE */
3587 : }
3588 : /* v has lower (or equal) priority as x's main variable */
3589 : static GEN
3590 119 : denompol(GEN x, long v)
3591 : {
3592 119 : long vx, tx = typ(x);
3593 119 : if (is_scalar_t(tx)) return gen_1;
3594 105 : switch(typ(x))
3595 : {
3596 14 : case t_SER:
3597 14 : if (varn(x) != v) return x;
3598 14 : vx = valser(x); return vx < 0? pol_xn(-vx, v): pol_1(v);
3599 63 : case t_RFRAC: x = gel(x,2); return varn(x) == v? x: pol_1(v);
3600 14 : case t_POL: return pol_1(v);
3601 14 : case t_VEC: case t_COL: case t_MAT: return vecdenompol(x, 1, lg(x)-1, v);
3602 : }
3603 0 : pari_err_TYPE("denom",x);
3604 : return NULL; /* LCOV_EXCL_LINE */
3605 : }
3606 : GEN
3607 228909 : denom(GEN x) { pari_sp av = avma; return gc_GEN(av, denom_i(x)); }
3608 :
3609 : static GEN
3610 287 : denominator_v(GEN x, long v)
3611 : {
3612 287 : long v0 = gvar(x);
3613 : GEN d;
3614 287 : if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
3615 105 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3616 105 : d = denompol(x, v0);
3617 105 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3618 105 : return d;
3619 : }
3620 : GEN
3621 128149 : denominator(GEN x, GEN D)
3622 : {
3623 128149 : pari_sp av = avma;
3624 : GEN d;
3625 128149 : if (!D) return denom(x);
3626 280 : if (isint1(D))
3627 : {
3628 140 : d = Q_denom_safe(x);
3629 140 : if (!d) { set_avma(av); return gen_1; }
3630 105 : return gc_GEN(av, d);
3631 : }
3632 140 : if (!gequalX(D)) pari_err_TYPE("denominator", D);
3633 140 : return gc_upto(av, denominator_v(x, varn(D)));
3634 : }
3635 : GEN
3636 8925 : numerator(GEN x, GEN D)
3637 : {
3638 8925 : pari_sp av = avma;
3639 : long v;
3640 8925 : if (!D) return numer(x);
3641 294 : if (isint1(D)) return Q_remove_denom(x,NULL);
3642 154 : if (!gequalX(D)) pari_err_TYPE("numerator", D);
3643 154 : v = varn(D); /* optimization */
3644 154 : if (typ(x) == t_RFRAC && varn(gel(x,2)) == v) return gcopy(gel(x,1));
3645 147 : return gc_upto(av, gmul(x, denominator_v(x,v)));
3646 : }
3647 : GEN
3648 131005 : content0(GEN x, GEN D)
3649 : {
3650 131005 : pari_sp av = avma;
3651 : long v, v0;
3652 : GEN d;
3653 131005 : if (!D) return content(x);
3654 294 : if (isint1(D))
3655 : {
3656 140 : d = Q_content_safe(x);
3657 140 : return d? d: gen_1;
3658 : }
3659 154 : if (!gequalX(D)) pari_err_TYPE("content", D);
3660 154 : v = varn(D);
3661 154 : v0 = gvar(x); if (v0 == NO_VARIABLE) return gen_1;
3662 56 : if (varncmp(v0,v) > 0)
3663 : {
3664 0 : v0 = gvar2(x);
3665 0 : return v0 == NO_VARIABLE? gen_1: pol_1(v0);
3666 : }
3667 56 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3668 56 : d = content(x);
3669 : /* gsubst is needed because of content([x]) = x */
3670 56 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3671 56 : return gc_upto(av, d);
3672 : }
3673 :
3674 : GEN
3675 9050681 : numer_i(GEN x)
3676 : {
3677 9050681 : switch(typ(x))
3678 : {
3679 1789310 : case t_INT:
3680 : case t_REAL:
3681 : case t_INTMOD:
3682 : case t_FFELT:
3683 : case t_PADIC:
3684 : case t_SER:
3685 : case t_VECSMALL:
3686 1789310 : case t_POL: return x;
3687 28 : case t_POLMOD: return mkpolmod(numer_i(gel(x,2)), gel(x,1));
3688 7261154 : case t_FRAC:
3689 7261154 : case t_RFRAC: return gel(x,1);
3690 175 : case t_COMPLEX:
3691 : case t_QUAD:
3692 : case t_VEC:
3693 : case t_COL:
3694 175 : case t_MAT: return gmul(denom_i(x),x);
3695 : }
3696 14 : pari_err_TYPE("numer",x);
3697 : return NULL; /* LCOV_EXCL_LINE */
3698 : }
3699 : GEN
3700 8631 : numer(GEN x) { pari_sp av = avma; return gc_GEN(av, numer_i(x)); }
3701 :
3702 : /* Lift only intmods if v does not occur in x, lift with respect to main
3703 : * variable of x if v < 0, with respect to variable v otherwise */
3704 : GEN
3705 2523007 : lift0(GEN x, long v)
3706 : {
3707 : GEN y;
3708 :
3709 2523007 : switch(typ(x))
3710 : {
3711 30569 : case t_INT: return icopy(x);
3712 2368353 : case t_INTMOD: return v < 0? icopy(gel(x,2)): gcopy(x);
3713 93642 : case t_POLMOD:
3714 93642 : if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
3715 14 : y = cgetg(3, t_POLMOD);
3716 14 : gel(y,1) = lift0(gel(x,1),v);
3717 14 : gel(y,2) = lift0(gel(x,2),v); return y;
3718 665 : case t_PADIC: return v < 0? padic_to_Q(x): gcopy(x);
3719 8813 : case t_POL:
3720 41916 : pari_APPLY_pol(lift0(gel(x,i), v));
3721 56 : case t_SER:
3722 56 : if (ser_isexactzero(x))
3723 : {
3724 14 : if (lg(x) == 2) return gcopy(x);
3725 14 : y = scalarser(lift0(gel(x,2),v), varn(x), 1);
3726 14 : setvalser(y, valser(x)); return y;
3727 : }
3728 434 : pari_APPLY_ser(lift0(gel(x,i), v));
3729 20720 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3730 : case t_VEC: case t_COL: case t_MAT:
3731 221970 : pari_APPLY_same(lift0(gel(x,i), v));
3732 189 : default: return gcopy(x);
3733 : }
3734 : }
3735 : /* same as lift, shallow */
3736 : GEN
3737 652418 : lift_shallow(GEN x)
3738 : {
3739 : GEN y;
3740 652418 : switch(typ(x))
3741 : {
3742 212379 : case t_INTMOD: case t_POLMOD: return gel(x,2);
3743 476 : case t_PADIC: return padic_to_Q(x);
3744 0 : case t_SER:
3745 0 : if (ser_isexactzero(x))
3746 : {
3747 0 : if (lg(x) == 2) return x;
3748 0 : y = scalarser(lift_shallow(gel(x,2)), varn(x), 1);
3749 0 : setvalser(y, valser(x)); return y;
3750 : }
3751 0 : pari_APPLY_ser(lift_shallow(gel(x,i)));
3752 56336 : case t_POL:
3753 316992 : pari_APPLY_pol(lift_shallow(gel(x,i)));
3754 11263 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3755 : case t_VEC: case t_COL: case t_MAT:
3756 275905 : pari_APPLY_same(lift_shallow(gel(x,i)));
3757 371964 : default: return x;
3758 : }
3759 : }
3760 : GEN
3761 2205748 : lift(GEN x) { return lift0(x,-1); }
3762 :
3763 : GEN
3764 2004499 : liftall_shallow(GEN x)
3765 : {
3766 : GEN y;
3767 2004499 : switch(typ(x))
3768 : {
3769 534100 : case t_INTMOD: return gel(x,2);
3770 547624 : case t_POLMOD:
3771 547624 : return liftall_shallow(gel(x,2));
3772 581 : case t_PADIC: return padic_to_Q(x);
3773 556073 : case t_POL:
3774 1357398 : pari_APPLY_pol(liftall_shallow(gel(x,i)));
3775 7 : case t_SER:
3776 7 : if (ser_isexactzero(x))
3777 : {
3778 0 : if (lg(x) == 2) return x;
3779 0 : y = scalarser(liftall_shallow(gel(x,2)), varn(x), 1);
3780 0 : setvalser(y, valser(x)); return y;
3781 : }
3782 35 : pari_APPLY_ser(liftall_shallow(gel(x,i)));
3783 132811 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3784 : case t_VEC: case t_COL: case t_MAT:
3785 760725 : pari_APPLY_same(liftall_shallow(gel(x,i)));
3786 233303 : default: return x;
3787 : }
3788 : }
3789 : GEN
3790 26264 : liftall(GEN x)
3791 26264 : { pari_sp av = avma; return gc_GEN(av, liftall_shallow(x)); }
3792 :
3793 : GEN
3794 546 : liftint_shallow(GEN x)
3795 : {
3796 : GEN y;
3797 546 : switch(typ(x))
3798 : {
3799 266 : case t_INTMOD: return gel(x,2);
3800 28 : case t_PADIC: return padic_to_Q(x);
3801 21 : case t_POL:
3802 70 : pari_APPLY_pol(liftint_shallow(gel(x,i)));
3803 14 : case t_SER:
3804 14 : if (ser_isexactzero(x))
3805 : {
3806 7 : if (lg(x) == 2) return x;
3807 7 : y = scalarser(liftint_shallow(gel(x,2)), varn(x), 1);
3808 7 : setvalser(y, valser(x)); return y;
3809 : }
3810 35 : pari_APPLY_ser(liftint_shallow(gel(x,i)));
3811 161 : case t_POLMOD: case t_COMPLEX: case t_QUAD: case t_RFRAC:
3812 : case t_VEC: case t_COL: case t_MAT:
3813 504 : pari_APPLY_same(liftint_shallow(gel(x,i)));
3814 56 : default: return x;
3815 : }
3816 : }
3817 : GEN
3818 119 : liftint(GEN x)
3819 119 : { pari_sp av = avma; return gc_GEN(av, liftint_shallow(x)); }
3820 :
3821 : GEN
3822 22577944 : liftpol_shallow(GEN x)
3823 : {
3824 : GEN y;
3825 22577944 : switch(typ(x))
3826 : {
3827 2357476 : case t_POLMOD:
3828 2357476 : return liftpol_shallow(gel(x,2));
3829 3179254 : case t_POL:
3830 12989643 : pari_APPLY_pol(liftpol_shallow(gel(x,i)));
3831 7 : case t_SER:
3832 7 : if (ser_isexactzero(x))
3833 : {
3834 0 : if (lg(x) == 2) return x;
3835 0 : y = scalarser(liftpol(gel(x,2)), varn(x), 1);
3836 0 : setvalser(y, valser(x)); return y;
3837 : }
3838 35 : pari_APPLY_ser(liftpol_shallow(gel(x,i)));
3839 145467 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
3840 1020411 : pari_APPLY_same(liftpol_shallow(gel(x,i)));
3841 16895740 : default: return x;
3842 : }
3843 : }
3844 : GEN
3845 5726 : liftpol(GEN x)
3846 5726 : { pari_sp av = avma; return gc_GEN(av, liftpol_shallow(x)); }
3847 :
3848 : static GEN
3849 42518 : centerliftii(GEN x, GEN y)
3850 : {
3851 42518 : pari_sp av = avma;
3852 42518 : long i = cmpii(shifti(x,1), y);
3853 42518 : set_avma(av); return (i > 0)? subii(x,y): icopy(x);
3854 : }
3855 :
3856 : /* see lift0 */
3857 : GEN
3858 707 : centerlift0(GEN x, long v)
3859 707 : { return v < 0? centerlift(x): lift0(x,v); }
3860 : GEN
3861 60473 : centerlift(GEN x)
3862 : {
3863 : GEN u, p, pd;
3864 : long v;
3865 60473 : switch(typ(x))
3866 : {
3867 784 : case t_INT: return icopy(x);
3868 784 : case t_INTMOD: return centerliftii(gel(x,2), gel(x,1));
3869 7 : case t_POLMOD: return gcopy(gel(x,2));
3870 1533 : case t_POL:
3871 9891 : pari_APPLY_pol(centerlift(gel(x,i)));
3872 7 : case t_SER:
3873 7 : if (ser_isexactzero(x)) return lift(x);
3874 35 : pari_APPLY_ser(centerlift(gel(x,i)));
3875 5551 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3876 : case t_VEC: case t_COL: case t_MAT:
3877 56210 : pari_APPLY_same(centerlift(gel(x,i)));
3878 51800 : case t_PADIC:
3879 51800 : if (!signe(padic_u(x))) return gen_0;
3880 41734 : v = valp(x); u = padic_u(x); p = padic_p(x); pd = padic_pd(x);
3881 41734 : if (v >= 0)
3882 : { /* here p^v is an integer */
3883 41727 : pari_sp av = avma;
3884 41727 : GEN z = centerliftii(u, pd);
3885 41727 : return v? gc_INT(av, mulii(powiu(p,v), z)): z;
3886 : }
3887 7 : retmkfrac(centerliftii(u, pd), powiu(p,-v));
3888 7 : default: return gcopy(x);
3889 : }
3890 : }
3891 :
3892 : /*******************************************************************/
3893 : /* */
3894 : /* REAL & IMAGINARY PARTS */
3895 : /* */
3896 : /*******************************************************************/
3897 :
3898 : static GEN
3899 204306921 : op_ReIm(GEN f(GEN), GEN x)
3900 : {
3901 204306921 : switch(typ(x))
3902 : {
3903 595427990 : case t_POL: pari_APPLY_pol(f(gel(x,i)));
3904 64463 : case t_SER: pari_APPLY_ser(f(gel(x,i)));
3905 42 : case t_RFRAC: {
3906 42 : pari_sp av = avma;
3907 42 : GEN n, d, dxb = conj_i(gel(x,2));
3908 42 : n = gmul(gel(x,1), dxb);
3909 42 : d = gmul(gel(x,2), dxb);
3910 42 : return gc_upto(av, gdiv(f(n), d));
3911 : }
3912 20609597 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(f(gel(x, i)));
3913 : }
3914 14 : pari_err_TYPE("greal/gimag",x);
3915 : return NULL; /* LCOV_EXCL_LINE */
3916 : }
3917 :
3918 : GEN
3919 328231247 : real_i(GEN x)
3920 : {
3921 328231247 : switch(typ(x))
3922 : {
3923 174618991 : case t_INT: case t_REAL: case t_FRAC:
3924 174618991 : return x;
3925 48704792 : case t_COMPLEX:
3926 48704792 : return gel(x,1);
3927 0 : case t_QUAD:
3928 0 : return gel(x,2);
3929 : }
3930 104907464 : return op_ReIm(real_i,x);
3931 : }
3932 : GEN
3933 301930374 : imag_i(GEN x)
3934 : {
3935 301930374 : switch(typ(x))
3936 : {
3937 167817560 : case t_INT: case t_REAL: case t_FRAC:
3938 167817560 : return gen_0;
3939 34714288 : case t_COMPLEX:
3940 34714288 : return gel(x,2);
3941 0 : case t_QUAD:
3942 0 : return gel(x,3);
3943 : }
3944 99398526 : return op_ReIm(imag_i,x);
3945 : }
3946 : GEN
3947 7000 : greal(GEN x)
3948 : {
3949 7000 : switch(typ(x))
3950 : {
3951 665 : case t_INT: case t_REAL:
3952 665 : return mpcopy(x);
3953 :
3954 7 : case t_FRAC:
3955 7 : return gcopy(x);
3956 :
3957 6062 : case t_COMPLEX:
3958 6062 : return gcopy(gel(x,1));
3959 :
3960 7 : case t_QUAD:
3961 7 : return gcopy(gel(x,2));
3962 : }
3963 259 : return op_ReIm(greal,x);
3964 : }
3965 : GEN
3966 29477 : gimag(GEN x)
3967 : {
3968 29477 : switch(typ(x))
3969 : {
3970 525 : case t_INT: case t_REAL: case t_FRAC:
3971 525 : return gen_0;
3972 :
3973 28273 : case t_COMPLEX:
3974 28273 : return gcopy(gel(x,2));
3975 :
3976 7 : case t_QUAD:
3977 7 : return gcopy(gel(x,3));
3978 : }
3979 672 : return op_ReIm(gimag,x);
3980 : }
3981 :
3982 : /* return Im(x * y), x and y scalars */
3983 : GEN
3984 103649 : mulimag(GEN x, GEN y)
3985 : {
3986 103649 : if (typ(x) == t_COMPLEX)
3987 : {
3988 79282 : if (typ(y) == t_COMPLEX)
3989 : {
3990 56042 : pari_sp av = avma;
3991 56042 : GEN a = gmul(gel(x,1), gel(y,2));
3992 56042 : GEN b = gmul(gel(x,2), gel(y,1));
3993 56042 : return gc_upto(av, gadd(a, b));
3994 : }
3995 23240 : x = gel(x,2);
3996 : }
3997 24367 : else if (typ(y) == t_COMPLEX) y = gel(y,2);
3998 8750 : else return gen_0;
3999 38857 : return gmul(x,y);
4000 : }
4001 :
4002 : /* return Re(x * y), x and y scalars */
4003 : GEN
4004 29809725 : mulreal(GEN x, GEN y)
4005 : {
4006 29809725 : if (typ(x) == t_COMPLEX)
4007 : {
4008 29658475 : if (typ(y) == t_COMPLEX)
4009 : {
4010 28419744 : pari_sp av = avma;
4011 28419744 : GEN a = gmul(gel(x,1), gel(y,1));
4012 28419744 : GEN b = gmul(gel(x,2), gel(y,2));
4013 28419744 : return gc_upto(av, gsub(a, b));
4014 : }
4015 1238731 : x = gel(x,1);
4016 : }
4017 : else
4018 151250 : if (typ(y) == t_COMPLEX) y = gel(y,1);
4019 1389981 : return gmul(x,y);
4020 : }
4021 : /* Compute Re(x * y), x and y matrices of compatible dimensions
4022 : * assume scalar entries */
4023 : GEN
4024 0 : RgM_mulreal(GEN x, GEN y)
4025 : {
4026 0 : long i, j, k, l, lx = lg(x), ly = lg(y);
4027 0 : GEN z = cgetg(ly,t_MAT);
4028 0 : l = (lx == 1)? 1: lgcols(x);
4029 0 : for (j=1; j<ly; j++)
4030 : {
4031 0 : GEN zj = cgetg(l,t_COL), yj = gel(y,j);
4032 0 : gel(z,j) = zj;
4033 0 : for (i=1; i<l; i++)
4034 : {
4035 0 : pari_sp av = avma;
4036 0 : GEN c = mulreal(gcoeff(x,i,1),gel(yj,1));
4037 0 : for (k=2; k<lx; k++) c = gadd(c, mulreal(gcoeff(x,i,k),gel(yj,k)));
4038 0 : gel(zj, i) = gc_upto(av, c);
4039 : }
4040 : }
4041 0 : return z;
4042 : }
4043 :
4044 : /* Compute Re(x * y), symmetric result, x and y vectors of compatible
4045 : * dimensions; assume scalar entries */
4046 : GEN
4047 21630 : RgC_RgV_mulrealsym(GEN x, GEN y)
4048 : {
4049 21630 : long i, j, l = lg(x);
4050 21630 : GEN q = cgetg(l, t_MAT);
4051 86520 : for (j = 1; j < l; j++)
4052 : {
4053 64890 : gel(q,j) = cgetg(l,t_COL);
4054 194670 : for (i = 1; i <= j; i++)
4055 129780 : gcoeff(q,i,j) = gcoeff(q,j,i) = mulreal(gel(x,i), gel(y,j));
4056 : }
4057 21630 : return q;
4058 : }
4059 :
4060 : /*******************************************************************/
4061 : /* */
4062 : /* LOGICAL OPERATIONS */
4063 : /* */
4064 : /*******************************************************************/
4065 : static long
4066 109043713 : _egal_i(GEN x, GEN y)
4067 : {
4068 109043713 : x = simplify_shallow(x);
4069 109043713 : y = simplify_shallow(y);
4070 109043713 : if (typ(y) == t_INT)
4071 : {
4072 108002899 : if (equali1(y)) return gequal1(x);
4073 62348244 : if (equalim1(y)) return gequalm1(x);
4074 : }
4075 1040814 : else if (typ(x) == t_INT)
4076 : {
4077 140 : if (equali1(x)) return gequal1(y);
4078 91 : if (equalim1(x)) return gequalm1(y);
4079 : }
4080 63257437 : return gequal(x, y);
4081 : }
4082 : static long
4083 109043713 : _egal(GEN x, GEN y)
4084 109043713 : { pari_sp av = avma; return gc_long(av, _egal_i(x, y)); }
4085 :
4086 : GEN
4087 6328264 : glt(GEN x, GEN y) { return gcmp(x,y)<0? gen_1: gen_0; }
4088 :
4089 : GEN
4090 7628237 : gle(GEN x, GEN y) { return gcmp(x,y)<=0? gen_1: gen_0; }
4091 :
4092 : GEN
4093 248443 : gge(GEN x, GEN y) { return gcmp(x,y)>=0? gen_1: gen_0; }
4094 :
4095 : GEN
4096 2375133 : ggt(GEN x, GEN y) { return gcmp(x,y)>0? gen_1: gen_0; }
4097 :
4098 : GEN
4099 47673975 : geq(GEN x, GEN y) { return _egal(x,y)? gen_1: gen_0; }
4100 :
4101 : GEN
4102 61369738 : gne(GEN x, GEN y) { return _egal(x,y)? gen_0: gen_1; }
4103 :
4104 : GEN
4105 606149 : gnot(GEN x) { return gequal0(x)? gen_1: gen_0; }
4106 :
4107 : /*******************************************************************/
4108 : /* */
4109 : /* FORMAL SIMPLIFICATIONS */
4110 : /* */
4111 : /*******************************************************************/
4112 :
4113 : GEN
4114 11020 : geval_gp(GEN x, GEN t)
4115 : {
4116 11020 : long lx, i, tx = typ(x);
4117 : pari_sp av;
4118 : GEN y, z;
4119 :
4120 11020 : if (is_const_t(tx) || tx==t_VECSMALL) return gcopy(x);
4121 10999 : switch(tx)
4122 : {
4123 10992 : case t_STR:
4124 10992 : return localvars_read_str(GSTR(x),t);
4125 :
4126 0 : case t_POLMOD:
4127 0 : av = avma;
4128 0 : return gc_upto(av, gmodulo(geval_gp(gel(x,2),t),
4129 0 : geval_gp(gel(x,1),t)));
4130 :
4131 7 : case t_POL:
4132 7 : lx=lg(x); if (lx==2) return gen_0;
4133 7 : z = fetch_var_value(varn(x),t);
4134 7 : if (!z) return RgX_copy(x);
4135 7 : av = avma; y = geval_gp(gel(x,lx-1),t);
4136 14 : for (i=lx-2; i>1; i--)
4137 7 : y = gadd(geval_gp(gel(x,i),t), gmul(z,y));
4138 7 : return gc_upto(av, y);
4139 :
4140 0 : case t_SER:
4141 0 : pari_err_IMPL( "evaluation of a power series");
4142 :
4143 0 : case t_RFRAC:
4144 0 : av = avma;
4145 0 : return gc_upto(av, gdiv(geval_gp(gel(x,1),t), geval_gp(gel(x,2),t)));
4146 :
4147 0 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
4148 0 : pari_APPLY_same(geval_gp(gel(x,i),t));
4149 :
4150 0 : case t_CLOSURE:
4151 0 : if (closure_arity(x) || closure_is_variadic(x))
4152 0 : pari_err_IMPL("eval on functions with parameters");
4153 0 : return closure_evalres(x);
4154 : }
4155 0 : pari_err_TYPE("geval",x);
4156 : return NULL; /* LCOV_EXCL_LINE */
4157 : }
4158 : GEN
4159 0 : geval(GEN x) { return geval_gp(x,NULL); }
4160 :
4161 : GEN
4162 523432520 : simplify_shallow(GEN x)
4163 : {
4164 : long v, lx;
4165 : GEN y, z;
4166 523432520 : if (!x) pari_err_BUG("simplify, NULL input");
4167 :
4168 523432520 : switch(typ(x))
4169 : {
4170 432705470 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
4171 : case t_PADIC: case t_QFB: case t_LIST: case t_STR: case t_VECSMALL:
4172 : case t_CLOSURE: case t_ERROR: case t_INFINITY:
4173 432705470 : return x;
4174 10178 : case t_COMPLEX: return isintzero(gel(x,2))? gel(x,1): x;
4175 399 : case t_QUAD: return isintzero(gel(x,3))? gel(x,2): x;
4176 :
4177 196134 : case t_POLMOD: y = cgetg(3,t_POLMOD);
4178 196134 : z = gel(x,1); v = varn(z); z = simplify_shallow(z);
4179 196134 : if (typ(z) != t_POL || varn(z) != v) z = scalarpol_shallow(z, v);
4180 196134 : gel(y,1) = z;
4181 196134 : gel(y,2) = simplify_shallow(gel(x,2)); return y;
4182 :
4183 83585860 : case t_POL:
4184 83585860 : lx = lg(x);
4185 83585860 : if (lx==2) return gen_0;
4186 72126985 : if (lx==3) return simplify_shallow(gel(x,2));
4187 230935998 : pari_APPLY_pol(simplify_shallow(gel(x,i)));
4188 :
4189 651 : case t_SER:
4190 566888 : pari_APPLY_ser(simplify_shallow(gel(x,i)));
4191 :
4192 648366 : case t_RFRAC: y = cgetg(3,t_RFRAC);
4193 648366 : gel(y,1) = simplify_shallow(gel(x,1));
4194 648366 : z = simplify_shallow(gel(x,2));
4195 648366 : if (typ(z) != t_POL) return gdiv(gel(y,1), z);
4196 648366 : gel(y,2) = z; return y;
4197 :
4198 6285462 : case t_VEC: case t_COL: case t_MAT:
4199 17910341 : pari_APPLY_same(simplify_shallow(gel(x,i)));
4200 : }
4201 0 : pari_err_BUG("simplify_shallow, type unknown");
4202 : return NULL; /* LCOV_EXCL_LINE */
4203 : }
4204 :
4205 : GEN
4206 11592339 : simplify(GEN x)
4207 : {
4208 11592339 : pari_sp av = avma;
4209 11592339 : GEN y = simplify_shallow(x);
4210 11592339 : return av == avma ? gcopy(y): gc_GEN(av, y);
4211 : }
4212 :
4213 : /*******************************************************************/
4214 : /* */
4215 : /* EVALUATION OF SOME SIMPLE OBJECTS */
4216 : /* */
4217 : /*******************************************************************/
4218 : /* q is a real symmetric matrix, x a RgV. Horner-type evaluation of q(x)
4219 : * using (n^2+3n-2)/2 mul */
4220 : GEN
4221 17023 : qfeval(GEN q, GEN x)
4222 : {
4223 17023 : pari_sp av = avma;
4224 17023 : long i, j, l = lg(q);
4225 : GEN z;
4226 17023 : if (lg(x) != l) pari_err_DIM("qfeval");
4227 17016 : if (l==1) return gen_0;
4228 17016 : if (lgcols(q) != l) pari_err_DIM("qfeval");
4229 : /* l = lg(x) = lg(q) > 1 */
4230 17009 : z = gmul(gcoeff(q,1,1), gsqr(gel(x,1)));
4231 74015 : for (i=2; i<l; i++)
4232 : {
4233 57006 : GEN c = gel(q,i), s;
4234 57006 : if (isintzero(gel(x,i))) continue;
4235 43356 : s = gmul(gel(c,1), gel(x,1));
4236 131810 : for (j=2; j<i; j++) s = gadd(s, gmul(gel(c,j),gel(x,j)));
4237 43356 : s = gadd(gshift(s,1), gmul(gel(c,i),gel(x,i)));
4238 43356 : z = gadd(z, gmul(gel(x,i), s));
4239 : }
4240 17009 : return gc_upto(av,z);
4241 : }
4242 :
4243 : static GEN
4244 347816 : qfbeval(GEN q, GEN z)
4245 : {
4246 347816 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3), x = gel(z,1), y = gel(z,2);
4247 347816 : pari_sp av = avma;
4248 347816 : A = gadd(gmul(x, gadd(gmul(a,x), gmul(b,y))), gmul(c, gsqr(y)));
4249 347816 : return gc_upto(av, A);
4250 : }
4251 : static GEN
4252 7 : qfbevalb(GEN q, GEN z, GEN z2)
4253 : {
4254 7 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3);
4255 7 : GEN x = gel(z,1), y = gel(z,2);
4256 7 : GEN X = gel(z2,1), Y = gel(z2,2);
4257 7 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4258 7 : pari_sp av = avma;
4259 : /* a2 x X + b (x Y + X y) + c2 y Y */
4260 7 : A = gadd(gmul(x, gadd(gmul(a2,X), gmul(b,Y))),
4261 : gmul(y, gadd(gmul(c2,Y), gmul(b,X))));
4262 7 : return gc_upto(av, gmul2n(A, -1));
4263 : }
4264 : static GEN
4265 21 : qfb_ZM_apply_i(GEN q, GEN M)
4266 : {
4267 21 : GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), d = gel(q,4);
4268 21 : GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
4269 21 : GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
4270 21 : GEN by = mulii(b,y), bt = mulii(b,t), bz = mulii(b,z);
4271 21 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4272 :
4273 21 : GEN A1 = mulii(x, addii(mulii(a,x), by));
4274 21 : GEN A2 = mulii(c, sqri(y));
4275 21 : GEN B1 = mulii(x, addii(mulii(a2,z), bt));
4276 21 : GEN B2 = mulii(y, addii(mulii(c2,t), bz));
4277 21 : GEN C1 = mulii(z, addii(mulii(a,z), bt));
4278 21 : GEN C2 = mulii(c, sqri(t));
4279 21 : GEN m = sqri(subii(mulii(x,t), mulii(y,z)));
4280 21 : if (signe(m)==0) pari_err_INV("qfapply",M);
4281 21 : retmkqfb(addii(A1,A2), addii(B1,B2), addii(C1, C2), mulii(d, m));
4282 : }
4283 :
4284 : GEN
4285 21 : qfb_ZM_apply(GEN q, GEN M)
4286 : {
4287 21 : pari_sp av = avma;
4288 21 : return gc_upto(av, qfb_ZM_apply_i(q, M));
4289 : }
4290 :
4291 : static GEN
4292 348026 : qfnorm0(GEN q, GEN x)
4293 : {
4294 348026 : if (!q) switch(typ(x))
4295 : {
4296 7 : case t_VEC: case t_COL: return RgV_dotsquare(x);
4297 7 : case t_MAT: return gram_matrix(x);
4298 7 : default: pari_err_TYPE("qfeval",x);
4299 : }
4300 348005 : switch(typ(q))
4301 : {
4302 161 : case t_MAT: break;
4303 347837 : case t_QFB:
4304 347837 : if (lg(x) == 3) switch(typ(x))
4305 : {
4306 347816 : case t_VEC:
4307 347816 : case t_COL: return qfbeval(q,x);
4308 21 : case t_MAT: if (RgM_is_ZM(x)) return qfb_ZM_apply(q,x);
4309 0 : default: pari_err_TYPE("qfeval",x);
4310 : }
4311 7 : default: pari_err_TYPE("qfeval",q);
4312 : }
4313 161 : switch(typ(x))
4314 : {
4315 154 : case t_VEC: case t_COL: break;
4316 7 : case t_MAT: return qf_RgM_apply(q, x);
4317 0 : default: pari_err_TYPE("qfeval",x);
4318 : }
4319 154 : return qfeval(q,x);
4320 : }
4321 : /* obsolete, use qfeval0 */
4322 : GEN
4323 0 : qfnorm(GEN x, GEN q) { return qfnorm0(q,x); }
4324 :
4325 : /* assume q is square, x~ * q * y using n^2+n mul */
4326 : GEN
4327 567 : qfevalb(GEN q, GEN x, GEN y)
4328 : {
4329 567 : pari_sp av = avma;
4330 567 : long l = lg(q);
4331 567 : if (lg(x) != l || lg(y) != l) pari_err_DIM("qfevalb");
4332 560 : return gc_upto(av, RgV_dotproduct(RgV_RgM_mul(x,q), y));
4333 : }
4334 :
4335 : /* obsolete, use qfeval0 */
4336 : GEN
4337 0 : qfbil(GEN x, GEN y, GEN q)
4338 : {
4339 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfbil",x);
4340 0 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfbil",y);
4341 0 : if (!q) {
4342 0 : if (lg(x) != lg(y)) pari_err_DIM("qfbil");
4343 0 : return RgV_dotproduct(x,y);
4344 : }
4345 0 : if (typ(q) != t_MAT) pari_err_TYPE("qfbil",q);
4346 0 : return qfevalb(q,x,y);
4347 : }
4348 : GEN
4349 348082 : qfeval0(GEN q, GEN x, GEN y)
4350 : {
4351 348082 : if (!y) return qfnorm0(q, x);
4352 56 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfeval",x);
4353 42 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfeval",y);
4354 42 : if (!q) {
4355 14 : if (lg(x) != lg(y)) pari_err_DIM("qfeval");
4356 7 : return RgV_dotproduct(x,y);
4357 : }
4358 28 : switch(typ(q))
4359 : {
4360 21 : case t_MAT: break;
4361 7 : case t_QFB:
4362 7 : if (lg(x) == 3 && lg(y) == 3) return qfbevalb(q,x,y);
4363 0 : default: pari_err_TYPE("qfeval",q);
4364 : }
4365 21 : return qfevalb(q,x,y);
4366 : }
4367 :
4368 : /* q a hermitian complex matrix, x a RgV */
4369 : GEN
4370 168 : hqfeval(GEN q, GEN x)
4371 : {
4372 168 : pari_sp av = avma;
4373 168 : long i, j, l = lg(q);
4374 : GEN z, xc;
4375 :
4376 168 : if (lg(x) != l) pari_err_DIM("hqfeval");
4377 168 : if (l==1) return gen_0;
4378 168 : if (lgcols(q) != l) pari_err_DIM("hqfeval");
4379 168 : if (l == 2) return gc_upto(av, gmul(gcoeff(q,1,1), gnorm(gel(x,1))));
4380 : /* l = lg(x) = lg(q) > 2 */
4381 168 : xc = conj_i(x);
4382 168 : z = mulreal(gcoeff(q,2,1), gmul(gel(x,2),gel(xc,1)));
4383 1092 : for (i=3;i<l;i++)
4384 5418 : for (j=1;j<i;j++)
4385 4494 : z = gadd(z, mulreal(gcoeff(q,i,j), gmul(gel(x,i),gel(xc,j))));
4386 168 : z = gshift(z,1);
4387 1428 : for (i=1;i<l;i++) z = gadd(z, gmul(gcoeff(q,i,i), gnorm(gel(x,i))));
4388 168 : return gc_upto(av,z);
4389 : }
4390 :
4391 : static void
4392 504718 : init_qf_apply(GEN q, GEN M, long *l)
4393 : {
4394 504718 : long k = lg(M);
4395 504718 : *l = lg(q);
4396 504718 : if (*l == 1) { if (k == 1) return; }
4397 504711 : else { if (k != 1 && lgcols(M) == *l) return; }
4398 0 : pari_err_DIM("qf_RgM_apply");
4399 : }
4400 : /* Return X = M'.q.M, assuming q is a symmetric matrix and M is a
4401 : * matrix of compatible dimensions. X_ij are X_ji identical, not copies */
4402 : GEN
4403 7663 : qf_RgM_apply(GEN q, GEN M)
4404 : {
4405 7663 : pari_sp av = avma;
4406 7663 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4407 7663 : return gc_upto(av, RgM_transmultosym(M, RgM_mul(q, M)));
4408 : }
4409 : GEN
4410 497055 : qf_ZM_apply(GEN q, GEN M)
4411 : {
4412 497055 : pari_sp av = avma;
4413 497055 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4414 : /* FIXME: ZM_transmultosym is asymptotically slow, so choose some random
4415 : * threshold defaulting to default implementation: this is suboptimal but
4416 : * has the right complexity in the dimension. Should implement M'*q*M as an
4417 : * atomic operation with the right complexity, see ZM_mul_i. */
4418 497048 : if (l > 20)
4419 182 : M = ZM_mul(shallowtrans(M), ZM_mul(q, M));
4420 : else
4421 496866 : M = ZM_transmultosym(M, ZM_mul(q, M));
4422 497048 : return gc_upto(av, M);
4423 : }
4424 :
4425 : GEN
4426 13415490 : poleval(GEN x, GEN y)
4427 : {
4428 13415490 : long i, j, imin, M = 1, tx = typ(x);
4429 13415490 : pari_sp av0 = avma, av;
4430 : GEN t, t2, r, s;
4431 :
4432 13415490 : if (is_scalar_t(tx)) return gcopy(x);
4433 13236067 : switch(tx)
4434 : {
4435 13089724 : case t_POL:
4436 13089724 : if (typ(y) == t_POL && varn(x) == varn(y))
4437 : {
4438 188188 : y = RgX_deflate_max(y, &M);
4439 188188 : if (degpol(y) == 1)
4440 : {
4441 126154 : GEN a = gel(y,3), b = gel(y,2); /* y = a t + b */
4442 126154 : if (isint1(a))
4443 : {
4444 95914 : y = RgX_Rg_translate(x, b);
4445 95914 : if (M == 1) return y;
4446 : }
4447 : else
4448 30240 : y = RgX_affine(x, a, b);
4449 112644 : if (M != 1) y = RgX_inflate(y, M);
4450 112644 : return gc_GEN(av0, y);
4451 : }
4452 : }
4453 :
4454 12963570 : i = lg(x)-1; imin = 2; break;
4455 :
4456 1568 : case t_RFRAC:
4457 1568 : return gc_upto(av0, gdiv(poleval(gel(x,1),y),
4458 1568 : poleval(gel(x,2),y)));
4459 :
4460 144775 : case t_VEC: case t_COL:
4461 144775 : i = lg(x)-1; imin = 1; break;
4462 0 : default: pari_err_TYPE("poleval",x);
4463 : return NULL; /* LCOV_EXCL_LINE */
4464 : }
4465 13108345 : if (i<=imin) return (i==imin)? gmul(gel(x,imin),Rg_get_1(y)): Rg_get_0(y);
4466 12699366 : if (isintzero(y)) return gcopy(gel(x,imin));
4467 :
4468 12691916 : t = gel(x,i); i--;
4469 12691916 : if (typ(y)!=t_COMPLEX)
4470 : {
4471 : #if 0 /* standard Horner's rule */
4472 : for ( ; i>=imin; i--)
4473 : t = gadd(gmul(t,y),gel(x,i));
4474 : #endif
4475 : /* specific attention to sparse polynomials */
4476 12586252 : av = avma;
4477 50006548 : for ( ; i>=imin; i=j-1)
4478 : {
4479 57739151 : for (j=i; isexactzero(gel(x,j)); j--)
4480 20318855 : if (j==imin)
4481 : {
4482 11073668 : if (i!=j) y = gpowgs(y, i-j+1);
4483 11073668 : y = gmul(t, y);
4484 11073668 : if (M == 1) return gc_upto(av0, y);
4485 0 : return gc_GEN(av0, RgX_inflate(y, M));
4486 : }
4487 37420296 : r = (i==j)? y: gpowgs(y, i-j+1);
4488 37420296 : t = gadd(gmul(t,r), gel(x,j));
4489 37420296 : if (gc_needed(av0,2))
4490 : {
4491 109 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4492 109 : t = gc_upto(av, t);
4493 : }
4494 : }
4495 1512584 : if (M == 1) return gc_upto(av0, t);
4496 7 : return gc_GEN(av0, RgX_inflate(t, M));
4497 : }
4498 :
4499 105664 : t2 = gel(x,i); i--; r = gtrace(y); s = gneg_i(gnorm(y));
4500 105664 : av = avma;
4501 5001242 : for ( ; i>=imin; i--)
4502 : {
4503 4895578 : GEN p = gadd(t2, gmul(r, t));
4504 4895578 : t2 = gadd(gel(x,i), gmul(s, t)); t = p;
4505 4895578 : if (gc_needed(av0,2))
4506 : {
4507 0 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4508 0 : (void)gc_all(av, 2, &t, &t2);
4509 : }
4510 : }
4511 105664 : return gc_upto(av0, gadd(t2, gmul(y, t)));
4512 : }
4513 :
4514 : /* Evaluate a polynomial using Horner. Unstable!
4515 : * If ui != NULL, ui = 1/u, evaluate P(1/u)*u^(deg P): useful for |u|>1 */
4516 : GEN
4517 2472913 : RgX_cxeval(GEN T, GEN u, GEN ui)
4518 : {
4519 2472913 : pari_sp ltop = avma;
4520 : GEN S;
4521 2472913 : long n, lim = lg(T)-1;
4522 2472913 : if (lim == 1) return gen_0;
4523 2472913 : if (lim == 2) return gcopy(gel(T,2));
4524 2471755 : if (!ui)
4525 : {
4526 1716288 : n = lim; S = gel(T,n);
4527 16866998 : for (n--; n >= 2; n--) S = gadd(gmul(u,S), gel(T,n));
4528 : }
4529 : else
4530 : {
4531 755467 : n = 2; S = gel(T,2);
4532 4623975 : for (n++; n <= lim; n++) S = gadd(gmul(ui,S), gel(T,n));
4533 755467 : S = gmul(gpowgs(u, lim-2), S);
4534 : }
4535 2471755 : return gc_upto(ltop, S);
4536 : }
4537 :
4538 : GEN
4539 63 : RgXY_cxevalx(GEN x, GEN u, GEN ui)
4540 : {
4541 427 : pari_APPLY_pol(typ(gel(x,i))==t_POL? RgX_cxeval(gel(x,i), u, ui): gel(x,i));
4542 : }
|
