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 26579 : recvar(hashtable *h, GEN x)
31 : {
32 26579 : long i = 1, lx = lg(x);
33 : void *v;
34 26579 : switch(typ(x))
35 : {
36 6090 : case t_POL: case t_SER:
37 6090 : v = (void*)varn(x);
38 6090 : if (!hash_search(h, v)) hash_insert(h, v, NULL);
39 6090 : i = 2; break;
40 1043 : case t_POLMOD: case t_RFRAC:
41 1043 : 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 19432 : default:
46 19432 : return;
47 : }
48 32669 : for (; i < lx; i++) recvar(h, gel(x,i));
49 : }
50 :
51 : GEN
52 1057 : variables_vecsmall(GEN x)
53 : {
54 1057 : hashtable *h = hash_create_ulong(100, 1);
55 1057 : recvar(h, x);
56 1057 : return vars_sort_inplace(hash_keys(h));
57 : }
58 :
59 : GEN
60 42 : variables_vec(GEN x)
61 42 : { return x? vars_to_RgXV(variables_vecsmall(x)): gpolvar(NULL); }
62 :
63 : long
64 128660678 : gvar(GEN x)
65 : {
66 : long i, v, w, lx;
67 128660678 : switch(typ(x))
68 : {
69 49679729 : case t_POL: case t_SER: return varn(x);
70 176066 : case t_POLMOD: return varn(gel(x,1));
71 14546564 : case t_RFRAC: return varn(gel(x,2));
72 3675896 : case t_VEC: case t_COL: case t_MAT:
73 3675896 : lx = lg(x); break;
74 14 : case t_LIST:
75 14 : x = list_data(x);
76 14 : lx = x? lg(x): 1; break;
77 60582409 : default:
78 60582409 : return NO_VARIABLE;
79 : }
80 3675910 : v = NO_VARIABLE;
81 32926441 : for (i=1; i < lx; i++) { w = gvar(gel(x,i)); if (varncmp(w,v) < 0) v = w; }
82 3675921 : 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 9968 : var2_aux(GEN T, GEN A)
88 : {
89 9968 : long a = gvar2(T);
90 9968 : long b = (typ(A) == t_POL && varn(A) == varn(T))? gvar2(A): gvar(A);
91 9968 : if (varncmp(a, b) > 0) a = b;
92 9968 : return a;
93 : }
94 : static long
95 7035 : var2_rfrac(GEN x) { return var2_aux(gel(x,2), gel(x,1)); }
96 : static long
97 2933 : 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 55160 : gvar9(GEN x)
103 55160 : { 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 20639976 : gvar2(GEN x)
108 : {
109 : long i, v, w;
110 20639976 : switch(typ(x))
111 : {
112 56 : case t_POLMOD:
113 56 : return var2_polmod(x);
114 16912 : case t_POL: case t_SER:
115 16912 : v = NO_VARIABLE;
116 71029 : for (i=2; i < lg(x); i++) {
117 54117 : w = gvar9(gel(x,i));
118 54117 : if (varncmp(w,v) < 0) v=w;
119 : }
120 16912 : return v;
121 7035 : case t_RFRAC:
122 7035 : 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 20615924 : return NO_VARIABLE;
132 : }
133 :
134 : /*******************************************************************/
135 : /* */
136 : /* PRECISION OF SCALAR OBJECTS */
137 : /* */
138 : /*******************************************************************/
139 : static long
140 10084769 : prec0(long e) { return (e < 0)? nbits2prec(-e): LOWDEFAULTPREC; }
141 : static long
142 934734048 : 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 1346570 : precrealexact(GEN x, GEN y)
146 : {
147 1346570 : long lx, ey = gexpo(y), ex, e;
148 1346571 : if (ey == -(long)HIGHEXPOBIT) return precREAL(x);
149 643893 : ex = expo(x);
150 643893 : e = ey - ex;
151 643893 : if (!signe(x)) return prec0((e >= 0)? -e: ex);
152 643844 : lx = realprec(x);
153 643844 : return (e > 0)? lx + nbits2extraprec(e): lx;
154 : }
155 : static long
156 22844516 : precCOMPLEX(GEN z)
157 : { /* ~ precision(|x| + |y|) */
158 22844516 : GEN x = gel(z,1), y = gel(z,2);
159 : long e, ex, ey, lz, lx, ly;
160 22844516 : if (typ(x) != t_REAL) {
161 2255174 : if (typ(y) != t_REAL) return 0;
162 1334295 : return precrealexact(y, x);
163 : }
164 20589342 : if (typ(y) != t_REAL) return precrealexact(x, y);
165 : /* x, y are t_REALs, cf addrr_sign */
166 20577063 : ex = expo(x);
167 20577063 : ey = expo(y);
168 20577063 : e = ey - ex;
169 20577063 : if (!signe(x)) {
170 554163 : if (!signe(y)) return prec0( minss(ex,ey) );
171 554065 : if (e <= 0) return prec0(ex);
172 553850 : lz = nbits2prec(e);
173 553846 : ly = realprec(y); if (lz > ly) lz = ly;
174 553846 : return lz;
175 : }
176 20022900 : if (!signe(y)) {
177 79757 : if (e >= 0) return prec0(ey);
178 79750 : lz = nbits2prec(-e);
179 79750 : lx = realprec(x); if (lz > lx) lz = lx;
180 79750 : return lz;
181 : }
182 19943143 : if (e < 0) { swap(x, y); e = -e; }
183 19943143 : lx = realprec(x);
184 19943143 : ly = realprec(y);
185 19943143 : if (e) {
186 16951724 : long d = nbits2extraprec(e), l = ly-d;
187 16951714 : return (l > lx)? lx + d: ly;
188 : }
189 2991419 : return minss(lx, ly);
190 : }
191 : long
192 947242944 : precision(GEN z)
193 : {
194 947242944 : switch(typ(z))
195 : {
196 930264826 : case t_REAL: return precREAL(z);
197 16935379 : case t_COMPLEX: return precCOMPLEX(z);
198 : }
199 42739 : return 0;
200 : }
201 :
202 : long
203 14441439 : gprecision(GEN x)
204 : {
205 : long i, k, l;
206 :
207 14441439 : switch(typ(x))
208 : {
209 3766936 : case t_REAL: return precREAL(x);
210 5909136 : case t_COMPLEX: return precCOMPLEX(x);
211 877624 : case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT:
212 : case t_PADIC: case t_QUAD: case t_POLMOD:
213 877624 : return 0;
214 :
215 154 : case t_POL: case t_SER:
216 154 : k = LONG_MAX;
217 525 : for (i=lg(x)-1; i>1; i--)
218 : {
219 371 : l = gprecision(gel(x,i));
220 371 : if (l && l<k) k = l;
221 : }
222 154 : return (k==LONG_MAX)? 0: k;
223 3887672 : case t_VEC: case t_COL: case t_MAT:
224 3887672 : k = LONG_MAX;
225 14349679 : for (i=lg(x)-1; i>0; i--)
226 : {
227 10461954 : l = gprecision(gel(x,i));
228 10462007 : if (l && l<k) k = l;
229 : }
230 3887725 : 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 0 : return 0;
242 : }
243 :
244 : static long
245 413 : vec_padicprec_relative(GEN x, long imin)
246 : {
247 : long s, t, i;
248 1288 : for (s=LONG_MAX, i=lg(x)-1; i>=imin; i--)
249 : {
250 875 : t = padicprec_relative(gel(x,i)); if (t<s) s = t;
251 : }
252 413 : 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 2359 : padicprec_relative(GEN x)
258 : {
259 2359 : switch(typ(x))
260 : {
261 427 : case t_INT: case t_FRAC:
262 427 : return LONG_MAX;
263 1519 : case t_PADIC:
264 1519 : return signe(gel(x,4))? precp(x): 0;
265 238 : case t_POLMOD: case t_VEC: case t_COL: case t_MAT:
266 238 : return vec_padicprec_relative(x, 1);
267 175 : case t_POL: case t_SER:
268 175 : 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(gel(x,2),p)) pari_err_MODULUS("padicprec", gel(x,2), 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 102268 : poldegree(GEN x, long v)
370 : {
371 102268 : const long DEGREE0 = -LONG_MAX;
372 102268 : long tx = typ(x), lx,w,i,d;
373 :
374 102268 : if (is_scalar_t(tx)) return gequal0(x)? DEGREE0: 0;
375 101914 : switch(tx)
376 : {
377 101823 : case t_POL:
378 101823 : if (!signe(x)) return DEGREE0;
379 101816 : w = varn(x);
380 101816 : if (v < 0 || v == w) return degpol(x);
381 144 : if (varncmp(v, w) < 0) return 0;
382 144 : lx = lg(x); d = DEGREE0;
383 684 : for (i=2; i<lx; i++)
384 : {
385 540 : long e = poldegree(gel(x,i), v);
386 540 : if (e > d) d = e;
387 : }
388 144 : return d;
389 :
390 91 : case t_RFRAC:
391 : {
392 91 : GEN a = gel(x,1), b = gel(x,2);
393 91 : if (gequal0(a)) return DEGREE0;
394 84 : if (v < 0)
395 : {
396 84 : v = varn(b); d = -degpol(b);
397 84 : if (typ(a) == t_POL && varn(a) == v) d += degpol(a);
398 84 : 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 30923 : gppoldegree(GEN x, long v)
408 : {
409 30923 : long d = poldegree(x,v);
410 30923 : 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 503764 : RgX_degree(GEN x, long v)
416 : {
417 503764 : long tx = typ(x), lx, w, i, d;
418 :
419 503764 : if (is_scalar_t(tx)) return gequal0(x)? -1: 0;
420 303191 : switch(tx)
421 : {
422 303191 : case t_POL:
423 303191 : if (!signe(x)) return -1;
424 303170 : w = varn(x);
425 303170 : if (v == w) return degpol(x);
426 107442 : if (varncmp(v, w) < 0) return 0;
427 107443 : lx = lg(x); d = -1;
428 485111 : for (i=2; i<lx; i++)
429 : {
430 377668 : long e = RgX_degree(gel(x,i), v);
431 377668 : if (e > d) d = e;
432 : }
433 107443 : 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 6790 : degree(GEN x)
447 : {
448 6790 : 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 gerepileupto(av, x);
484 : }
485 :
486 : /* returns 1 if there's a real component in the structure, 0 otherwise */
487 : int
488 14357 : isinexactreal(GEN x)
489 : {
490 : long i;
491 14357 : switch(typ(x))
492 : {
493 1246 : case t_REAL: return 1;
494 2597 : case t_COMPLEX: return (typ(gel(x,1))==t_REAL || typ(gel(x,2))==t_REAL);
495 :
496 9926 : case t_INT: case t_INTMOD: case t_FRAC:
497 : case t_FFELT: case t_PADIC: case t_QUAD:
498 9926 : 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 1859986 : isrealappr(GEN x, long e)
518 : {
519 : long i;
520 1859986 : switch(typ(x))
521 : {
522 695551 : case t_INT: case t_REAL: case t_FRAC:
523 695551 : return 1;
524 1164438 : case t_COMPLEX:
525 1164438 : 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 132117986 : isinexact(GEN x)
547 : {
548 : long lx, i;
549 :
550 132117986 : switch(typ(x))
551 : {
552 583010 : case t_REAL: case t_PADIC: case t_SER:
553 583010 : return 1;
554 89102323 : case t_INT: case t_INTMOD: case t_FFELT: case t_FRAC:
555 : case t_QFB:
556 89102323 : return 0;
557 2411120 : case t_COMPLEX: case t_QUAD: case t_RFRAC: case t_POLMOD:
558 2411120 : return isinexact(gel(x,1)) || isinexact(gel(x,2));
559 40003039 : case t_POL:
560 125419379 : for (i=lg(x)-1; i>1; i--)
561 85586492 : if (isinexact(gel(x,i))) return 1;
562 39832887 : return 0;
563 18494 : case t_VEC: case t_COL: case t_MAT:
564 22134 : for (i=lg(x)-1; i>0; i--)
565 21742 : if (isinexact(gel(x,i))) return 1;
566 392 : 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 177345 : is_realext(GEN x)
613 177345 : { long t = typ(x);
614 177345 : 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 28 : _quotrs(GEN x, long y)
637 : {
638 28 : GEN q = divrs(x,y), f = floorr(q);
639 28 : if (y < 0 && signe(subir(f,q))) f = addiu(f, 1);
640 28 : 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 gerepileupto(av, _quot(x, y)); }
680 : static GEN
681 14 : quotrs(GEN x, long y)
682 14 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotrs(x,y)); }
683 : static GEN
684 301 : quotfs(GEN x, long s)
685 301 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotfs(x,s)); }
686 : static GEN
687 35 : quotsr(long x, GEN y)
688 35 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotsr(x, y)); }
689 : static GEN
690 35 : quotsf(long x, GEN y)
691 35 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotsf(x, y)); }
692 : static GEN
693 35 : quotsq(long x, GEN y)
694 35 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotsq(x, y)); }
695 : static GEN
696 14 : quotqs(GEN x, long y)
697 14 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotqs(x, y)); }
698 : static GEN
699 7 : quotfi(GEN x, GEN y)
700 7 : { pari_sp av = avma; return gerepileuptoleaf(av, _quotfi(x, y)); }
701 : static GEN
702 7 : quotri(GEN x, GEN y)
703 7 : { pari_sp av = avma; return gerepileuptoleaf(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 gerepileuptoleaf(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 gerepileuptoleaf(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 gerepileupto(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 4846595 : modRr_i(GEN x, GEN y, GEN iy)
746 : {
747 : GEN q, f;
748 : long e;
749 4846595 : if (isintzero(x)) return gen_0;
750 4846591 : q = gmul(x, iy); /* t_REAL */
751 :
752 4846614 : e = expo(q);
753 4846614 : if (e >= 0 && nbits2prec(e+1) > realprec(q)) return NULL;
754 4846612 : f = floorr(q);
755 4846524 : if (signe(y) < 0 && signe(subri(q,f))) f = addiu(f, 1);
756 4846587 : return signe(f)? gsub(x, mulir(f,y)): x;
757 : }
758 :
759 : GEN
760 46136635 : gmod(GEN x, GEN y)
761 : {
762 : pari_sp av;
763 : long ty, tx;
764 : GEN z;
765 :
766 46136635 : tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gmodsg(itos(x),y);
767 1146398 : ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gmodgs(x,itos(y));
768 1745502 : if (is_matvec_t(tx)) pari_APPLY_same(gmod(gel(x,i), y));
769 808182 : 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 : 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 gerepileupto(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 gerepileupto(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 22028249 : gmodgs(GEN x, long y)
801 : {
802 : ulong u;
803 22028249 : long i, tx = typ(x);
804 : GEN z;
805 43822471 : if (is_matvec_t(tx)) pari_APPLY_same(gmodgs(gel(x,i), y));
806 19649951 : if (!y) pari_err_INV("gmodgs",gen_0);
807 19649951 : switch(tx)
808 : {
809 19566486 : 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 82028 : case t_FRAC:
819 82028 : u = (ulong)labs(y);
820 82028 : return utoi( Fl_div(umodiu(gel(x,1), u),
821 82028 : 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 gerepileupto(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 44990237 : gmodsg(long x, GEN y)
837 : {
838 44990237 : switch(typ(y))
839 : {
840 44989880 : 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 gerepileupto(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 15531 : gdvd(GEN x, GEN y)
859 : {
860 15531 : pari_sp av = avma;
861 15531 : return gc_bool(av, gequal0( gmod(x,y) ));
862 : }
863 :
864 : GEN
865 774631 : gmodulss(long x, long y)
866 : {
867 774631 : if (!y) pari_err_INV("%",gen_0);
868 774624 : y = labs(y);
869 774624 : retmkintmod(utoi(umodsu(x, y)), utoipos(y));
870 : }
871 : GEN
872 1077111 : gmodulsg(long x, GEN y)
873 : {
874 1077111 : switch(typ(y))
875 : {
876 828014 : case t_INT:
877 828014 : if (!is_bigint(y)) return gmodulss(x,itos(y));
878 62388 : retmkintmod(modsi(x,y), absi(y));
879 249090 : case t_POL:
880 249090 : if (!signe(y)) pari_err_INV("%", y);
881 249083 : 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 1600113 : gmodulo(GEN x,GEN y)
888 : {
889 1600113 : long tx = typ(x), vx, vy;
890 1600113 : if (tx == t_INT && !is_bigint(x)) return gmodulsg(itos(x), y);
891 1069344 : if (is_matvec_t(tx)) pari_APPLY_same(gmodulo(gel(x,i), y));
892 518030 : switch(typ(y))
893 : {
894 98212 : case t_INT:
895 98212 : if (!is_const_t(tx)) return gmul(x, gmodulsg(1,y));
896 98163 : if (tx == t_INTMOD) return gmod(x,y);
897 98156 : retmkintmod(Rg_to_Fp(x,y), absi(y));
898 419818 : case t_POL:
899 419818 : vx = gvar(x); vy = varn(y);
900 419818 : if (varncmp(vy, vx) > 0) return gmul(x, gmodulsg(1,y));
901 416066 : if (vx == vy && tx == t_POLMOD) return grem(x,y);
902 402668 : 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 6200418 : gdivent(GEN x, GEN y)
915 : {
916 : long tx, ty;
917 6200418 : 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 : 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 292780 : gdiventgs(GEN x, long y)
945 : {
946 292780 : switch(typ(x))
947 : {
948 244019 : case t_INT: return truedivis(x,y);
949 14 : 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 282969 : 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 6198230 : gdiventsg(long x, GEN y)
961 : {
962 6198230 : switch(typ(y))
963 : {
964 6197775 : 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 = gerepileupto(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 : 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;
1030 : GEN q, r;
1031 1057 : if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
1032 7 : vx = varn(x); if (vx != v) x = swap_vars(x,v);
1033 7 : vy = varn(y); if (vy != v) y = swap_vars(y,v);
1034 7 : q = RgX_divrem(x,y, &r);
1035 7 : if (v && (vx != v || vy != v))
1036 : {
1037 7 : GEN X = pol_x(v);
1038 7 : q = gsubst(q, v, X); /* poleval broken for t_RFRAC, subst is safe */
1039 7 : r = gsubst(r, v, X);
1040 : }
1041 7 : return gerepilecopy(av, mkcol2(q, r));
1042 : }
1043 :
1044 : GEN
1045 63195934 : diviiround(GEN x, GEN y)
1046 : {
1047 63195934 : pari_sp av1, av = avma;
1048 : GEN q,r;
1049 : int fl;
1050 :
1051 63195934 : q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
1052 63187772 : if (r==gen_0) return q;
1053 33497684 : av1 = avma;
1054 33497684 : fl = abscmpii(shifti(r,1),y);
1055 33501326 : set_avma(av1); cgiv(r);
1056 33510612 : if (fl >= 0) /* If 2*|r| >= |y| */
1057 : {
1058 18063651 : long sz = signe(x)*signe(y);
1059 18063651 : if (fl || sz > 0) q = gerepileuptoint(av, addis(q,sz));
1060 : }
1061 33512074 : return q;
1062 : }
1063 :
1064 : static GEN
1065 518 : _abs(GEN x)
1066 : {
1067 518 : if (typ(x) == t_QUAD) return (gsigne(x) < 0)? gneg(x): x;
1068 364 : return R_abs_shallow(x);
1069 : }
1070 :
1071 : /* If x and y are not both scalars, same as gdivent.
1072 : * Otherwise, compute the quotient x/y, rounded to the nearest integer
1073 : * (towards +oo in case of tie). */
1074 : GEN
1075 1465604 : gdivround(GEN x, GEN y)
1076 : {
1077 : pari_sp av;
1078 1465604 : long tx = typ(x), ty = typ(y);
1079 : GEN q, r;
1080 :
1081 1465604 : if (tx == t_INT && ty == t_INT) return diviiround(x,y);
1082 176147 : av = avma;
1083 176147 : if (is_realext(x) && is_realext(y))
1084 : { /* same as diviiround, less efficient */
1085 : pari_sp av1;
1086 : int fl;
1087 259 : q = quotrem(x,y,&r); av1 = avma;
1088 259 : fl = gcmp(gmul2n(_abs(r),1), _abs(y));
1089 259 : set_avma(av1); cgiv(r);
1090 259 : if (fl >= 0) /* If 2*|r| >= |y| */
1091 : {
1092 84 : long sz = gsigne(y);
1093 84 : if (fl || sz > 0) q = gerepileupto(av, gaddgs(q, sz));
1094 : }
1095 259 : return q;
1096 : }
1097 1587130 : if (is_matvec_t(tx)) pari_APPLY_same(gdivround(gel(x,i),y));
1098 931 : return gdivent(x,y);
1099 : }
1100 :
1101 : GEN
1102 0 : gdivmod(GEN x, GEN y, GEN *pr)
1103 : {
1104 0 : switch(typ(x))
1105 : {
1106 0 : case t_INT:
1107 0 : switch(typ(y))
1108 : {
1109 0 : case t_INT: return dvmdii(x,y,pr);
1110 0 : case t_POL: *pr=icopy(x); return gen_0;
1111 : }
1112 0 : break;
1113 0 : case t_POL: return poldivrem(x,y,pr);
1114 : }
1115 0 : pari_err_TYPE2("gdivmod",x,y);
1116 : return NULL;/*LCOV_EXCL_LINE*/
1117 : }
1118 :
1119 : /*******************************************************************/
1120 : /* */
1121 : /* SHIFT */
1122 : /* */
1123 : /*******************************************************************/
1124 :
1125 : /* Shift tronque si n<0 (multiplication tronquee par 2^n) */
1126 :
1127 : GEN
1128 47665529 : gshift(GEN x, long n)
1129 : {
1130 47665529 : switch(typ(x))
1131 : {
1132 38928071 : case t_INT: return shifti(x,n);
1133 8039709 : case t_REAL:return shiftr(x,n);
1134 2181619 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gshift(gel(x,i),n));
1135 : }
1136 141166 : return gmul2n(x,n);
1137 : }
1138 :
1139 : /*******************************************************************/
1140 : /* */
1141 : /* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
1142 : /* */
1143 : /*******************************************************************/
1144 :
1145 : /* Convert t_SER --> t_POL, ignoring valser. INTERNAL ! */
1146 : GEN
1147 10123592 : ser2pol_i(GEN x, long lx)
1148 : {
1149 10123592 : long i = lx-1;
1150 : GEN y;
1151 13981194 : while (i > 1 && isrationalzero(gel(x,i))) i--;
1152 10123592 : if (!signe(x))
1153 : { /* danger */
1154 119 : if (i == 1) return zeropol(varn(x));
1155 119 : y = cgetg(3,t_POL); y[1] = x[1] & ~VALSERBITS;
1156 119 : gel(y,2) = gel(x,2); return y;
1157 : }
1158 10123473 : y = cgetg(i+1, t_POL); y[1] = x[1] & ~VALSERBITS;
1159 42851236 : for ( ; i > 1; i--) gel(y,i) = gel(x,i);
1160 10123473 : return y;
1161 : }
1162 :
1163 : GEN
1164 756803 : ser2pol_i_normalize(GEN x, long l, long *v)
1165 : {
1166 756803 : long i = 2, j = l-1, k;
1167 : GEN y;
1168 756838 : while (i < l && gequal0(gel(x,i))) i++;
1169 756803 : *v = i - 2; if (i == l) return zeropol(varn(x));
1170 1001137 : while (j > i && gequal0(gel(x,j))) j--;
1171 756789 : l = j - *v + 1;
1172 756789 : y = cgetg(l, t_POL); y[1] = x[1] & ~VALSERBITS;
1173 3911849 : k = l; while (k > 2) gel(y, --k) = gel(x,j--);
1174 756789 : return y;
1175 : }
1176 :
1177 : GEN
1178 45542 : ser_inv(GEN b)
1179 : {
1180 45542 : pari_sp av = avma;
1181 45542 : long e, l = lg(b);
1182 : GEN x, y;
1183 45542 : y = ser2pol_i_normalize(b, l, &e);
1184 45542 : if (e)
1185 : {
1186 0 : pari_warn(warner,"normalizing a series with 0 leading term");
1187 0 : l -= e; if (l <= 2) pari_err_INV("inv_ser", b);
1188 : }
1189 45542 : y = RgXn_inv_i(y, l-2);
1190 45535 : x = RgX_to_ser(y, l); setvalser(x, - valser(b) - e);
1191 45535 : return gerepilecopy(av, x);
1192 : }
1193 :
1194 : /* T t_POL in var v, mod out by T components of x which are t_POL/t_RFRAC in v.
1195 : * Recursively. Make sure that resulting polynomials of degree 0 in v are
1196 : * simplified (map K[X]_0 to K) */
1197 : static GEN
1198 196 : mod_r(GEN x, long v, GEN T)
1199 : {
1200 196 : long w, tx = typ(x);
1201 : GEN y;
1202 :
1203 196 : if (is_const_t(tx)) return x;
1204 175 : switch(tx)
1205 : {
1206 7 : case t_POLMOD:
1207 7 : w = varn(gel(x,1));
1208 7 : if (w == v) pari_err_PRIORITY("subst", gel(x,1), "=", v);
1209 7 : if (varncmp(v, w) < 0) return x;
1210 7 : return gmodulo(mod_r(gel(x,2),v,T), mod_r(gel(x,1),v,T));
1211 7 : case t_SER:
1212 7 : w = varn(x);
1213 7 : if (w == v) break; /* fail */
1214 7 : if (varncmp(v, w) < 0 || ser_isexactzero(x)) return x;
1215 21 : pari_APPLY_ser(mod_r(gel(x,i),v,T));
1216 133 : case t_POL:
1217 133 : w = varn(x);
1218 133 : if (w == v)
1219 : {
1220 105 : x = RgX_rem(x, T);
1221 105 : if (!degpol(x)) x = gel(x,2);
1222 105 : return x;
1223 : }
1224 28 : if (varncmp(v, w) < 0) return x;
1225 98 : pari_APPLY_pol(mod_r(gel(x,i),v,T));
1226 14 : case t_RFRAC:
1227 14 : x = gdiv(mod_r(gel(x,1),v,T), mod_r(gel(x,2),v,T));
1228 14 : if (typ(x) == t_POL && varn(x) == v && lg(x) == 3) x = gel(x,2);
1229 14 : return x;
1230 7 : case t_VEC: case t_COL: case t_MAT:
1231 21 : pari_APPLY_same(mod_r(gel(x,i),v,T));
1232 7 : case t_LIST:
1233 7 : y = mklist();
1234 7 : list_data(y) = list_data(x)? mod_r(list_data(x),v,T): NULL;
1235 7 : return y;
1236 : }
1237 0 : pari_err_TYPE("substpol",x);
1238 : return NULL;/*LCOV_EXCL_LINE*/
1239 : }
1240 : GEN
1241 8708 : gsubstpol(GEN x, GEN T, GEN y)
1242 : {
1243 8708 : pari_sp av = avma;
1244 : long v;
1245 : GEN z;
1246 8708 : if (typ(T) == t_POL && RgX_is_monomial(T) && gequal1(leading_coeff(T)))
1247 : { /* T = t^d */
1248 8687 : long d = degpol(T);
1249 8687 : v = varn(T); z = (d==1)? x: gdeflate(x, v, d);
1250 8673 : if (z) return gerepileupto(av, gsubst(z, v, y));
1251 : }
1252 49 : v = fetch_var(); T = simplify_shallow(T);
1253 49 : if (typ(T) == t_RFRAC)
1254 21 : z = gsub(gel(T,1), gmul(pol_x(v), gel(T,2)));
1255 : else
1256 28 : z = gsub(T, pol_x(v));
1257 49 : z = mod_r(x, gvar(T), z);
1258 49 : z = gsubst(z, v, y); (void)delete_var();
1259 49 : return gerepileupto(av, z);
1260 : }
1261 :
1262 : long
1263 852711 : RgX_deflate_order(GEN x)
1264 : {
1265 852711 : ulong d = 0, i, lx = (ulong)lg(x);
1266 1769409 : for (i=3; i<lx; i++)
1267 1570501 : if (!gequal0(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1268 198908 : return d? (long)d: 1;
1269 : }
1270 : long
1271 512375 : ZX_deflate_order(GEN x)
1272 : {
1273 512375 : ulong d = 0, i, lx = (ulong)lg(x);
1274 1604961 : for (i=3; i<lx; i++)
1275 1413958 : if (signe(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1276 191003 : return d? (long)d: 1;
1277 : }
1278 :
1279 : /* deflate (non-leaf) x recursively */
1280 : static GEN
1281 63 : vdeflate(GEN x, long v, long d)
1282 : {
1283 63 : long i = lontyp[typ(x)], lx;
1284 63 : GEN z = cgetg_copy(x, &lx);
1285 63 : if (i == 2) z[1] = x[1];
1286 154 : for (; i<lx; i++)
1287 : {
1288 133 : gel(z,i) = gdeflate(gel(x,i),v,d);
1289 133 : if (!z[i]) return NULL;
1290 : }
1291 21 : return z;
1292 : }
1293 :
1294 : /* don't return NULL if substitution fails (fallback won't be able to handle
1295 : * t_SER anyway), fail with a meaningful message */
1296 : static GEN
1297 5768 : serdeflate(GEN x, long v, long d)
1298 : {
1299 5768 : long V, dy, lx, vx = varn(x);
1300 : pari_sp av;
1301 : GEN y;
1302 5768 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1303 5761 : if (varncmp(vx, v) > 0) return gcopy(x);
1304 5761 : av = avma;
1305 5761 : V = valser(x);
1306 5761 : lx = lg(x);
1307 5761 : if (lx == 2) return zeroser(v, V / d);
1308 5761 : y = ser2pol_i(x, lx);
1309 5761 : dy = degpol(y);
1310 5761 : if (V % d != 0 || (dy > 0 && RgX_deflate_order(y) % d != 0))
1311 : {
1312 14 : const char *s = stack_sprintf("valuation(x) %% %ld", d);
1313 14 : pari_err_DOMAIN("gdeflate", s, "!=", gen_0,x);
1314 : }
1315 5747 : if (dy > 0) y = RgX_deflate(y, d);
1316 5747 : y = RgX_to_ser(y, 3 + (lx-3)/d);
1317 5747 : setvalser(y, V/d); return gerepilecopy(av, y);
1318 : }
1319 : static GEN
1320 8722 : poldeflate(GEN x, long v, long d)
1321 : {
1322 8722 : long vx = varn(x);
1323 : pari_sp av;
1324 8722 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1325 8694 : if (varncmp(vx, v) > 0 || degpol(x) <= 0) return gcopy(x);
1326 8659 : av = avma;
1327 : /* x nonconstant */
1328 8659 : if (RgX_deflate_order(x) % d != 0) return NULL;
1329 8631 : return gerepilecopy(av, RgX_deflate(x,d));
1330 : }
1331 : static GEN
1332 21 : listdeflate(GEN x, long v, long d)
1333 : {
1334 21 : GEN y = NULL, z = mklist();
1335 21 : if (list_data(x))
1336 : {
1337 14 : y = vdeflate(list_data(x),v,d);
1338 14 : if (!y) return NULL;
1339 : }
1340 14 : list_data(z) = y; return z;
1341 : }
1342 : /* return NULL if substitution fails */
1343 : GEN
1344 14553 : gdeflate(GEN x, long v, long d)
1345 : {
1346 14553 : if (d <= 0) pari_err_DOMAIN("gdeflate", "degree", "<=", gen_0,stoi(d));
1347 14553 : switch(typ(x))
1348 : {
1349 28 : case t_INT:
1350 : case t_REAL:
1351 : case t_INTMOD:
1352 : case t_FRAC:
1353 : case t_FFELT:
1354 : case t_COMPLEX:
1355 : case t_PADIC:
1356 28 : case t_QUAD: return gcopy(x);
1357 8722 : case t_POL: return poldeflate(x,v,d);
1358 5768 : case t_SER: return serdeflate(x,v,d);
1359 7 : case t_POLMOD:
1360 7 : if (varncmp(varn(gel(x,1)), v) >= 0) return gcopy(x);
1361 : /* fall through */
1362 : case t_RFRAC:
1363 : case t_VEC:
1364 : case t_COL:
1365 14 : case t_MAT: return vdeflate(x,v,d);
1366 21 : case t_LIST: return listdeflate(x,v,d);
1367 : }
1368 0 : pari_err_TYPE("gdeflate",x);
1369 : return NULL; /* LCOV_EXCL_LINE */
1370 : }
1371 :
1372 : /* set *m to the largest d such that x0 = A(X^d); return A */
1373 : GEN
1374 530759 : RgX_deflate_max(GEN x, long *m)
1375 : {
1376 530759 : *m = RgX_deflate_order(x);
1377 530771 : return RgX_deflate(x, *m);
1378 : }
1379 : GEN
1380 316873 : ZX_deflate_max(GEN x, long *m)
1381 : {
1382 316873 : *m = ZX_deflate_order(x);
1383 316872 : return RgX_deflate(x, *m);
1384 : }
1385 :
1386 : static int
1387 21161 : serequalXk(GEN x)
1388 : {
1389 21161 : long i, l = lg(x);
1390 21161 : if (l == 2 || !isint1(gel(x,2))) return 0;
1391 9786 : for (i = 3; i < l; i++)
1392 7833 : if (!isintzero(gel(x,i))) return 0;
1393 1953 : return 1;
1394 : }
1395 :
1396 : static GEN
1397 84 : gsubst_v(GEN e, long v, GEN x)
1398 245 : { pari_APPLY_same(gsubst(e, v, gel(x,i))); }
1399 :
1400 : static GEN
1401 14 : constmat(GEN z, long n)
1402 : {
1403 14 : GEN y = cgetg(n+1, t_MAT), c = const_col(n, gcopy(z));
1404 : long i;
1405 35 : for (i = 1; i <= n; i++) gel(y, i) = c;
1406 14 : return y;
1407 : }
1408 : static GEN
1409 56 : scalarmat2(GEN o, GEN z, long n)
1410 : {
1411 : GEN y;
1412 : long i;
1413 56 : if (n == 0) return cgetg(1, t_MAT);
1414 56 : if (n == 1) retmkmat(mkcol(gcopy(o)));
1415 35 : y = cgetg(n+1, t_MAT); z = gcopy(z); o = gcopy(o);
1416 105 : for (i = 1; i <= n; i++) { gel(y, i) = const_col(n, z); gcoeff(y,i,i) = o; }
1417 35 : return y;
1418 : }
1419 : /* x * y^0, n = dim(y) if t_MAT, else -1 */
1420 : static GEN
1421 751163 : subst_higher(GEN x, GEN y, long n)
1422 : {
1423 751163 : GEN o = Rg_get_1(y);
1424 751163 : if (o == gen_1) return n < 0? gcopy(x): scalarmat(x,n);
1425 63 : x = gmul(x,o); return n < 0? x: scalarmat2(x, Rg_get_0(y), n);
1426 : }
1427 :
1428 : /* x t_POLMOD, v strictly lower priority than var(x) */
1429 : static GEN
1430 10969 : subst_polmod(GEN x, long v, GEN y)
1431 : {
1432 : long l, i;
1433 10969 : GEN a = gsubst(gel(x,2),v,y), b = gsubst(gel(x,1),v,y), z;
1434 :
1435 10969 : if (typ(b) != t_POL) pari_err_TYPE2("substitution",x,y);
1436 10969 : if (typ(a) != t_POL || varncmp(varn(a), varn(b)) >= 0) return gmodulo(a, b);
1437 532 : l = lg(a); z = cgetg(l,t_POL); z[1] = a[1];
1438 4067 : for (i = 2; i < l; i++) gel(z,i) = gmodulo(gel(a,i),b);
1439 532 : return normalizepol_lg(z, l);
1440 : }
1441 : /* Trunc to n terms; x + O(t^(n + v(x))). FIXME: export ? */
1442 : static GEN
1443 70 : sertrunc(GEN x, long n)
1444 : {
1445 70 : long i, l = n + 2;
1446 : GEN y;
1447 70 : if (l >= lg(x)) return x;
1448 14 : if (n <= 0) return zeroser(varn(x), n + valser(x));
1449 14 : y = cgetg(l, t_SER);
1450 28 : for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
1451 14 : y[1] = x[1]; return y;
1452 : }
1453 : /* FIXME: export ? */
1454 : static GEN
1455 1960 : sertrunc_copy(GEN x, long n)
1456 : {
1457 1960 : long i, l = minss(n + 2, lg(x));
1458 1960 : GEN y = cgetg(l, t_SER);
1459 13349 : for (i = 2; i < l; i++) gel(y,i) = gcopy(gel(x,i));
1460 1960 : y[1] = x[1]; return y;
1461 : }
1462 :
1463 : GEN
1464 2726355 : gsubst(GEN x, long v, GEN y)
1465 : {
1466 2726355 : long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
1467 : long l, vx, vy, ex, ey, i, j, k, jb, matn;
1468 : pari_sp av, av2;
1469 : GEN X, t, z;
1470 :
1471 2726355 : switch(ty)
1472 : {
1473 84 : case t_VEC: case t_COL:
1474 84 : return gsubst_v(x, v, y);
1475 175 : case t_MAT:
1476 175 : if (ly==1) return cgetg(1,t_MAT);
1477 168 : if (ly == lgcols(y)) { matn = ly - 1; break; }
1478 : /* fall through */
1479 : case t_QFB:
1480 7 : pari_err_TYPE2("substitution",x,y);
1481 2726096 : default: matn = -1;
1482 : }
1483 2726257 : if (is_scalar_t(tx))
1484 : {
1485 413760 : if (tx == t_POLMOD && varncmp(v, varn(gel(x,1))) > 0)
1486 : {
1487 10969 : av = avma;
1488 10969 : return gerepileupto(av, subst_polmod(x, v, y));
1489 : }
1490 402791 : return subst_higher(x, y, matn);
1491 : }
1492 :
1493 2312497 : switch(tx)
1494 : {
1495 2072175 : case t_POL:
1496 2072175 : vx = varn(x);
1497 2072175 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1498 2070075 : if (varncmp(vx, v) < 0)
1499 : {
1500 177143 : av = avma; z = cgetg(lx, t_POL); z[1] = x[1];
1501 177143 : if (lx == 2) return z;
1502 823549 : for (i = 2; i < lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
1503 176702 : z = normalizepol_lg(z, lx); lx = lg(z);
1504 176702 : if (lx == 2) { set_avma(av); return zeropol(vx); }
1505 176688 : if (lx == 3) return gerepileupto(av, gmul(pol_1(vx), gel(z,2)));
1506 154449 : return gerepileupto(av, poleval(z, pol_x(vx)));
1507 : }
1508 : /* v = vx */
1509 1892932 : if (lx == 2)
1510 : {
1511 27867 : GEN z = Rg_get_0(y);
1512 27867 : return matn >= 0? constmat(z, matn): z;
1513 : }
1514 1865065 : if (lx == 3)
1515 : {
1516 346265 : x = subst_higher(gel(x,2), y, matn);
1517 346265 : if (matn >= 0) return x;
1518 346251 : vy = gvar(y);
1519 346251 : return (vy == NO_VARIABLE)? x: gmul(x, pol_1(vy));
1520 : }
1521 1518800 : return matn >= 0? RgX_RgM_eval(x, y): poleval(x,y);
1522 :
1523 26621 : case t_SER:
1524 26621 : vx = varn(x);
1525 26621 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1526 26614 : ex = valser(x);
1527 26614 : if (varncmp(vx, v) < 0)
1528 : {
1529 56 : if (lx == 2) return matn >= 0? scalarmat(x, matn): gcopy(x);
1530 56 : av = avma; X = pol_x(vx);
1531 56 : av2 = avma;
1532 56 : z = gadd(gsubst(gel(x,lx-1),v,y), zeroser(vx,1));
1533 224 : for (i = lx-2; i>=2; i--)
1534 : {
1535 168 : z = gadd(gmul(z,X), gsubst(gel(x,i),v,y));
1536 168 : if (gc_needed(av2,1))
1537 : {
1538 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1539 0 : z = gerepileupto(av2, z);
1540 : }
1541 : }
1542 56 : if (ex) z = gmul(z, pol_xnall(ex,vx));
1543 56 : return gerepileupto(av, z);
1544 : }
1545 26558 : switch(ty) /* here vx == v */
1546 : {
1547 21273 : case t_SER:
1548 21273 : vy = varn(y); ey = valser(y);
1549 21273 : if (ey < 1 || lx == 2) return zeroser(vy, ey*(ex+lx-2));
1550 21273 : if (ey == 1 && serequalXk(y)
1551 1953 : && (varncmp(vx,vy) >= 0 || varncmp(gvar2(x), vy) >= 0))
1552 : { /* y = t + O(t^N) */
1553 1953 : if (lx > ly)
1554 : { /* correct number of significant terms */
1555 1624 : l = ly;
1556 1624 : if (!ex)
1557 1603 : for (i = 3; i < lx; i++)
1558 1603 : if (++l >= lx || !gequal0(gel(x,i))) break;
1559 1624 : lx = l;
1560 : }
1561 1953 : z = sertrunc_copy(x, lx - 2); if (vx != vy) setvarn(z,vy);
1562 1953 : return z;
1563 : }
1564 19320 : if (vy != vx)
1565 : {
1566 28 : long nx = lx - 2, n = minss(ey * nx, ly - 2);
1567 28 : av = avma; z = gel(x, nx+1);
1568 91 : for (i = nx; i > 1; i--)
1569 : {
1570 63 : z = gadd(gmul(y,z), gel(x,i));
1571 63 : if (gc_needed(av,1))
1572 : {
1573 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1574 0 : z = gerepileupto(av, z);
1575 : }
1576 : }
1577 28 : if (ex)
1578 : {
1579 21 : if (ex < 0) { y = ginv(y); ex = -ex; }
1580 21 : z = gmul(z, gpowgs(sertrunc(y, n), ex));
1581 : }
1582 28 : if (lg(z)-2 > n) z = sertrunc_copy(z, n);
1583 28 : return gerepileupto(av,z);
1584 : }
1585 19292 : l = (lx-2)*ey + 2;
1586 19292 : if (ex) { if (l>ly) l = ly; }
1587 19243 : else if (lx != 3)
1588 : {
1589 19257 : for (i = 3; i < lx; i++)
1590 19257 : if (!isexactzero(gel(x,i))) break;
1591 19243 : l = minss(l, (i-2)*ey + (gequal0(y)? 2 : ly));
1592 : }
1593 19292 : av = avma; t = leafcopy(y);
1594 19292 : if (l < ly) setlg(t, l);
1595 19292 : z = scalarser(gen_1, varn(y), l-2);
1596 19292 : gel(z,2) = gel(x,2); /* ensure lg(z) = l even if x[2] = 0 */
1597 77224 : for (i = 3, jb = ey; jb <= l-2; i++,jb += ey)
1598 : {
1599 57939 : if (i < lx) {
1600 132286 : for (j = jb+2; j < minss(l, jb+ly); j++)
1601 74424 : gel(z,j) = gadd(gel(z,j), gmul(gel(x,i),gel(t,j-jb)));
1602 : }
1603 93275 : for (j = minss(ly-1, l-1-jb-ey); j > 1; j--)
1604 : {
1605 35343 : GEN a = gmul(gel(t,2), gel(y,j));
1606 84483 : for (k=2; k<j; k++) a = gadd(a, gmul(gel(t,j-k+2), gel(y,k)));
1607 35343 : gel(t,j) = a;
1608 : }
1609 57932 : if (gc_needed(av,1))
1610 : {
1611 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst");
1612 0 : gerepileall(av,2, &z,&t);
1613 : }
1614 : }
1615 19285 : if (!ex) return gerepilecopy(av,z);
1616 49 : if (ex < 0) { ex = -ex; y = ginv(y); }
1617 49 : return gerepileupto(av, gmul(z, gpowgs(sertrunc(y, l-2), ex)));
1618 :
1619 5243 : case t_POL: case t_RFRAC:
1620 : {
1621 5243 : long N, n = lx-2;
1622 5243 : vy = gvar(y); ey = gval(y,vy);
1623 5243 : if (ey == LONG_MAX)
1624 : { /* y = 0 */
1625 49 : if (ex < 0) pari_err_INV("gsubst",y);
1626 35 : if (!n) return gcopy(x);
1627 28 : if (ex > 0) return Rg_get_0(ty == t_RFRAC? gel(y,2): y);
1628 14 : y = Rg_get_1(ty == t_RFRAC? gel(y,2): y);
1629 14 : return gmul(y, gel(x,2));
1630 : }
1631 5194 : if (ey < 1 || n == 0) return zeroser(vy, ey*(ex+n));
1632 5187 : av = avma;
1633 5187 : n *= ey;
1634 5187 : N = ex? n: maxss(n-ey,1);
1635 5187 : y = (ty == t_RFRAC)? rfrac_to_ser_i(y, N+2): RgX_to_ser(y, N+2);
1636 5187 : if (lg(y)-2 > n) setlg(y, n+2);
1637 5187 : x = ser2pol_i(x, lx);
1638 5187 : if (varncmp(vy,vx) > 0)
1639 49 : z = gadd(poleval(x, y), zeroser(vy,n));
1640 : else
1641 : {
1642 5138 : z = RgXn_eval(x, ser2rfrac_i(y), n);
1643 5138 : if (varn(z) == vy) z = RgX_to_ser(z, n+2);
1644 0 : else z = scalarser(z, vy, n);
1645 : }
1646 5187 : if (!ex) return gerepilecopy(av, z);
1647 5082 : return gerepileupto(av, gmul(z, gpowgs(y,ex)));
1648 : }
1649 :
1650 42 : default:
1651 42 : if (isexactzero(y))
1652 : {
1653 35 : if (ex < 0) pari_err_INV("gsubst",y);
1654 14 : if (ex > 0) return gcopy(y);
1655 7 : if (lx > 2) return gadd(gel(x,2), y); /*add maps to correct ring*/
1656 : }
1657 7 : pari_err_TYPE2("substitution",x,y);
1658 : }
1659 0 : break;
1660 :
1661 1253 : case t_RFRAC:
1662 : {
1663 1253 : GEN a = gel(x,1), b = gel(x,2);
1664 1253 : av = avma;
1665 1253 : a = gsubst(a, v, y);
1666 1253 : b = gsubst(b, v, y); return gerepileupto(av, gdiv(a, b));
1667 : }
1668 :
1669 660101 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gsubst(gel(x,i),v,y));
1670 56 : case t_LIST:
1671 56 : z = mklist();
1672 56 : list_data(z) = list_data(x)? gsubst(list_data(x),v,y): NULL;
1673 56 : return z;
1674 : }
1675 0 : return gcopy(x);
1676 : }
1677 :
1678 : /* Return P(x * h), not memory clean */
1679 : GEN
1680 4193 : ser_unscale(GEN P, GEN h)
1681 : {
1682 4193 : long l = lg(P);
1683 4193 : GEN Q = cgetg(l,t_SER);
1684 4193 : Q[1] = P[1];
1685 4193 : if (l != 2)
1686 : {
1687 4193 : long i = 2;
1688 4193 : GEN hi = gpowgs(h, valser(P));
1689 4193 : gel(Q,i) = gmul(gel(P,i), hi);
1690 200508 : for (i++; i<l; i++) { hi = gmul(hi,h); gel(Q,i) = gmul(gel(P,i), hi); }
1691 : }
1692 4193 : return Q;
1693 : }
1694 :
1695 : GEN
1696 959 : gsubstvec(GEN e, GEN v, GEN r)
1697 : {
1698 959 : pari_sp av = avma;
1699 959 : long i, j, k, l = lg(v);
1700 : GEN w, z, R;
1701 959 : if ( !is_vec_t(typ(v)) ) pari_err_TYPE("substvec",v);
1702 959 : if ( !is_vec_t(typ(r)) ) pari_err_TYPE("substvec",r);
1703 959 : if (lg(r)!=l) pari_err_DIM("substvec");
1704 959 : w = cgetg(l, t_VECSMALL);
1705 959 : z = cgetg(l, t_VECSMALL);
1706 959 : R = cgetg(l, t_VEC); k = 0;
1707 4235 : for(i = j = 1; i < l; i++)
1708 : {
1709 3276 : GEN T = gel(v,i), ri = gel(r,i);
1710 3276 : if (!gequalX(T)) pari_err_TYPE("substvec [not a variable]", T);
1711 3276 : if (gvar(ri) == NO_VARIABLE) /* no need to take precautions */
1712 : {
1713 1855 : e = gsubst(e, varn(T), ri);
1714 1855 : if (is_vec_t(typ(ri)) && k++) e = shallowconcat1(e);
1715 : }
1716 : else
1717 : {
1718 1421 : w[j] = varn(T);
1719 1421 : z[j] = fetch_var();
1720 1421 : gel(R,j) = ri; j++;
1721 : }
1722 : }
1723 2380 : for(i = 1; i < j; i++) e = gsubst(e,w[i],pol_x(z[i]));
1724 2380 : for(i = 1; i < j; i++)
1725 : {
1726 1421 : e = gsubst(e,z[i],gel(R,i));
1727 1421 : if (is_vec_t(typ(gel(R,i))) && k++) e = shallowconcat1(e);
1728 : }
1729 2380 : for(i = 1; i < j; i++) (void)delete_var();
1730 959 : return k > 1? gerepilecopy(av, e): gerepileupto(av, e);
1731 : }
1732 :
1733 : /*******************************************************************/
1734 : /* */
1735 : /* SERIE RECIPROQUE D'UNE SERIE */
1736 : /* */
1737 : /*******************************************************************/
1738 :
1739 : GEN
1740 98 : serreverse(GEN x)
1741 : {
1742 98 : long v=varn(x), lx = lg(x), i, mi;
1743 98 : pari_sp av0 = avma, av;
1744 : GEN a, y, u;
1745 :
1746 98 : if (typ(x)!=t_SER) pari_err_TYPE("serreverse",x);
1747 98 : if (valser(x)!=1) pari_err_DOMAIN("serreverse", "valuation", "!=", gen_1,x);
1748 91 : if (lx < 3) pari_err_DOMAIN("serreverse", "x", "=", gen_0,x);
1749 91 : y = ser_normalize(x);
1750 91 : if (y == x) a = NULL; else { a = gel(x,2); x = y; }
1751 91 : av = avma;
1752 252 : mi = lx-1; while (mi>=3 && gequal0(gel(x,mi))) mi--;
1753 91 : u = cgetg(lx,t_SER);
1754 91 : y = cgetg(lx,t_SER);
1755 91 : u[1] = y[1] = evalsigne(1) | _evalvalser(1) | evalvarn(v);
1756 91 : gel(u,2) = gel(y,2) = gen_1;
1757 91 : if (lx > 3)
1758 : {
1759 84 : gel(u,3) = gmulsg(-2,gel(x,3));
1760 84 : gel(y,3) = gneg(gel(x,3));
1761 : }
1762 1113 : for (i=3; i<lx-1; )
1763 : {
1764 : pari_sp av2;
1765 : GEN p1;
1766 1022 : long j, k, K = minss(i,mi);
1767 8456 : for (j=3; j<i+1; j++)
1768 : {
1769 7434 : av2 = avma; p1 = gel(x,j);
1770 39291 : for (k = maxss(3,j+2-mi); k < j; k++)
1771 31857 : p1 = gadd(p1, gmul(gel(u,k),gel(x,j-k+2)));
1772 7434 : p1 = gneg(p1);
1773 7434 : gel(u,j) = gerepileupto(av2, gadd(gel(u,j), p1));
1774 : }
1775 1022 : av2 = avma;
1776 1022 : p1 = gmulsg(i,gel(x,i+1));
1777 8309 : for (k = 2; k < K; k++)
1778 : {
1779 7287 : GEN p2 = gmul(gel(x,k+1),gel(u,i-k+2));
1780 7287 : p1 = gadd(p1, gmulsg(k,p2));
1781 : }
1782 1022 : i++;
1783 1022 : gel(u,i) = gerepileupto(av2, gneg(p1));
1784 1022 : gel(y,i) = gdivgu(gel(u,i), i-1);
1785 1022 : if (gc_needed(av,2))
1786 : {
1787 0 : GEN dummy = cgetg(1,t_VEC);
1788 0 : if(DEBUGMEM>1) pari_warn(warnmem,"serreverse");
1789 0 : for(k=i+1; k<lx; k++) gel(u,k) = gel(y,k) = dummy;
1790 0 : gerepileall(av,2, &u,&y);
1791 : }
1792 : }
1793 91 : if (a) y = ser_unscale(y, ginv(a));
1794 91 : return gerepilecopy(av0,y);
1795 : }
1796 :
1797 : /*******************************************************************/
1798 : /* */
1799 : /* DERIVATION ET INTEGRATION */
1800 : /* */
1801 : /*******************************************************************/
1802 : GEN
1803 25410 : derivser(GEN x)
1804 : {
1805 25410 : long i, vx = varn(x), e = valser(x), lx = lg(x);
1806 : GEN y;
1807 25410 : if (ser_isexactzero(x))
1808 : {
1809 7 : x = gcopy(x);
1810 7 : if (e) setvalser(x,e-1);
1811 7 : return x;
1812 : }
1813 25403 : if (e)
1814 : {
1815 602 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|evalvalser(e-1) | evalvarn(vx);
1816 22960 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i+e-2,gel(x,i));
1817 : } else {
1818 24801 : if (lx == 3) return zeroser(vx, 0);
1819 20951 : lx--;
1820 20951 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|_evalvalser(0) | evalvarn(vx);
1821 67382 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1822 : }
1823 21553 : return normalizeser(y);
1824 : }
1825 :
1826 : static GEN
1827 56 : rfrac_deriv(GEN x, long v)
1828 : {
1829 56 : pari_sp av = avma;
1830 56 : GEN y = cgetg(3,t_RFRAC), a = gel(x,1), b = gel(x,2), bp, b0, t, T;
1831 56 : long vx = varn(b);
1832 :
1833 56 : bp = deriv(b, v);
1834 56 : t = simplify_shallow(RgX_gcd(bp, b));
1835 56 : if (typ(t) != t_POL || varn(t) != vx)
1836 : {
1837 35 : if (gequal1(t)) b0 = b;
1838 : else
1839 : {
1840 0 : b0 = RgX_Rg_div(b, t);
1841 0 : bp = RgX_Rg_div(bp, t);
1842 : }
1843 35 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1844 35 : if (isexactzero(a)) return gerepileupto(av, a);
1845 35 : if (b0 == b)
1846 : {
1847 35 : gel(y,1) = gerepileupto((pari_sp)y, a);
1848 35 : gel(y,2) = RgX_sqr(b);
1849 : }
1850 : else
1851 : {
1852 0 : gel(y,1) = a;
1853 0 : gel(y,2) = RgX_Rg_mul(RgX_sqr(b0), t);
1854 0 : y = gerepilecopy(av, y);
1855 : }
1856 35 : return y;
1857 : }
1858 21 : b0 = gdivexact(b, t);
1859 21 : bp = gdivexact(bp,t);
1860 21 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1861 21 : if (isexactzero(a)) return gerepileupto(av, a);
1862 14 : T = RgX_gcd(a, t);
1863 14 : if (typ(T) != t_POL || varn(T) != vx)
1864 : {
1865 0 : a = gdiv(a, T);
1866 0 : t = gdiv(t, T);
1867 : }
1868 14 : else if (!gequal1(T))
1869 : {
1870 0 : a = gdivexact(a, T);
1871 0 : t = gdivexact(t, T);
1872 : }
1873 14 : gel(y,1) = a;
1874 14 : gel(y,2) = gmul(RgX_sqr(b0), t);
1875 14 : return gerepilecopy(av, y);
1876 : }
1877 :
1878 : GEN
1879 114128 : deriv(GEN x, long v)
1880 : {
1881 114128 : long tx = typ(x);
1882 114128 : if (is_const_t(tx))
1883 39795 : switch(tx)
1884 : {
1885 14 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
1886 14 : case t_FFELT: return FF_zero(x);
1887 39767 : default: return gen_0;
1888 : }
1889 74333 : if (v < 0)
1890 : {
1891 49 : if (tx == t_CLOSURE) return closure_deriv(x);
1892 49 : v = gvar9(x);
1893 : }
1894 74333 : switch(tx)
1895 : {
1896 14 : case t_POLMOD:
1897 : {
1898 14 : GEN a = gel(x,2), b = gel(x,1);
1899 14 : if (v == varn(b)) return Rg_get_0(b);
1900 7 : retmkpolmod(deriv(a,v), RgX_copy(b));
1901 : }
1902 74074 : case t_POL:
1903 74074 : switch(varncmp(varn(x), v))
1904 : {
1905 0 : case 1: return Rg_get_0(x);
1906 66479 : case 0: return RgX_deriv(x);
1907 : }
1908 113505 : pari_APPLY_pol(deriv(gel(x,i),v));
1909 :
1910 147 : case t_SER:
1911 147 : switch(varncmp(varn(x), v))
1912 : {
1913 0 : case 1: return Rg_get_0(x);
1914 133 : case 0: return derivser(x);
1915 : }
1916 14 : if (ser_isexactzero(x)) return gcopy(x);
1917 28 : pari_APPLY_ser(deriv(gel(x,i),v));
1918 :
1919 56 : case t_RFRAC:
1920 56 : return rfrac_deriv(x,v);
1921 :
1922 42 : case t_VEC: case t_COL: case t_MAT:
1923 84 : pari_APPLY_same(deriv(gel(x,i),v));
1924 : }
1925 0 : pari_err_TYPE("deriv",x);
1926 : return NULL; /* LCOV_EXCL_LINE */
1927 : }
1928 :
1929 : /* n-th derivative of t_SER x, n > 0 */
1930 : static GEN
1931 189 : derivnser(GEN x, long n)
1932 : {
1933 189 : long i, vx = varn(x), e = valser(x), lx = lg(x);
1934 : GEN y;
1935 189 : if (ser_isexactzero(x))
1936 : {
1937 7 : x = gcopy(x);
1938 7 : if (e) setvalser(x,e-n);
1939 7 : return x;
1940 : }
1941 182 : if (e < 0 || e >= n)
1942 : {
1943 154 : y = cgetg(lx,t_SER);
1944 154 : y[1] = evalsigne(1)| evalvalser(e-n) | evalvarn(vx);
1945 714 : for (i=0; i<lx-2; i++)
1946 560 : gel(y,i+2) = gmul(muls_interval(i+e-n+1,i+e), gel(x,i+2));
1947 : } else {
1948 28 : if (lx <= n+2) return zeroser(vx, 0);
1949 28 : lx -= n;
1950 28 : y = cgetg(lx,t_SER);
1951 28 : y[1] = evalsigne(1)|_evalvalser(0) | evalvarn(vx);
1952 91 : for (i=0; i<lx-2; i++)
1953 63 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n-e));
1954 : }
1955 182 : return normalizeser(y);
1956 : }
1957 :
1958 : /* n-th derivative of t_POL x, n > 0 */
1959 : static GEN
1960 833 : RgX_derivn(GEN x, long n)
1961 : {
1962 833 : long i, vx = varn(x), lx = lg(x);
1963 : GEN y;
1964 833 : if (lx <= n+2) return pol_0(vx);
1965 749 : lx -= n;
1966 749 : y = cgetg(lx,t_POL);
1967 749 : y[1] = evalsigne(1)| evalvarn(vx);
1968 50904 : for (i=0; i<lx-2; i++)
1969 50155 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n));
1970 749 : return normalizepol_lg(y, lx);
1971 : }
1972 :
1973 : static GEN
1974 42 : rfrac_derivn(GEN x, long n, long vs)
1975 : {
1976 42 : pari_sp av = avma;
1977 42 : GEN u = gel(x,1), v = gel(x,2);
1978 42 : GEN dv = deriv(v, vs);
1979 : long i;
1980 112 : for (i=1; i<=n; i++)
1981 : {
1982 70 : GEN du = deriv(u, vs);
1983 70 : u = gadd(gmul(du, v), gmulsg (-i, gmul(dv, u)));
1984 : }
1985 42 : v = gpowgs(v, n+1);
1986 42 : return gerepileupto(av, gdiv(u, v));
1987 : }
1988 :
1989 : GEN
1990 1351 : derivn(GEN x, long n, long v)
1991 : {
1992 : long tx;
1993 1351 : if (n < 0) pari_err_DOMAIN("derivn","n","<", gen_0, stoi(n));
1994 1344 : if (n == 0) return gcopy(x);
1995 1344 : tx = typ(x);
1996 1344 : if (is_const_t(tx))
1997 49 : switch(tx)
1998 : {
1999 21 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
2000 21 : case t_FFELT: return FF_zero(x);
2001 7 : default: return gen_0;
2002 : }
2003 1295 : if (v < 0)
2004 : {
2005 1057 : if (tx == t_CLOSURE) return closure_derivn(x, n);
2006 952 : v = gvar9(x);
2007 : }
2008 1190 : switch(tx)
2009 : {
2010 21 : case t_POLMOD:
2011 : {
2012 21 : GEN a = gel(x,2), b = gel(x,1);
2013 21 : if (v == varn(b)) return Rg_get_0(b);
2014 14 : retmkpolmod(derivn(a,n,v), RgX_copy(b));
2015 : }
2016 861 : case t_POL:
2017 861 : switch(varncmp(varn(x), v))
2018 : {
2019 0 : case 1: return Rg_get_0(x);
2020 833 : case 0: return RgX_derivn(x,n);
2021 : }
2022 84 : pari_APPLY_pol(derivn(gel(x,i),n,v));
2023 :
2024 196 : case t_SER:
2025 196 : switch(varncmp(varn(x), v))
2026 : {
2027 0 : case 1: return Rg_get_0(x);
2028 189 : case 0: return derivnser(x, n);
2029 : }
2030 7 : if (ser_isexactzero(x)) return gcopy(x);
2031 28 : pari_APPLY_ser(derivn(gel(x,i),n,v));
2032 :
2033 42 : case t_RFRAC:
2034 42 : return rfrac_derivn(x, n, v);
2035 :
2036 63 : case t_VEC: case t_COL: case t_MAT:
2037 126 : pari_APPLY_same(derivn(gel(x,i),n,v));
2038 : }
2039 7 : pari_err_TYPE("derivn",x);
2040 : return NULL; /* LCOV_EXCL_LINE */
2041 : }
2042 :
2043 : static long
2044 833 : lookup(GEN v, long vx)
2045 : {
2046 833 : long i ,l = lg(v);
2047 1491 : for(i=1; i<l; i++)
2048 1253 : if (varn(gel(v,i)) == vx) return i;
2049 238 : return 0;
2050 : }
2051 :
2052 : GEN
2053 3535 : diffop(GEN x, GEN v, GEN dv)
2054 : {
2055 : pari_sp av;
2056 3535 : long i, idx, lx, tx = typ(x), vx;
2057 : GEN y;
2058 3535 : if (!is_vec_t(typ(v))) pari_err_TYPE("diffop",v);
2059 3535 : if (!is_vec_t(typ(dv))) pari_err_TYPE("diffop",dv);
2060 3535 : if (lg(v)!=lg(dv)) pari_err_DIM("diffop");
2061 3535 : if (is_const_t(tx)) return gen_0;
2062 1148 : switch(tx)
2063 : {
2064 84 : case t_POLMOD:
2065 84 : av = avma;
2066 84 : vx = varn(gel(x,1)); idx = lookup(v,vx);
2067 84 : if (idx) /*Assume the users now what they are doing */
2068 0 : y = gmodulo(diffop(gel(x,2),v,dv), gel(x,1));
2069 : else
2070 : {
2071 84 : GEN m = gel(x,1), pol=gel(x,2);
2072 84 : GEN u = gneg(gdiv(diffop(m,v,dv),RgX_deriv(m)));
2073 84 : y = diffop(pol,v,dv);
2074 84 : if (typ(pol)==t_POL && varn(pol)==varn(m))
2075 70 : y = gadd(y, gmul(u,RgX_deriv(pol)));
2076 84 : y = gmodulo(y, gel(x,1));
2077 : }
2078 84 : return gerepileupto(av, y);
2079 952 : case t_POL:
2080 952 : if (signe(x)==0) return gen_0;
2081 742 : vx = varn(x); idx = lookup(v,vx);
2082 742 : av = avma; lx = lg(x);
2083 742 : y = diffop(gel(x,lx-1),v,dv);
2084 2842 : for (i=lx-2; i>=2; i--) y = gadd(gmul(y,pol_x(vx)),diffop(gel(x,i),v,dv));
2085 742 : if (idx) y = gadd(y, gmul(gel(dv,idx),RgX_deriv(x)));
2086 742 : return gerepileupto(av, y);
2087 :
2088 7 : case t_SER:
2089 7 : if (signe(x)==0) return gen_0;
2090 7 : vx = varn(x); idx = lookup(v,vx);
2091 7 : if (!idx) return gen_0;
2092 7 : av = avma;
2093 7 : if (ser_isexactzero(x)) y = x;
2094 : else
2095 : {
2096 7 : y = cgetg_copy(x, &lx); y[1] = x[1];
2097 119 : for (i=2; i<lx; i++) gel(y,i) = diffop(gel(x,i),v,dv);
2098 7 : y = normalizeser(y); /* y is probably invalid */
2099 7 : y = gsubst(y, vx, pol_x(vx)); /* Fix that */
2100 : }
2101 7 : y = gadd(y, gmul(gel(dv,idx),derivser(x)));
2102 7 : return gerepileupto(av, y);
2103 :
2104 105 : case t_RFRAC: {
2105 105 : GEN a = gel(x,1), b = gel(x,2), ap, bp;
2106 105 : av = avma;
2107 105 : ap = diffop(a, v, dv); bp = diffop(b, v, dv);
2108 105 : y = gsub(gdiv(ap,b),gdiv(gmul(a,bp),gsqr(b)));
2109 105 : return gerepileupto(av, y);
2110 : }
2111 :
2112 0 : case t_VEC: case t_COL: case t_MAT:
2113 0 : pari_APPLY_same(diffop(gel(x,i),v,dv));
2114 : }
2115 0 : pari_err_TYPE("diffop",x);
2116 : return NULL; /* LCOV_EXCL_LINE */
2117 : }
2118 :
2119 : GEN
2120 42 : diffop0(GEN x, GEN v, GEN dv, long n)
2121 : {
2122 42 : pari_sp av=avma;
2123 : long i;
2124 245 : for(i=1; i<=n; i++)
2125 203 : x = gerepileupto(av, diffop(x,v,dv));
2126 42 : return x;
2127 : }
2128 :
2129 : /********************************************************************/
2130 : /** **/
2131 : /** TAYLOR SERIES **/
2132 : /** **/
2133 : /********************************************************************/
2134 : /* swap vars (vx,v) in x (assume vx < v, vx main variable in x), then call
2135 : * act(data, v, x). FIXME: use in other places */
2136 : static GEN
2137 28 : swapvar_act(GEN x, long vx, long v, GEN (*act)(void*, long, GEN), void *data)
2138 : {
2139 28 : long v0 = fetch_var();
2140 28 : GEN y = act(data, v, gsubst(x,vx,pol_x(v0)));
2141 21 : y = gsubst(y,v0,pol_x(vx));
2142 21 : (void)delete_var(); return y;
2143 : }
2144 : /* x + O(v^data) */
2145 : static GEN
2146 14 : tayl_act(void *data, long v, GEN x) { return gadd(zeroser(v, (long)data), x); }
2147 : static GEN
2148 14 : integ_act(void *data, long v, GEN x) { (void)data; return integ(x,v); }
2149 :
2150 : GEN
2151 21 : tayl(GEN x, long v, long precS)
2152 : {
2153 21 : long vx = gvar9(x);
2154 : pari_sp av;
2155 :
2156 21 : if (varncmp(v, vx) <= 0) return gadd(zeroser(v,precS), x);
2157 14 : av = avma;
2158 14 : return gerepileupto(av, swapvar_act(x, vx, v, tayl_act, (void*)precS));
2159 : }
2160 :
2161 : GEN
2162 6979 : ggrando(GEN x, long n)
2163 : {
2164 : long m, v;
2165 :
2166 6979 : switch(typ(x))
2167 : {
2168 4032 : case t_INT:/* bug 3 + O(1) */
2169 4032 : if (signe(x) <= 0) pari_err_DOMAIN("O", "x", "<=", gen_0, x);
2170 4032 : if (!is_pm1(x)) return zeropadic(x,n);
2171 : /* +/-1 = x^0 */
2172 91 : v = m = 0; break;
2173 2940 : case t_POL:
2174 2940 : if (!signe(x)) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2175 2940 : v = varn(x);
2176 2940 : m = n * RgX_val(x); break;
2177 7 : case t_RFRAC:
2178 7 : if (gequal0(gel(x,1))) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2179 7 : v = gvar(x);
2180 7 : m = n * gval(x,v); break;
2181 0 : default: pari_err_TYPE("O", x);
2182 : v = m = 0; /* LCOV_EXCL_LINE */
2183 : }
2184 3038 : return zeroser(v,m);
2185 : }
2186 :
2187 : /*******************************************************************/
2188 : /* */
2189 : /* FORMAL INTEGRATION */
2190 : /* */
2191 : /*******************************************************************/
2192 : GEN
2193 105 : RgX_integ(GEN x)
2194 : {
2195 105 : long i, lx = lg(x);
2196 : GEN y;
2197 105 : if (lx == 2) return RgX_copy(x);
2198 91 : y = cgetg(lx+1, t_POL); y[1] = x[1]; gel(y,2) = gen_0;
2199 273 : for (i=3; i<=lx; i++) gel(y,i) = gdivgu(gel(x,i-1),i-2);
2200 91 : return y;
2201 : }
2202 :
2203 : static void
2204 35 : err_intformal(GEN x)
2205 35 : { pari_err_DOMAIN("intformal", "residue(series, pole)", "!=", gen_0, x); }
2206 :
2207 : GEN
2208 26040 : integser(GEN x)
2209 : {
2210 26040 : long i, lx = lg(x), vx = varn(x), e = valser(x);
2211 : GEN y;
2212 26040 : if (lx == 2) return zeroser(vx, e+1);
2213 22169 : y = cgetg(lx, t_SER);
2214 96432 : for (i=2; i<lx; i++)
2215 : {
2216 74270 : long j = i+e-1;
2217 74270 : GEN c = gel(x,i);
2218 74270 : if (j)
2219 73955 : c = gdivgs(c, j);
2220 : else
2221 : { /* should be isexactzero, but try to avoid error */
2222 315 : if (!gequal0(c)) err_intformal(x);
2223 308 : c = gen_0;
2224 : }
2225 74263 : gel(y,i) = c;
2226 : }
2227 22162 : y[1] = evalsigne(1) | evalvarn(vx) | evalvalser(e+1); return y;
2228 : }
2229 :
2230 : GEN
2231 350 : integ(GEN x, long v)
2232 : {
2233 350 : long tx = typ(x), vx, n;
2234 350 : pari_sp av = avma;
2235 : GEN y, p1;
2236 :
2237 350 : if (v < 0) { v = gvar9(x); if (v == NO_VARIABLE) v = 0; }
2238 350 : if (is_scalar_t(tx))
2239 : {
2240 63 : if (tx == t_POLMOD)
2241 : {
2242 14 : GEN a = gel(x,2), b = gel(x,1);
2243 14 : vx = varn(b);
2244 14 : if (varncmp(v, vx) > 0) retmkpolmod(integ(a,v), RgX_copy(b));
2245 7 : if (v == vx) pari_err_PRIORITY("intformal",x,"=",v);
2246 : }
2247 49 : return deg1pol(x, gen_0, v);
2248 : }
2249 :
2250 287 : switch(tx)
2251 : {
2252 112 : case t_POL:
2253 112 : vx = varn(x);
2254 112 : if (v == vx) return RgX_integ(x);
2255 42 : if (lg(x) == 2) {
2256 14 : if (varncmp(vx, v) < 0) v = vx;
2257 14 : return zeropol(v);
2258 : }
2259 28 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2260 84 : pari_APPLY_pol(integ(gel(x,i),v));
2261 :
2262 77 : case t_SER:
2263 77 : vx = varn(x);
2264 77 : if (v == vx) return integser(x);
2265 21 : if (lg(x) == 2) {
2266 14 : if (varncmp(vx, v) < 0) v = vx;
2267 14 : return zeroser(v, valser(x));
2268 : }
2269 7 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2270 28 : pari_APPLY_ser(integ(gel(x,i),v));
2271 :
2272 56 : case t_RFRAC:
2273 : {
2274 56 : GEN a = gel(x,1), b = gel(x,2), c, d, s;
2275 56 : vx = varn(b);
2276 56 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2277 49 : if (varncmp(vx, v) < 0)
2278 14 : return gerepileupto(av, swapvar_act(x, vx, v, integ_act, NULL));
2279 :
2280 35 : n = degpol(b);
2281 35 : if (typ(a) == t_POL && varn(a) == vx) n += degpol(a);
2282 35 : y = integ(gadd(x, zeroser(v,n + 2)), v);
2283 35 : y = gdiv(gtrunc(gmul(b, y)), b);
2284 35 : if (typ(y) != t_RFRAC) pari_err_BUG("intformal(t_RFRAC)");
2285 35 : c = gel(y,1); d = gel(y,2);
2286 35 : s = gsub(gmul(deriv(c,v),d), gmul(c,deriv(d,v)));
2287 : /* (c'd-cd')/d^2 = y' = x = a/b ? */
2288 35 : if (!gequal(gmul(s,b), gmul(a,gsqr(d)))) err_intformal(x);
2289 7 : if (typ(y)==t_RFRAC && lg(gel(y,1)) == lg(gel(y,2)))
2290 : {
2291 7 : GEN p2 = leading_coeff(gel(y,2));
2292 7 : p1 = gel(y,1);
2293 7 : if (typ(p1) == t_POL && varn(p1) == vx) p1 = leading_coeff(p1);
2294 7 : y = gsub(y, gdiv(p1,p2));
2295 : }
2296 7 : return gerepileupto(av,y);
2297 : }
2298 :
2299 42 : case t_VEC: case t_COL: case t_MAT:
2300 84 : pari_APPLY_same(integ(gel(x,i),v));
2301 : }
2302 0 : pari_err_TYPE("integ",x);
2303 : return NULL; /* LCOV_EXCL_LINE */
2304 : }
2305 :
2306 : /*******************************************************************/
2307 : /* */
2308 : /* FLOOR */
2309 : /* */
2310 : /*******************************************************************/
2311 : static GEN
2312 518 : quad_floor(GEN x)
2313 : {
2314 518 : GEN Q = gel(x,1), D = quad_disc(x), u, v, b, d, z;
2315 518 : if (signe(D) < 0) return NULL;
2316 490 : x = Q_remove_denom(x, &d);
2317 490 : u = gel(x,2);
2318 490 : v = gel(x,3); b = gel(Q,3);
2319 490 : if (typ(u) != t_INT || typ(v) != t_INT) return NULL;
2320 : /* x0 = (2u + v*(-b + sqrt(D))) / (2d) */
2321 483 : z = sqrti(mulii(D, sqri(v)));
2322 483 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2323 : /* z = floor(v * sqrt(D)) */
2324 483 : z = addii(subii(shifti(u,1), mulii(v,b)), z);
2325 483 : return truedivii(z, d? shifti(d,1): gen_2);
2326 : }
2327 : GEN
2328 5353733 : gfloor(GEN x)
2329 : {
2330 5353733 : switch(typ(x))
2331 : {
2332 5298378 : case t_INT: return icopy(x);
2333 35 : case t_POL: return RgX_copy(x);
2334 37841 : case t_REAL: return floorr(x);
2335 16653 : case t_FRAC: return truedivii(gel(x,1),gel(x,2));
2336 511 : case t_QUAD:
2337 : {
2338 511 : pari_sp av = avma;
2339 : GEN y;
2340 511 : if (!(y = quad_floor(x))) break;
2341 476 : return gerepileuptoint(av, y);
2342 : }
2343 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2344 98 : case t_VEC: case t_COL: case t_MAT:
2345 1533 : pari_APPLY_same(gfloor(gel(x,i)));
2346 : }
2347 238 : pari_err_TYPE("gfloor",x);
2348 : return NULL; /* LCOV_EXCL_LINE */
2349 : }
2350 :
2351 : GEN
2352 426311 : gfrac(GEN x)
2353 : {
2354 : pari_sp av;
2355 : GEN y;
2356 426311 : switch(typ(x))
2357 : {
2358 23529 : case t_INT: return gen_0;
2359 7 : case t_POL: return pol_0(varn(x));
2360 164416 : case t_REAL: av = avma; return gerepileuptoleaf(av, subri(x, floorr(x)));
2361 234629 : case t_FRAC: retmkfrac(modii(gel(x,1),gel(x,2)), icopy(gel(x,2)));
2362 7 : case t_QUAD:
2363 7 : av = avma; if (!(y = quad_floor(x))) break;
2364 7 : return gerepileupto(av, gsub(x, y));
2365 7 : case t_RFRAC: retmkrfrac(grem(gel(x,1),gel(x,2)), gcopy(gel(x,2)));
2366 3688 : case t_VEC: case t_COL: case t_MAT:
2367 34423 : pari_APPLY_same(gfrac(gel(x,i)));
2368 : }
2369 28 : pari_err_TYPE("gfrac",x);
2370 : return NULL; /* LCOV_EXCL_LINE */
2371 : }
2372 :
2373 : /* assume x t_REAL */
2374 : GEN
2375 2454 : ceilr(GEN x)
2376 : {
2377 2454 : pari_sp av = avma;
2378 2454 : GEN y = floorr(x);
2379 2454 : if (cmpri(x, y)) return gerepileuptoint(av, addui(1,y));
2380 29 : return y;
2381 : }
2382 :
2383 : GEN
2384 235223 : gceil(GEN x)
2385 : {
2386 : pari_sp av;
2387 : GEN y;
2388 235223 : switch(typ(x))
2389 : {
2390 17724 : case t_INT: return icopy(x);
2391 21 : case t_POL: return RgX_copy(x);
2392 2376 : case t_REAL: return ceilr(x);
2393 214976 : case t_FRAC:
2394 214976 : av = avma; y = divii(gel(x,1),gel(x,2));
2395 214976 : if (signe(gel(x,1)) > 0) y = gerepileuptoint(av, addui(1,y));
2396 214976 : return y;
2397 49 : case t_QUAD:
2398 49 : if (!is_realquad(x)) break;
2399 42 : if (gequal0(gel(x,3))) return gceil(gel(x,2));
2400 35 : av = avma; return gerepileupto(av, addiu(gfloor(x), 1));
2401 7 : case t_RFRAC:
2402 7 : return gdeuc(gel(x,1),gel(x,2));
2403 35 : case t_VEC: case t_COL: case t_MAT:
2404 105 : pari_APPLY_same(gceil(gel(x,i)));
2405 : }
2406 42 : pari_err_TYPE("gceil",x);
2407 : return NULL; /* LCOV_EXCL_LINE */
2408 : }
2409 :
2410 : GEN
2411 6027 : round0(GEN x, GEN *pte)
2412 : {
2413 6027 : if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
2414 6020 : return ground(x);
2415 : }
2416 :
2417 : /* x t_REAL, return q=floor(x+1/2), set e = expo(x-q) */
2418 : static GEN
2419 53166652 : round_i(GEN x, long *pe)
2420 : {
2421 : long e;
2422 53166652 : GEN B, q,r, m = mantissa_real(x, &e); /* x = m/2^e */
2423 53161097 : if (e <= 0)
2424 : {
2425 13345665 : if (e) m = shifti(m,-e);
2426 13345583 : if (pe) *pe = -e;
2427 13345583 : return m;
2428 : }
2429 39815432 : B = int2n(e-1);
2430 39814612 : m = addii(m, B);
2431 39814691 : q = shifti(m, -e);
2432 39814020 : r = remi2n(m, e);
2433 : /* 2^e (x+1/2) = m = 2^e q + r, sgn(r)=sgn(m), |r|<2^e */
2434 39816753 : if (!signe(r))
2435 24510 : { if (pe) *pe = -1; }
2436 : else
2437 : {
2438 39792243 : if (signe(m) < 0)
2439 : {
2440 15762137 : q = subiu(q,1);
2441 15762125 : r = addii(r, B);
2442 : }
2443 : else
2444 24030106 : r = subii(r, B);
2445 : /* |x - q| = |r| / 2^e */
2446 39792563 : if (pe) *pe = signe(r)? expi(r) - e: -e;
2447 39792412 : cgiv(r);
2448 : }
2449 39819709 : return q;
2450 : }
2451 : /* assume x a t_REAL */
2452 : GEN
2453 3029524 : roundr(GEN x)
2454 : {
2455 3029524 : long ex, s = signe(x);
2456 : pari_sp av;
2457 3029524 : if (!s || (ex=expo(x)) < -1) return gen_0;
2458 2350395 : if (ex == -1) return s>0? gen_1:
2459 182951 : absrnz_equal2n(x)? gen_0: gen_m1;
2460 1758857 : av = avma; x = round_i(x, &ex);
2461 1758829 : if (ex >= 0) pari_err_PREC( "roundr (precision loss in truncation)");
2462 1758829 : return gerepileuptoint(av, x);
2463 : }
2464 : GEN
2465 272620 : roundr_safe(GEN x)
2466 : {
2467 272620 : long ex, s = signe(x);
2468 : pari_sp av;
2469 :
2470 272620 : if (!s || (ex = expo(x)) < -1) return gen_0;
2471 272576 : if (ex == -1) return s>0? gen_1:
2472 0 : absrnz_equal2n(x)? gen_0: gen_m1;
2473 272549 : av = avma; x = round_i(x, NULL);
2474 272549 : return gerepileuptoint(av, x);
2475 : }
2476 :
2477 : GEN
2478 2797646 : ground(GEN x)
2479 : {
2480 : pari_sp av;
2481 : GEN y;
2482 :
2483 2797646 : switch(typ(x))
2484 : {
2485 593942 : case t_INT: return icopy(x);
2486 14 : case t_INTMOD: return gcopy(x);
2487 1642742 : case t_REAL: return roundr(x);
2488 60625 : case t_FRAC: return diviiround(gel(x,1), gel(x,2));
2489 49 : case t_QUAD:
2490 : {
2491 49 : GEN Q = gel(x,1), u, v, b, d, z;
2492 49 : av = avma;
2493 49 : if (is_realquad(x)) return gerepileupto(av, gfloor(gadd(x, ghalf)));
2494 7 : u = gel(x,2);
2495 7 : v = gel(x,3); b = gel(Q,3);
2496 7 : u = ground(gsub(u, gmul2n(gmul(v,b),-1)));
2497 7 : v = Q_remove_denom(v, &d);
2498 7 : if (!d) d = gen_1;
2499 : /* Im x = v sqrt(|D|) / (2d),
2500 : * Im(round(x)) = floor((d + v sqrt(|D|)) / (2d))
2501 : * = floor(floor(d + v sqrt(|D|)) / (2d)) */
2502 7 : z = sqrti(mulii(sqri(v), quad_disc(x)));
2503 7 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2504 : /* z = floor(v * sqrt(|D|)) */
2505 7 : v = truedivii(addii(z, d), shifti(d,1));
2506 7 : return gerepilecopy(av, signe(v)? mkcomplex(u,v): u);
2507 : }
2508 14 : case t_POLMOD:
2509 14 : retmkpolmod(ground(gel(x,2)), RgX_copy(gel(x,1)));
2510 :
2511 4257 : case t_COMPLEX:
2512 4257 : av = avma; y = cgetg(3, t_COMPLEX);
2513 4257 : gel(y,2) = ground(gel(x,2));
2514 4257 : if (!signe(gel(y,2))) { set_avma(av); return ground(gel(x,1)); }
2515 217 : gel(y,1) = ground(gel(x,1)); return y;
2516 :
2517 91 : case t_POL:
2518 4011 : pari_APPLY_pol(ground(gel(x,i)));
2519 182 : case t_SER:
2520 182 : if (ser_isexactzero(x)) return gcopy(x);
2521 42 : pari_APPLY_ser(ground(gel(x,i)));
2522 21 : case t_RFRAC:
2523 21 : av = avma;
2524 21 : return gerepileupto(av, gdiv(ground(gel(x,1)), ground(gel(x,2))));
2525 495706 : case t_VEC: case t_COL: case t_MAT:
2526 2410782 : pari_APPLY_same(ground(gel(x,i)));
2527 : }
2528 6 : pari_err_TYPE("ground",x);
2529 : return NULL; /* LCOV_EXCL_LINE */
2530 : }
2531 :
2532 : /* e = number of error bits on integral part */
2533 : GEN
2534 89041671 : grndtoi(GEN x, long *e)
2535 : {
2536 : GEN y;
2537 : long i, lx, e1;
2538 : pari_sp av;
2539 :
2540 89041671 : if (e) *e = -(long)HIGHEXPOBIT;
2541 89041671 : switch(typ(x))
2542 : {
2543 10024365 : case t_INT: return icopy(x);
2544 55811251 : case t_REAL: {
2545 55811251 : long ex = expo(x);
2546 55811251 : if (!signe(x) || ex < -1)
2547 : {
2548 4675717 : if (e) *e = ex;
2549 4675717 : return gen_0;
2550 : }
2551 51135534 : av = avma; x = round_i(x, e);
2552 51133272 : return gerepileuptoint(av, x);
2553 : }
2554 27653 : case t_FRAC:
2555 27653 : y = diviiround(gel(x,1), gel(x,2)); if (!e) return y;
2556 26498 : av = avma; *e = gexpo(gsub(x, y)); set_avma(av);
2557 26498 : return y;
2558 7 : case t_INTMOD: return gcopy(x);
2559 7 : case t_QUAD:
2560 7 : y = ground(x); av = avma; if (!e) return y;
2561 7 : *e = gexpo(gsub(x,y)); set_avma(avma); return y;
2562 1634127 : case t_COMPLEX:
2563 1634127 : av = avma; y = cgetg(3, t_COMPLEX);
2564 1634127 : gel(y,2) = grndtoi(gel(x,2), e);
2565 1634127 : if (!signe(gel(y,2))) {
2566 244554 : set_avma(av);
2567 244554 : y = grndtoi(gel(x,1), e? &e1: NULL);
2568 : }
2569 : else
2570 1389573 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL);
2571 1634127 : if (e && e1 > *e) *e = e1;
2572 1634127 : return y;
2573 :
2574 7 : case t_POLMOD:
2575 7 : retmkpolmod(grndtoi(gel(x,2), e), RgX_copy(gel(x,1)));
2576 :
2577 318837 : case t_POL:
2578 318837 : y = cgetg_copy(x, &lx); y[1] = x[1];
2579 3246854 : for (i=2; i<lx; i++)
2580 : {
2581 2928018 : gel(y,i) = grndtoi(gel(x,i), &e1);
2582 2928017 : if (e && e1 > *e) *e = e1;
2583 : }
2584 318836 : return normalizepol_lg(y, lx);
2585 168 : case t_SER:
2586 168 : if (ser_isexactzero(x)) return gcopy(x);
2587 154 : y = cgetg_copy(x, &lx); y[1] = x[1];
2588 462 : for (i=2; i<lx; i++)
2589 : {
2590 308 : gel(y,i) = grndtoi(gel(x,i), &e1);
2591 308 : if (e && e1 > *e) *e = e1;
2592 : }
2593 154 : return normalizeser(y);
2594 7 : case t_RFRAC:
2595 7 : y = cgetg(3,t_RFRAC);
2596 7 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2597 7 : gel(y,2) = grndtoi(gel(x,2), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2598 7 : return y;
2599 21226091 : case t_VEC: case t_COL: case t_MAT:
2600 21226091 : y = cgetg_copy(x, &lx);
2601 75134395 : for (i=1; i<lx; i++)
2602 : {
2603 53908722 : gel(y,i) = grndtoi(gel(x,i), e? &e1: NULL);
2604 53908305 : if (e && e1 > *e) *e = e1;
2605 : }
2606 21225673 : return y;
2607 : }
2608 6 : pari_err_TYPE("grndtoi",x);
2609 : return NULL; /* LCOV_EXCL_LINE */
2610 : }
2611 :
2612 : /* trunc(x * 2^s). lindep() sanity checks rely on this function to return a
2613 : * t_INT or fail when fed a non-t_COMPLEX input; so do not make this one
2614 : * recursive [ or change the lindep call ] */
2615 : GEN
2616 116700868 : gtrunc2n(GEN x, long s)
2617 : {
2618 : GEN z;
2619 116700868 : switch(typ(x))
2620 : {
2621 37295431 : case t_INT: return shifti(x, s);
2622 62244876 : case t_REAL: return trunc2nr(x, s);
2623 483 : case t_FRAC: {
2624 : pari_sp av;
2625 483 : GEN a = gel(x,1), b = gel(x,2), q;
2626 483 : if (s == 0) return divii(a, b);
2627 483 : av = avma;
2628 483 : if (s < 0) q = divii(shifti(a, s), b);
2629 : else {
2630 : GEN r;
2631 483 : q = dvmdii(a, b, &r);
2632 483 : q = addii(shifti(q,s), divii(shifti(r,s), b));
2633 : }
2634 483 : return gerepileuptoint(av, q);
2635 : }
2636 17231778 : case t_COMPLEX:
2637 17231778 : z = cgetg(3, t_COMPLEX);
2638 17234848 : gel(z,2) = gtrunc2n(gel(x,2), s);
2639 17232384 : if (!signe(gel(z,2))) {
2640 1595067 : set_avma((pari_sp)(z + 3));
2641 1595061 : return gtrunc2n(gel(x,1), s);
2642 : }
2643 15637317 : gel(z,1) = gtrunc2n(gel(x,1), s);
2644 15635737 : return z;
2645 0 : default: pari_err_TYPE("gtrunc2n",x);
2646 : return NULL; /* LCOV_EXCL_LINE */
2647 : }
2648 : }
2649 :
2650 : /* e = number of error bits on integral part */
2651 : GEN
2652 1115174 : gcvtoi(GEN x, long *e)
2653 : {
2654 1115174 : long tx = typ(x), lx, e1;
2655 : GEN y;
2656 :
2657 1115174 : if (tx == t_REAL)
2658 : {
2659 1115041 : long ex = expo(x); if (ex < 0) { *e = ex; return gen_0; }
2660 1114936 : e1 = ex - bit_prec(x) + 1;
2661 1114937 : y = mantissa2nr(x, e1);
2662 1114933 : if (e1 <= 0) { pari_sp av = avma; e1 = expo(subri(x,y)); set_avma(av); }
2663 1114928 : *e = e1; return y;
2664 : }
2665 133 : *e = -(long)HIGHEXPOBIT;
2666 133 : if (is_matvec_t(tx))
2667 : {
2668 : long i;
2669 28 : y = cgetg_copy(x, &lx);
2670 84 : for (i=1; i<lx; i++)
2671 : {
2672 56 : gel(y,i) = gcvtoi(gel(x,i),&e1);
2673 56 : if (e1 > *e) *e = e1;
2674 : }
2675 28 : return y;
2676 : }
2677 105 : return gtrunc(x);
2678 : }
2679 :
2680 : int
2681 445876 : isint(GEN n, GEN *ptk)
2682 : {
2683 445876 : switch(typ(n))
2684 : {
2685 396792 : case t_INT: *ptk = n; return 1;
2686 1330 : case t_REAL: {
2687 1330 : pari_sp av0 = avma;
2688 1330 : GEN z = floorr(n);
2689 1330 : pari_sp av = avma;
2690 1330 : long s = signe(subri(n, z));
2691 1330 : if (s) return gc_bool(av0,0);
2692 21 : *ptk = z; return gc_bool(av,1);
2693 : }
2694 30604 : case t_FRAC: return 0;
2695 16961 : case t_COMPLEX: return gequal0(gel(n,2)) && isint(gel(n,1),ptk);
2696 0 : case t_QUAD: return gequal0(gel(n,3)) && isint(gel(n,2),ptk);
2697 189 : default: pari_err_TYPE("isint",n);
2698 : return 0; /* LCOV_EXCL_LINE */
2699 : }
2700 : }
2701 :
2702 : int
2703 326463 : issmall(GEN n, long *ptk)
2704 : {
2705 326463 : pari_sp av = avma;
2706 : GEN z;
2707 : long k;
2708 326463 : if (!isint(n, &z)) return 0;
2709 324881 : k = itos_or_0(z); set_avma(av);
2710 324881 : if (k || lgefint(z) == 2) { *ptk = k; return 1; }
2711 0 : return 0;
2712 : }
2713 :
2714 : /* smallest integer greater than any incarnations of the real x
2715 : * Avoid "precision loss in truncation" */
2716 : GEN
2717 1041765 : ceil_safe(GEN x)
2718 : {
2719 1041765 : pari_sp av = avma;
2720 1041765 : long e, tx = typ(x);
2721 : GEN y;
2722 :
2723 1041765 : if (is_rational_t(tx)) return gceil(x);
2724 1041470 : y = gcvtoi(x,&e);
2725 1041465 : if (gsigne(x) >= 0)
2726 : {
2727 1040967 : if (e < 0) e = 0;
2728 1040967 : y = addii(y, int2n(e));
2729 : }
2730 1041467 : return gerepileuptoint(av, y);
2731 : }
2732 : /* largest integer smaller than any incarnations of the real x
2733 : * Avoid "precision loss in truncation" */
2734 : GEN
2735 61927 : floor_safe(GEN x)
2736 : {
2737 61927 : pari_sp av = avma;
2738 61927 : long e, tx = typ(x);
2739 : GEN y;
2740 :
2741 61927 : if (is_rational_t(tx)) return gfloor(x);
2742 61927 : y = gcvtoi(x,&e);
2743 61926 : if (gsigne(x) <= 0)
2744 : {
2745 21 : if (e < 0) e = 0;
2746 21 : y = subii(y, int2n(e));
2747 : }
2748 61926 : return gerepileuptoint(av, y);
2749 : }
2750 :
2751 : GEN
2752 10801 : ser2rfrac_i(GEN x)
2753 : {
2754 10801 : long e = valser(x);
2755 10801 : GEN a = ser2pol_i(x, lg(x));
2756 10801 : if (e) {
2757 6370 : if (e > 0) a = RgX_shift_shallow(a, e);
2758 0 : else a = gred_rfrac_simple(a, pol_xn(-e, varn(a)));
2759 : }
2760 10801 : return a;
2761 : }
2762 :
2763 : static GEN
2764 539 : ser2rfrac(GEN x)
2765 : {
2766 539 : pari_sp av = avma;
2767 539 : return gerepilecopy(av, ser2rfrac_i(x));
2768 : }
2769 :
2770 : /* x t_PADIC, truncate to rational (t_INT/t_FRAC) */
2771 : GEN
2772 689401 : padic_to_Q(GEN x)
2773 : {
2774 689401 : GEN u = gel(x,4), p;
2775 : long v;
2776 689401 : if (!signe(u)) return gen_0;
2777 683339 : v = valp(x);
2778 683339 : if (!v) return icopy(u);
2779 258499 : p = gel(x,2);
2780 258499 : if (v>0)
2781 : {
2782 258383 : pari_sp av = avma;
2783 258383 : return gerepileuptoint(av, mulii(u, powiu(p,v)));
2784 : }
2785 116 : retmkfrac(icopy(u), powiu(p,-v));
2786 : }
2787 : GEN
2788 14 : padic_to_Q_shallow(GEN x)
2789 : {
2790 14 : GEN u = gel(x,4), p, q, q2;
2791 : long v;
2792 14 : if (!signe(u)) return gen_0;
2793 14 : q = gel(x,3); q2 = shifti(q,-1);
2794 14 : v = valp(x);
2795 14 : u = Fp_center_i(u, q, q2);
2796 14 : if (!v) return u;
2797 7 : p = gel(x,2);
2798 7 : if (v>0) return mulii(powiu(p,v), u);
2799 7 : return mkfrac(u, powiu(p,-v));
2800 : }
2801 : GEN
2802 196 : QpV_to_QV(GEN v)
2803 : {
2804 : long i, l;
2805 196 : GEN w = cgetg_copy(v, &l);
2806 931 : for (i = 1; i < l; i++)
2807 : {
2808 735 : GEN c = gel(v,i);
2809 735 : switch(typ(c))
2810 : {
2811 721 : case t_INT:
2812 721 : case t_FRAC: break;
2813 14 : case t_PADIC: c = padic_to_Q_shallow(c); break;
2814 0 : default: pari_err_TYPE("padic_to_Q", v);
2815 : }
2816 735 : gel(w,i) = c;
2817 : }
2818 196 : return w;
2819 : }
2820 :
2821 : GEN
2822 597616 : gtrunc(GEN x)
2823 : {
2824 597616 : switch(typ(x))
2825 : {
2826 84175 : case t_INT: return icopy(x);
2827 49098 : case t_REAL: return truncr(x);
2828 56 : case t_FRAC: return divii(gel(x,1),gel(x,2));
2829 432731 : case t_PADIC: return padic_to_Q(x);
2830 42 : case t_POL: return RgX_copy(x);
2831 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2832 511 : case t_SER: return ser2rfrac(x);
2833 187215 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gtrunc(gel(x,i)));
2834 : }
2835 56 : pari_err_TYPE("gtrunc",x);
2836 : return NULL; /* LCOV_EXCL_LINE */
2837 : }
2838 :
2839 : GEN
2840 217 : trunc0(GEN x, GEN *pte)
2841 : {
2842 217 : if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
2843 189 : return gtrunc(x);
2844 : }
2845 : /*******************************************************************/
2846 : /* */
2847 : /* CONVERSIONS --> INT, POL & SER */
2848 : /* */
2849 : /*******************************************************************/
2850 :
2851 : /* return a_(n-1) B^(n-1) + ... + a_0, where B = 2^32.
2852 : * The a_i are 32bits integers */
2853 : GEN
2854 24076 : mkintn(long n, ...)
2855 : {
2856 : va_list ap;
2857 : GEN x, y;
2858 : long i;
2859 : #ifdef LONG_IS_64BIT
2860 20670 : long e = (n&1);
2861 20670 : n = (n+1) >> 1;
2862 : #endif
2863 24076 : va_start(ap,n);
2864 24076 : x = cgetipos(n+2);
2865 24076 : y = int_MSW(x);
2866 85138 : for (i=0; i <n; i++)
2867 : {
2868 : #ifdef LONG_IS_64BIT
2869 47700 : ulong a = (e && !i)? 0: (ulong) va_arg(ap, unsigned int);
2870 47700 : ulong b = (ulong) va_arg(ap, unsigned int);
2871 47700 : *y = (a << 32) | b;
2872 : #else
2873 13362 : *y = (ulong) va_arg(ap, unsigned int);
2874 : #endif
2875 61062 : y = int_precW(y);
2876 : }
2877 24076 : va_end(ap);
2878 24076 : return int_normalize(x, 0);
2879 : }
2880 :
2881 : /* 2^32 a + b */
2882 : GEN
2883 243465 : uu32toi(ulong a, ulong b)
2884 : {
2885 : #ifdef LONG_IS_64BIT
2886 196156 : return utoi((a<<32) | b);
2887 : #else
2888 47309 : return uutoi(a, b);
2889 : #endif
2890 : }
2891 : /* - (2^32 a + b), assume a or b != 0 */
2892 : GEN
2893 72146 : uu32toineg(ulong a, ulong b)
2894 : {
2895 : #ifdef LONG_IS_64BIT
2896 61721 : return utoineg((a<<32) | b);
2897 : #else
2898 10425 : return uutoineg(a, b);
2899 : #endif
2900 : }
2901 :
2902 : /* return a_(n-1) x^(n-1) + ... + a_0 */
2903 : GEN
2904 5553047 : mkpoln(long n, ...)
2905 : {
2906 : va_list ap;
2907 : GEN x, y;
2908 5553047 : long i, l = n+2;
2909 5553047 : va_start(ap,n);
2910 5553047 : x = cgetg(l, t_POL); y = x + 2;
2911 5554408 : x[1] = evalvarn(0);
2912 24437469 : for (i=n-1; i >= 0; i--) gel(y,i) = va_arg(ap, GEN);
2913 5558428 : va_end(ap); return normalizepol_lg(x, l);
2914 : }
2915 :
2916 : /* return [a_1, ..., a_n] */
2917 : GEN
2918 1138120 : mkvecn(long n, ...)
2919 : {
2920 : va_list ap;
2921 : GEN x;
2922 : long i;
2923 1138120 : va_start(ap,n);
2924 1138120 : x = cgetg(n+1, t_VEC);
2925 7295252 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2926 1138120 : va_end(ap); return x;
2927 : }
2928 :
2929 : GEN
2930 1379 : mkcoln(long n, ...)
2931 : {
2932 : va_list ap;
2933 : GEN x;
2934 : long i;
2935 1379 : va_start(ap,n);
2936 1379 : x = cgetg(n+1, t_COL);
2937 11032 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2938 1379 : va_end(ap); return x;
2939 : }
2940 :
2941 : GEN
2942 126181 : mkvecsmalln(long n, ...)
2943 : {
2944 : va_list ap;
2945 : GEN x;
2946 : long i;
2947 126181 : va_start(ap,n);
2948 126181 : x = cgetg(n+1, t_VECSMALL);
2949 878540 : for (i=1; i <= n; i++) x[i] = va_arg(ap, long);
2950 126181 : va_end(ap); return x;
2951 : }
2952 :
2953 : GEN
2954 16902078 : scalarpol(GEN x, long v)
2955 : {
2956 : GEN y;
2957 16902078 : if (isrationalzero(x)) return zeropol(v);
2958 5140839 : y = cgetg(3,t_POL);
2959 5140865 : y[1] = gequal0(x)? evalvarn(v)
2960 5140862 : : evalvarn(v) | evalsigne(1);
2961 5140862 : gel(y,2) = gcopy(x); return y;
2962 : }
2963 : GEN
2964 3673350 : scalarpol_shallow(GEN x, long v)
2965 : {
2966 : GEN y;
2967 3673350 : if (isrationalzero(x)) return zeropol(v);
2968 3564210 : y = cgetg(3,t_POL);
2969 3564212 : y[1] = gequal0(x)? evalvarn(v)
2970 3564212 : : evalvarn(v) | evalsigne(1);
2971 3564212 : gel(y,2) = x; return y;
2972 : }
2973 :
2974 : /* x0 + x1*T, do not assume x1 != 0 */
2975 : GEN
2976 559968 : deg1pol(GEN x1, GEN x0,long v)
2977 : {
2978 559968 : GEN x = cgetg(4,t_POL);
2979 559967 : x[1] = evalsigne(1) | evalvarn(v);
2980 559967 : gel(x,2) = x0 == gen_0? x0: gcopy(x0); /* gen_0 frequent */
2981 559967 : gel(x,3) = gcopy(x1); return normalizepol_lg(x,4);
2982 : }
2983 :
2984 : /* same, no copy */
2985 : GEN
2986 8516454 : deg1pol_shallow(GEN x1, GEN x0,long v)
2987 : {
2988 8516454 : GEN x = cgetg(4,t_POL);
2989 8516531 : x[1] = evalsigne(1) | evalvarn(v);
2990 8516531 : gel(x,2) = x0;
2991 8516531 : gel(x,3) = x1; return normalizepol_lg(x,4);
2992 : }
2993 :
2994 : GEN
2995 319062 : deg2pol_shallow(GEN x2, GEN x1, GEN x0, long v)
2996 : {
2997 319062 : GEN x = cgetg(5,t_POL);
2998 319062 : x[1] = evalsigne(1) | evalvarn(v);
2999 319062 : gel(x,2) = x0;
3000 319062 : gel(x,3) = x1;
3001 319062 : gel(x,4) = x2;
3002 319062 : return normalizepol_lg(x,5);
3003 : }
3004 :
3005 : static GEN
3006 3469535 : _gtopoly(GEN x, long v, int reverse)
3007 : {
3008 3469535 : long tx = typ(x);
3009 : GEN y;
3010 :
3011 3469535 : if (v<0) v = 0;
3012 3469535 : switch(tx)
3013 : {
3014 21 : case t_POL:
3015 21 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3016 21 : y = RgX_copy(x); break;
3017 28 : case t_SER:
3018 28 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3019 28 : y = ser2rfrac(x);
3020 28 : if (typ(y) != t_POL)
3021 0 : pari_err_DOMAIN("gtopoly", "valuation", "<", gen_0, x);
3022 28 : break;
3023 42 : case t_RFRAC:
3024 : {
3025 42 : GEN a = gel(x,1), b = gel(x,2);
3026 42 : long vb = varn(b);
3027 42 : if (varncmp(vb, v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3028 42 : if (typ(a) != t_POL || varn(a) != vb) return zeropol(v);
3029 21 : y = RgX_div(a,b); break;
3030 : }
3031 337043 : case t_VECSMALL: x = zv_to_ZV(x); /* fall through */
3032 3468877 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
3033 : {
3034 3468877 : long j, k, lx = lg(x);
3035 : GEN z;
3036 3468877 : if (tx == t_QFB) lx--;
3037 3468877 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("gtopoly", x, "<=", v);
3038 3468883 : y = cgetg(lx+1, t_POL);
3039 3468881 : y[1] = evalvarn(v);
3040 3468881 : if (reverse) {
3041 2434012 : x--;
3042 26224492 : for (j=2; j<=lx; j++) gel(y,j) = gel(x,j);
3043 : } else {
3044 5297813 : for (j=2, k=lx; k>=2; j++) gel(y,j) = gel(x,--k);
3045 : }
3046 3468881 : z = RgX_copy(normalizepol_lg(y,lx+1));
3047 3468884 : settyp(y, t_VECSMALL);/* left on stack */
3048 3468884 : return z;
3049 : }
3050 567 : default:
3051 567 : if (is_scalar_t(tx)) return scalarpol(x,v);
3052 7 : pari_err_TYPE("gtopoly",x);
3053 : return NULL; /* LCOV_EXCL_LINE */
3054 : }
3055 70 : setvarn(y,v); return y;
3056 : }
3057 :
3058 : GEN
3059 2434047 : gtopolyrev(GEN x, long v) { return _gtopoly(x,v,1); }
3060 :
3061 : GEN
3062 1035486 : gtopoly(GEN x, long v) { return _gtopoly(x,v,0); }
3063 :
3064 : static GEN
3065 1092 : gtovecpost(GEN x, long n)
3066 : {
3067 1092 : long i, imax, lx, tx = typ(x);
3068 1092 : GEN y = zerovec(n);
3069 :
3070 1092 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,1) = gcopy(x); return y; }
3071 343 : switch(tx)
3072 : {
3073 56 : case t_POL:
3074 56 : lx = lg(x); imax = minss(lx-2, n);
3075 224 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,lx-i));
3076 56 : return y;
3077 28 : case t_SER:
3078 28 : lx = lg(x); imax = minss(lx-2, n); x++;
3079 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3080 28 : return y;
3081 28 : case t_QFB:
3082 28 : lx = lg(x)-1; /* remove discriminant */
3083 28 : imax = minss(lx-1, n);
3084 112 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3085 28 : return y;
3086 28 : case t_VEC: case t_COL:
3087 28 : lx = lg(x); imax = minss(lx-1, n);
3088 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3089 28 : return y;
3090 63 : case t_LIST:
3091 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3092 56 : x = list_data(x); lx = x? lg(x): 1;
3093 56 : imax = minss(lx-1, n);
3094 252 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3095 56 : return y;
3096 56 : case t_STR:
3097 : {
3098 56 : char *s = GSTR(x);
3099 56 : imax = minss(strlen(s), n); s--;
3100 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3101 56 : return y;
3102 : }
3103 28 : case t_VECSMALL:
3104 28 : lx = lg(x); imax = minss(lx-1, n);
3105 84 : for (i=1; i<=imax; i++) gel(y,i) = stoi(x[i]);
3106 28 : return y;
3107 56 : default: pari_err_TYPE("gtovec",x);
3108 : return NULL; /*LCOV_EXCL_LINE*/
3109 : }
3110 : }
3111 :
3112 : static GEN
3113 3563 : init_vectopre(long a, long n, GEN y, long *imax)
3114 : {
3115 3563 : if (n <= a) { *imax = n; return y; }
3116 2765 : *imax = a; return y + n - a;
3117 : }
3118 : /* assume n > 0 */
3119 : static GEN
3120 3619 : gtovecpre(GEN x, long n)
3121 : {
3122 3619 : long a, i, imax, lx, tx = typ(x);
3123 3619 : GEN y = zerovec(n), y0;
3124 :
3125 3619 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,n) = gcopy(x); return y; }
3126 3563 : switch(tx)
3127 : {
3128 56 : case t_POL:
3129 56 : lx = lg(x); a = lx-2;
3130 56 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x -= a-imax;
3131 224 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,lx-i));
3132 56 : return y;
3133 3248 : case t_SER:
3134 3248 : a = lg(x)-2; x++;
3135 3248 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3136 131425 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3137 3248 : return y;
3138 28 : case t_QFB:
3139 28 : a = lg(x)-2; /* remove discriminant */
3140 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3141 112 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3142 28 : return y;
3143 28 : case t_VEC: case t_COL:
3144 28 : a = lg(x)-1;
3145 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3146 84 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3147 28 : return y;
3148 63 : case t_LIST:
3149 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3150 56 : x = list_data(x); a = x? lg(x)-1: 0;
3151 56 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3152 252 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3153 56 : return y;
3154 56 : case t_STR:
3155 : {
3156 56 : char *s = GSTR(x);
3157 56 : a = strlen(s);
3158 56 : y0 = init_vectopre(a, n, y, &imax); s--; if (imax == n) s += a-imax;
3159 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3160 56 : return y;
3161 : }
3162 28 : case t_VECSMALL:
3163 28 : a = lg(x)-1;
3164 28 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3165 84 : for (i=1; i<=imax; i++) gel(y0,i) = stoi(x[i]);
3166 28 : return y;
3167 56 : default: pari_err_TYPE("gtovec",x);
3168 : return NULL; /*LCOV_EXCL_LINE*/
3169 : }
3170 : }
3171 : GEN
3172 106133 : gtovec0(GEN x, long n)
3173 : {
3174 106133 : if (!n) return gtovec(x);
3175 4711 : if (n > 0) return gtovecpost(x, n);
3176 3619 : return gtovecpre(x, -n);
3177 : }
3178 :
3179 : GEN
3180 101772 : gtovec(GEN x)
3181 : {
3182 101772 : long i, lx, tx = typ(x);
3183 : GEN y;
3184 :
3185 101772 : if (is_scalar_t(tx)) return mkveccopy(x);
3186 101639 : switch(tx)
3187 : {
3188 15232 : case t_POL:
3189 15232 : lx=lg(x); y=cgetg(lx-1,t_VEC);
3190 1498371 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,lx-i));
3191 15232 : return y;
3192 385 : case t_SER:
3193 385 : lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
3194 12264 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
3195 385 : return y;
3196 28 : case t_RFRAC: return mkveccopy(x);
3197 70049 : case t_QFB:
3198 70049 : retmkvec3(icopy(gel(x,1)), icopy(gel(x,2)), icopy(gel(x,3)));
3199 14203 : case t_VEC: case t_COL: case t_MAT:
3200 14203 : lx=lg(x); y=cgetg(lx,t_VEC);
3201 1467354 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3202 14203 : return y;
3203 76 : case t_LIST:
3204 76 : if (list_typ(x) == t_LIST_MAP) return mapdomain(x);
3205 62 : x = list_data(x); lx = x? lg(x): 1;
3206 62 : y = cgetg(lx, t_VEC);
3207 304 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3208 62 : return y;
3209 105 : case t_STR:
3210 : {
3211 105 : char *s = GSTR(x);
3212 105 : lx = strlen(s)+1; y = cgetg(lx, t_VEC);
3213 105 : s--;
3214 340239 : for (i=1; i<lx; i++) gel(y,i) = chartoGENstr(s[i]);
3215 105 : return y;
3216 : }
3217 1498 : case t_VECSMALL:
3218 1498 : return vecsmall_to_vec(x);
3219 63 : case t_ERROR:
3220 63 : lx=lg(x); y = cgetg(lx,t_VEC);
3221 63 : gel(y,1) = errname(x);
3222 168 : for (i=2; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3223 63 : return y;
3224 0 : default: pari_err_TYPE("gtovec",x);
3225 : return NULL; /*LCOV_EXCL_LINE*/
3226 : }
3227 : }
3228 :
3229 : GEN
3230 308 : gtovecrev0(GEN x, long n)
3231 : {
3232 308 : GEN y = gtovec0(x, -n);
3233 280 : vecreverse_inplace(y);
3234 280 : return y;
3235 : }
3236 : GEN
3237 0 : gtovecrev(GEN x) { return gtovecrev0(x, 0); }
3238 :
3239 : GEN
3240 3668 : gtocol0(GEN x, long n)
3241 : {
3242 : GEN y;
3243 3668 : if (!n) return gtocol(x);
3244 3458 : y = gtovec0(x, n);
3245 3402 : settyp(y, t_COL); return y;
3246 : }
3247 : GEN
3248 210 : gtocol(GEN x)
3249 : {
3250 : long lx, tx, i, j, h;
3251 : GEN y;
3252 210 : tx = typ(x);
3253 210 : if (tx != t_MAT) { y = gtovec(x); settyp(y, t_COL); return y; }
3254 14 : lx = lg(x); if (lx == 1) return cgetg(1, t_COL);
3255 14 : h = lgcols(x); y = cgetg(h, t_COL);
3256 42 : for (i = 1 ; i < h; i++) {
3257 28 : gel(y,i) = cgetg(lx, t_VEC);
3258 112 : for (j = 1; j < lx; j++) gmael(y,i,j) = gcopy(gcoeff(x,i,j));
3259 : }
3260 14 : return y;
3261 : }
3262 :
3263 : GEN
3264 294 : gtocolrev0(GEN x, long n)
3265 : {
3266 294 : GEN y = gtocol0(x, -n);
3267 266 : long ly = lg(y), lim = ly>>1, i;
3268 763 : for (i = 1; i <= lim; i++) swap(gel(y,i), gel(y,ly-i));
3269 266 : return y;
3270 : }
3271 : GEN
3272 0 : gtocolrev(GEN x) { return gtocolrev0(x, 0); }
3273 :
3274 : static long
3275 19005 : Itos(GEN x)
3276 : {
3277 19005 : if (typ(x) != t_INT) pari_err_TYPE("vectosmall",x);
3278 19005 : return itos(x);
3279 : }
3280 :
3281 : static GEN
3282 98 : gtovecsmallpost(GEN x, long n)
3283 : {
3284 : long i, imax, lx;
3285 98 : GEN y = zero_Flv(n);
3286 :
3287 98 : switch(typ(x))
3288 : {
3289 7 : case t_INT:
3290 7 : y[1] = itos(x); return y;
3291 14 : case t_POL:
3292 14 : lx=lg(x); imax = minss(lx-2, n);
3293 56 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,lx-i));
3294 14 : return y;
3295 7 : case t_SER:
3296 7 : lx=lg(x); imax = minss(lx-2, n); x++;
3297 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3298 7 : return y;
3299 7 : case t_VEC: case t_COL:
3300 7 : lx=lg(x); imax = minss(lx-1, n);
3301 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3302 7 : return y;
3303 14 : case t_LIST:
3304 14 : x = list_data(x); lx = x? lg(x): 1;
3305 14 : imax = minss(lx-1, n);
3306 63 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3307 14 : return y;
3308 7 : case t_VECSMALL:
3309 7 : lx=lg(x);
3310 7 : imax = minss(lx-1, n);
3311 21 : for (i=1; i<=imax; i++) y[i] = x[i];
3312 7 : return y;
3313 14 : case t_STR:
3314 : {
3315 14 : unsigned char *s = (unsigned char*)GSTR(x);
3316 14 : imax = minss(strlen((const char *)s), n); s--;
3317 56 : for (i=1; i<=imax; i++) y[i] = (long)s[i];
3318 14 : return y;
3319 : }
3320 28 : default: pari_err_TYPE("gtovecsmall",x);
3321 : return NULL; /*LCOV_EXCL_LINE*/
3322 : }
3323 : }
3324 : static GEN
3325 98 : gtovecsmallpre(GEN x, long n)
3326 : {
3327 98 : GEN y = zero_Flv(n), y0;
3328 : long a, i, imax, lx;
3329 :
3330 98 : switch(typ(x))
3331 : {
3332 7 : case t_INT:
3333 7 : y[n] = itos(x); return y;
3334 14 : case t_POL:
3335 14 : lx = lg(x); a = lx-2;
3336 14 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x -= a-imax;
3337 56 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,lx-i));
3338 14 : return y;
3339 7 : case t_SER:
3340 7 : a = lg(x)-2; x++;
3341 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3342 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3343 7 : return y;
3344 7 : case t_VEC: case t_COL:
3345 7 : a = lg(x)-1;
3346 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3347 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3348 7 : return y;
3349 14 : case t_LIST:
3350 14 : x = list_data(x); a = x? lg(x)-1: 0;
3351 14 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3352 63 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3353 14 : return y;
3354 7 : case t_VECSMALL:
3355 7 : a = lg(x)-1;
3356 7 : y0 = init_vectopre(a, n, y, &imax); if (imax == n) x += a-imax;
3357 21 : for (i=1; i<=imax; i++) y0[i] = x[i];
3358 7 : return y;
3359 14 : case t_STR:
3360 : {
3361 14 : unsigned char *s = (unsigned char*)GSTR(x);
3362 14 : a = strlen((const char *)s);
3363 14 : y0 = init_vectopre(a, n, y, &imax); s--; if (imax == n) s += a-imax;
3364 56 : for (i=1; i<=imax; i++) y0[i] = (long)s[i];
3365 14 : return y;
3366 : }
3367 28 : default: pari_err_TYPE("gtovecsmall",x);
3368 : return NULL; /*LCOV_EXCL_LINE*/
3369 : }
3370 : }
3371 :
3372 : GEN
3373 7728 : gtovecsmall0(GEN x, long n)
3374 : {
3375 7728 : if (!n) return gtovecsmall(x);
3376 196 : if (n > 0) return gtovecsmallpost(x, n);
3377 98 : return gtovecsmallpre(x, -n);
3378 : }
3379 :
3380 : GEN
3381 18970 : gtovecsmall(GEN x)
3382 : {
3383 : GEN V;
3384 : long l, i;
3385 :
3386 18970 : switch(typ(x))
3387 : {
3388 112 : case t_INT: return mkvecsmall(itos(x));
3389 28 : case t_STR:
3390 : {
3391 28 : unsigned char *s = (unsigned char*)GSTR(x);
3392 28 : l = strlen((const char *)s);
3393 28 : V = cgetg(l+1, t_VECSMALL);
3394 28 : s--;
3395 1953 : for (i=1; i<=l; i++) V[i] = (long)s[i];
3396 28 : return V;
3397 : }
3398 13139 : case t_VECSMALL: return leafcopy(x);
3399 14 : case t_LIST:
3400 14 : x = list_data(x);
3401 14 : if (!x) return cgetg(1, t_VECSMALL);
3402 : /* fall through */
3403 : case t_VEC: case t_COL:
3404 5642 : l = lg(x); V = cgetg(l,t_VECSMALL);
3405 24346 : for(i=1; i<l; i++) V[i] = Itos(gel(x,i));
3406 5642 : return V;
3407 14 : case t_POL:
3408 14 : l = lg(x); V = cgetg(l-1,t_VECSMALL);
3409 63 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,l-i));
3410 14 : return V;
3411 7 : case t_SER:
3412 7 : l = lg(x); V = cgetg(l-1,t_VECSMALL); x++;
3413 21 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,i));
3414 7 : return V;
3415 28 : default:
3416 28 : pari_err_TYPE("vectosmall",x);
3417 : return NULL; /* LCOV_EXCL_LINE */
3418 : }
3419 : }
3420 :
3421 : GEN
3422 327 : compo(GEN x, long n)
3423 : {
3424 327 : long tx = typ(x);
3425 327 : ulong l, lx = (ulong)lg(x);
3426 :
3427 327 : if (!is_recursive_t(tx))
3428 : {
3429 63 : if (tx == t_VECSMALL)
3430 : {
3431 21 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3432 21 : if ((ulong)n >= lx) pari_err_COMPONENT("", ">", utoi(lx-1), stoi(n));
3433 7 : return stoi(x[n]);
3434 : }
3435 42 : pari_err_TYPE("component [leaf]", x);
3436 : }
3437 264 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3438 257 : if (tx == t_LIST) {
3439 28 : x = list_data(x); tx = t_VEC;
3440 28 : lx = (ulong)(x? lg(x): 1);
3441 : }
3442 257 : l = (ulong)lontyp[tx] + (ulong)n-1; /* beware overflow */
3443 257 : if (l >= lx) pari_err_COMPONENT("", ">", utoi(lx-lontyp[tx]), stoi(n));
3444 173 : return gcopy(gel(x,l));
3445 : }
3446 :
3447 : /* assume x a t_POL */
3448 : static GEN
3449 2253279 : _polcoef(GEN x, long n, long v)
3450 : {
3451 2253279 : long i, w, lx = lg(x), dx = lx-3;
3452 : GEN z;
3453 2253279 : if (dx < 0) return gen_0;
3454 1612065 : if (v < 0 || v == (w=varn(x)))
3455 1238953 : return (n < 0 || n > dx)? gen_0: gel(x,n+2);
3456 373112 : if (varncmp(w,v) > 0) return n? gen_0: x;
3457 : /* w < v */
3458 372273 : z = cgetg(lx, t_POL); z[1] = x[1];
3459 1795912 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3460 372271 : z = normalizepol_lg(z, lx);
3461 372271 : switch(lg(z))
3462 : {
3463 27952 : case 2: z = gen_0; break;
3464 125063 : case 3: z = gel(z,2); break;
3465 : }
3466 372271 : return z;
3467 : }
3468 :
3469 : /* assume x a t_SER */
3470 : static GEN
3471 111902 : _sercoef(GEN x, long n, long v)
3472 : {
3473 111902 : long i, w = varn(x), lx = lg(x), dx = lx-3, N;
3474 : GEN z;
3475 111902 : if (v < 0) v = w;
3476 111902 : N = v == w? n - valser(x): n;
3477 111902 : if (dx < 0)
3478 : {
3479 21 : if (N >= 0) pari_err_DOMAIN("polcoef", "t_SER", "=", x, x);
3480 14 : return gen_0;
3481 : }
3482 111881 : if (v == w)
3483 : {
3484 111839 : if (!dx && !signe(x) && !isinexact(gel(x,2))) dx = -1;
3485 111839 : if (N > dx)
3486 28 : pari_err_DOMAIN("polcoef", "degree", ">", stoi(dx+valser(x)), stoi(n));
3487 111811 : return (N < 0)? gen_0: gel(x,N+2);
3488 : }
3489 42 : if (varncmp(w,v) > 0) return N? gen_0: x;
3490 : /* w < v */
3491 28 : z = cgetg(lx, t_SER); z[1] = x[1];
3492 91 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3493 28 : return normalizeser(z);
3494 : }
3495 :
3496 : /* assume x a t_RFRAC(n) */
3497 : static GEN
3498 21 : _rfraccoef(GEN x, long n, long v)
3499 : {
3500 21 : GEN P,Q, p = gel(x,1), q = gel(x,2);
3501 21 : long vp = gvar(p), vq = gvar(q);
3502 21 : if (v < 0) v = varncmp(vp, vq) < 0? vp: vq;
3503 21 : P = (vp == v)? p: swap_vars(p, v);
3504 21 : Q = (vq == v)? q: swap_vars(q, v);
3505 21 : if (!RgX_is_monomial(Q)) pari_err_TYPE("polcoef", x);
3506 21 : n += degpol(Q);
3507 21 : return gdiv(_polcoef(P, n, v), leading_coeff(Q));
3508 : }
3509 :
3510 : GEN
3511 2954033 : polcoef_i(GEN x, long n, long v)
3512 : {
3513 2954033 : long tx = typ(x);
3514 2954033 : switch(tx)
3515 : {
3516 2253258 : case t_POL: return _polcoef(x,n,v);
3517 111902 : case t_SER: return _sercoef(x,n,v);
3518 21 : case t_RFRAC: return _rfraccoef(x,n,v);
3519 : }
3520 588852 : if (!is_scalar_t(tx)) pari_err_TYPE("polcoef", x);
3521 588678 : return n? gen_0: x;
3522 : }
3523 :
3524 : /* with respect to the main variable if v<0, with respect to the variable v
3525 : * otherwise. v ignored if x is not a polynomial/series. */
3526 : GEN
3527 727580 : polcoef(GEN x, long n, long v)
3528 : {
3529 727580 : pari_sp av = avma;
3530 727580 : x = polcoef_i(x,n,v);
3531 727349 : if (x == gen_0) return x;
3532 130725 : return (avma == av)? gcopy(x): gerepilecopy(av, x);
3533 : }
3534 :
3535 : static GEN
3536 111541 : vecdenom(GEN v, long imin, long imax)
3537 : {
3538 111541 : long i = imin;
3539 : GEN s;
3540 111541 : if (imin > imax) return gen_1;
3541 111541 : s = denom_i(gel(v,i));
3542 223663 : for (i++; i<=imax; i++)
3543 : {
3544 112122 : GEN t = denom_i(gel(v,i));
3545 112122 : if (t != gen_1) s = glcm(s,t);
3546 : }
3547 111541 : return s;
3548 : }
3549 : static GEN denompol(GEN x, long v);
3550 : static GEN
3551 14 : vecdenompol(GEN v, long imin, long imax, long vx)
3552 : {
3553 14 : long i = imin;
3554 : GEN s;
3555 14 : if (imin > imax) return gen_1;
3556 14 : s = denompol(gel(v,i), vx);
3557 14 : for (i++; i<=imax; i++)
3558 : {
3559 0 : GEN t = denompol(gel(v,i), vx);
3560 0 : if (t != gen_1) s = glcm(s,t);
3561 : }
3562 14 : return s;
3563 : }
3564 : GEN
3565 10349480 : denom_i(GEN x)
3566 : {
3567 10349480 : switch(typ(x))
3568 : {
3569 2891418 : case t_INT:
3570 : case t_REAL:
3571 : case t_INTMOD:
3572 : case t_FFELT:
3573 : case t_PADIC:
3574 : case t_SER:
3575 2891418 : case t_VECSMALL: return gen_1;
3576 81972 : case t_FRAC: return gel(x,2);
3577 294 : case t_COMPLEX: return vecdenom(x,1,2);
3578 69069 : case t_QUAD: return vecdenom(x,2,3);
3579 42 : case t_POLMOD: return denom_i(gel(x,2));
3580 7263520 : case t_RFRAC: return gel(x,2);
3581 973 : case t_POL: return pol_1(varn(x));
3582 42178 : case t_VEC: case t_COL: case t_MAT: return vecdenom(x, 1, lg(x)-1);
3583 : }
3584 14 : pari_err_TYPE("denom",x);
3585 : return NULL; /* LCOV_EXCL_LINE */
3586 : }
3587 : /* v has lower (or equal) priority as x's main variable */
3588 : static GEN
3589 119 : denompol(GEN x, long v)
3590 : {
3591 119 : long vx, tx = typ(x);
3592 119 : if (is_scalar_t(tx)) return gen_1;
3593 105 : switch(typ(x))
3594 : {
3595 14 : case t_SER:
3596 14 : if (varn(x) != v) return x;
3597 14 : vx = valser(x); return vx < 0? pol_xn(-vx, v): pol_1(v);
3598 63 : case t_RFRAC: x = gel(x,2); return varn(x) == v? x: pol_1(v);
3599 14 : case t_POL: return pol_1(v);
3600 14 : case t_VEC: case t_COL: case t_MAT: return vecdenompol(x, 1, lg(x)-1, v);
3601 : }
3602 0 : pari_err_TYPE("denom",x);
3603 : return NULL; /* LCOV_EXCL_LINE */
3604 : }
3605 : GEN
3606 225024 : denom(GEN x) { pari_sp av = avma; return gerepilecopy(av, denom_i(x)); }
3607 :
3608 : static GEN
3609 287 : denominator_v(GEN x, long v)
3610 : {
3611 287 : long v0 = gvar(x);
3612 : GEN d;
3613 287 : if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
3614 105 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3615 105 : d = denompol(x, v0);
3616 105 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3617 105 : return d;
3618 : }
3619 : GEN
3620 124264 : denominator(GEN x, GEN D)
3621 : {
3622 124264 : pari_sp av = avma;
3623 : GEN d;
3624 124264 : if (!D) return denom(x);
3625 280 : if (isint1(D))
3626 : {
3627 140 : d = Q_denom_safe(x);
3628 140 : if (!d) { set_avma(av); return gen_1; }
3629 105 : return gerepilecopy(av, d);
3630 : }
3631 140 : if (!gequalX(D)) pari_err_TYPE("denominator", D);
3632 140 : return gerepileupto(av, denominator_v(x, varn(D)));
3633 : }
3634 : GEN
3635 8925 : numerator(GEN x, GEN D)
3636 : {
3637 8925 : pari_sp av = avma;
3638 : long v;
3639 8925 : if (!D) return numer(x);
3640 294 : if (isint1(D)) return Q_remove_denom(x,NULL);
3641 154 : if (!gequalX(D)) pari_err_TYPE("numerator", D);
3642 154 : v = varn(D); /* optimization */
3643 154 : if (typ(x) == t_RFRAC && varn(gel(x,2)) == v) return gcopy(gel(x,1));
3644 147 : return gerepileupto(av, gmul(x, denominator_v(x,v)));
3645 : }
3646 : GEN
3647 131005 : content0(GEN x, GEN D)
3648 : {
3649 131005 : pari_sp av = avma;
3650 : long v, v0;
3651 : GEN d;
3652 131005 : if (!D) return content(x);
3653 294 : if (isint1(D))
3654 : {
3655 140 : d = Q_content_safe(x);
3656 140 : return d? d: gen_1;
3657 : }
3658 154 : if (!gequalX(D)) pari_err_TYPE("content", D);
3659 154 : v = varn(D);
3660 154 : v0 = gvar(x); if (v0 == NO_VARIABLE) return gen_1;
3661 56 : if (varncmp(v0,v) > 0)
3662 : {
3663 0 : v0 = gvar2(x);
3664 0 : return v0 == NO_VARIABLE? gen_1: pol_1(v0);
3665 : }
3666 56 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3667 56 : d = content(x);
3668 : /* gsubst is needed because of content([x]) = x */
3669 56 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3670 56 : return gerepileupto(av, d);
3671 : }
3672 :
3673 : GEN
3674 8977687 : numer_i(GEN x)
3675 : {
3676 8977687 : switch(typ(x))
3677 : {
3678 1715609 : case t_INT:
3679 : case t_REAL:
3680 : case t_INTMOD:
3681 : case t_FFELT:
3682 : case t_PADIC:
3683 : case t_SER:
3684 : case t_VECSMALL:
3685 1715609 : case t_POL: return x;
3686 28 : case t_POLMOD: return mkpolmod(numer_i(gel(x,2)), gel(x,1));
3687 7261861 : case t_FRAC:
3688 7261861 : case t_RFRAC: return gel(x,1);
3689 175 : case t_COMPLEX:
3690 : case t_QUAD:
3691 : case t_VEC:
3692 : case t_COL:
3693 175 : case t_MAT: return gmul(denom_i(x),x);
3694 : }
3695 14 : pari_err_TYPE("numer",x);
3696 : return NULL; /* LCOV_EXCL_LINE */
3697 : }
3698 : GEN
3699 8631 : numer(GEN x) { pari_sp av = avma; return gerepilecopy(av, numer_i(x)); }
3700 :
3701 : /* Lift only intmods if v does not occur in x, lift with respect to main
3702 : * variable of x if v < 0, with respect to variable v otherwise */
3703 : GEN
3704 2306975 : lift0(GEN x, long v)
3705 : {
3706 : GEN y;
3707 :
3708 2306975 : switch(typ(x))
3709 : {
3710 30380 : case t_INT: return icopy(x);
3711 2173324 : case t_INTMOD: return v < 0? icopy(gel(x,2)): gcopy(x);
3712 91742 : case t_POLMOD:
3713 91742 : if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
3714 14 : y = cgetg(3, t_POLMOD);
3715 14 : gel(y,1) = lift0(gel(x,1),v);
3716 14 : gel(y,2) = lift0(gel(x,2),v); return y;
3717 665 : case t_PADIC: return v < 0? padic_to_Q(x): gcopy(x);
3718 8715 : case t_POL:
3719 41405 : pari_APPLY_pol(lift0(gel(x,i), v));
3720 56 : case t_SER:
3721 56 : if (ser_isexactzero(x))
3722 : {
3723 14 : if (lg(x) == 2) return gcopy(x);
3724 14 : y = scalarser(lift0(gel(x,2),v), varn(x), 1);
3725 14 : setvalser(y, valser(x)); return y;
3726 : }
3727 434 : pari_APPLY_ser(lift0(gel(x,i), v));
3728 1904 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3729 : case t_VEC: case t_COL: case t_MAT:
3730 41139 : pari_APPLY_same(lift0(gel(x,i), v));
3731 189 : default: return gcopy(x);
3732 : }
3733 : }
3734 : /* same as lift, shallow */
3735 : GEN
3736 592911 : lift_shallow(GEN x)
3737 : {
3738 : GEN y;
3739 592911 : switch(typ(x))
3740 : {
3741 197140 : case t_INTMOD: case t_POLMOD: return gel(x,2);
3742 476 : case t_PADIC: return padic_to_Q(x);
3743 0 : case t_SER:
3744 0 : if (ser_isexactzero(x))
3745 : {
3746 0 : if (lg(x) == 2) return x;
3747 0 : y = scalarser(lift_shallow(gel(x,2)), varn(x), 1);
3748 0 : setvalser(y, valser(x)); return y;
3749 : }
3750 0 : pari_APPLY_ser(lift_shallow(gel(x,i)));
3751 47852 : case t_POL:
3752 272703 : pari_APPLY_pol(lift_shallow(gel(x,i)));
3753 11158 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3754 : case t_VEC: case t_COL: case t_MAT:
3755 275597 : pari_APPLY_same(lift_shallow(gel(x,i)));
3756 336285 : default: return x;
3757 : }
3758 : }
3759 : GEN
3760 2156092 : lift(GEN x) { return lift0(x,-1); }
3761 :
3762 : GEN
3763 2003435 : liftall_shallow(GEN x)
3764 : {
3765 : GEN y;
3766 2003435 : switch(typ(x))
3767 : {
3768 533778 : case t_INTMOD: return gel(x,2);
3769 547519 : case t_POLMOD:
3770 547519 : return liftall_shallow(gel(x,2));
3771 581 : case t_PADIC: return padic_to_Q(x);
3772 555912 : case t_POL:
3773 1356460 : pari_APPLY_pol(liftall_shallow(gel(x,i)));
3774 7 : case t_SER:
3775 7 : if (ser_isexactzero(x))
3776 : {
3777 0 : if (lg(x) == 2) return x;
3778 0 : y = scalarser(liftall_shallow(gel(x,2)), varn(x), 1);
3779 0 : setvalser(y, valser(x)); return y;
3780 : }
3781 35 : pari_APPLY_ser(liftall_shallow(gel(x,i)));
3782 132762 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3783 : case t_VEC: case t_COL: case t_MAT:
3784 760515 : pari_APPLY_same(liftall_shallow(gel(x,i)));
3785 232876 : default: return x;
3786 : }
3787 : }
3788 : GEN
3789 26243 : liftall(GEN x)
3790 26243 : { pari_sp av = avma; return gerepilecopy(av, liftall_shallow(x)); }
3791 :
3792 : GEN
3793 546 : liftint_shallow(GEN x)
3794 : {
3795 : GEN y;
3796 546 : switch(typ(x))
3797 : {
3798 266 : case t_INTMOD: return gel(x,2);
3799 28 : case t_PADIC: return padic_to_Q(x);
3800 21 : case t_POL:
3801 70 : pari_APPLY_pol(liftint_shallow(gel(x,i)));
3802 14 : case t_SER:
3803 14 : if (ser_isexactzero(x))
3804 : {
3805 7 : if (lg(x) == 2) return x;
3806 7 : y = scalarser(liftint_shallow(gel(x,2)), varn(x), 1);
3807 7 : setvalser(y, valser(x)); return y;
3808 : }
3809 35 : pari_APPLY_ser(liftint_shallow(gel(x,i)));
3810 161 : case t_POLMOD: case t_COMPLEX: case t_QUAD: case t_RFRAC:
3811 : case t_VEC: case t_COL: case t_MAT:
3812 504 : pari_APPLY_same(liftint_shallow(gel(x,i)));
3813 56 : default: return x;
3814 : }
3815 : }
3816 : GEN
3817 119 : liftint(GEN x)
3818 119 : { pari_sp av = avma; return gerepilecopy(av, liftint_shallow(x)); }
3819 :
3820 : GEN
3821 21236006 : liftpol_shallow(GEN x)
3822 : {
3823 : GEN y;
3824 21236006 : switch(typ(x))
3825 : {
3826 2048352 : case t_POLMOD:
3827 2048352 : return liftpol_shallow(gel(x,2));
3828 2828794 : case t_POL:
3829 11759477 : pari_APPLY_pol(liftpol_shallow(gel(x,i)));
3830 7 : case t_SER:
3831 7 : if (ser_isexactzero(x))
3832 : {
3833 0 : if (lg(x) == 2) return x;
3834 0 : y = scalarser(liftpol(gel(x,2)), varn(x), 1);
3835 0 : setvalser(y, valser(x)); return y;
3836 : }
3837 35 : pari_APPLY_ser(liftpol_shallow(gel(x,i)));
3838 132013 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
3839 946204 : pari_APPLY_same(liftpol_shallow(gel(x,i)));
3840 16226840 : default: return x;
3841 : }
3842 : }
3843 : GEN
3844 5656 : liftpol(GEN x)
3845 5656 : { pari_sp av = avma; return gerepilecopy(av, liftpol_shallow(x)); }
3846 :
3847 : static GEN
3848 42518 : centerliftii(GEN x, GEN y)
3849 : {
3850 42518 : pari_sp av = avma;
3851 42518 : long i = cmpii(shifti(x,1), y);
3852 42518 : set_avma(av); return (i > 0)? subii(x,y): icopy(x);
3853 : }
3854 :
3855 : /* see lift0 */
3856 : GEN
3857 707 : centerlift0(GEN x, long v)
3858 707 : { return v < 0? centerlift(x): lift0(x,v); }
3859 : GEN
3860 60473 : centerlift(GEN x)
3861 : {
3862 : long v;
3863 : GEN y;
3864 60473 : switch(typ(x))
3865 : {
3866 784 : case t_INT: return icopy(x);
3867 784 : case t_INTMOD: return centerliftii(gel(x,2), gel(x,1));
3868 7 : case t_POLMOD: return gcopy(gel(x,2));
3869 1554 : case t_POL:
3870 9912 : pari_APPLY_pol(centerlift(gel(x,i)));
3871 7 : case t_SER:
3872 7 : if (ser_isexactzero(x)) return lift(x);
3873 35 : pari_APPLY_ser(centerlift(gel(x,i)));
3874 5551 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3875 : case t_VEC: case t_COL: case t_MAT:
3876 56210 : pari_APPLY_same(centerlift(gel(x,i)));
3877 51779 : case t_PADIC:
3878 51779 : if (!signe(gel(x,4))) return gen_0;
3879 41734 : v = valp(x);
3880 41734 : if (v>=0)
3881 : { /* here p^v is an integer */
3882 41727 : GEN z = centerliftii(gel(x,4), gel(x,3));
3883 : pari_sp av;
3884 41727 : if (!v) return z;
3885 27027 : av = avma; y = powiu(gel(x,2),v);
3886 27027 : return gerepileuptoint(av, mulii(y,z));
3887 : }
3888 7 : y = cgetg(3,t_FRAC);
3889 7 : gel(y,1) = centerliftii(gel(x,4), gel(x,3));
3890 7 : gel(y,2) = powiu(gel(x,2),-v);
3891 7 : return y;
3892 7 : default: return gcopy(x);
3893 : }
3894 : }
3895 :
3896 : /*******************************************************************/
3897 : /* */
3898 : /* REAL & IMAGINARY PARTS */
3899 : /* */
3900 : /*******************************************************************/
3901 :
3902 : static GEN
3903 202095578 : op_ReIm(GEN f(GEN), GEN x)
3904 : {
3905 202095578 : switch(typ(x))
3906 : {
3907 589631365 : case t_POL: pari_APPLY_pol(f(gel(x,i)));
3908 64463 : case t_SER: pari_APPLY_ser(f(gel(x,i)));
3909 42 : case t_RFRAC: {
3910 42 : pari_sp av = avma;
3911 42 : GEN n, d, dxb = conj_i(gel(x,2));
3912 42 : n = gmul(gel(x,1), dxb);
3913 42 : d = gmul(gel(x,2), dxb);
3914 42 : return gerepileupto(av, gdiv(f(n), d));
3915 : }
3916 19289240 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(f(gel(x, i)));
3917 : }
3918 12 : pari_err_TYPE("greal/gimag",x);
3919 : return NULL; /* LCOV_EXCL_LINE */
3920 : }
3921 :
3922 : GEN
3923 318951648 : real_i(GEN x)
3924 : {
3925 318951648 : switch(typ(x))
3926 : {
3927 171244542 : case t_INT: case t_REAL: case t_FRAC:
3928 171244542 : return x;
3929 44109003 : case t_COMPLEX:
3930 44109003 : return gel(x,1);
3931 0 : case t_QUAD:
3932 0 : return gel(x,2);
3933 : }
3934 103598103 : return op_ReIm(real_i,x);
3935 : }
3936 : GEN
3937 297089965 : imag_i(GEN x)
3938 : {
3939 297089965 : switch(typ(x))
3940 : {
3941 165904864 : case t_INT: case t_REAL: case t_FRAC:
3942 165904864 : return gen_0;
3943 32750124 : case t_COMPLEX:
3944 32750124 : return gel(x,2);
3945 0 : case t_QUAD:
3946 0 : return gel(x,3);
3947 : }
3948 98434977 : return op_ReIm(imag_i,x);
3949 : }
3950 : GEN
3951 7014 : greal(GEN x)
3952 : {
3953 7014 : switch(typ(x))
3954 : {
3955 707 : case t_INT: case t_REAL:
3956 707 : return mpcopy(x);
3957 :
3958 7 : case t_FRAC:
3959 7 : return gcopy(x);
3960 :
3961 6055 : case t_COMPLEX:
3962 6055 : return gcopy(gel(x,1));
3963 :
3964 7 : case t_QUAD:
3965 7 : return gcopy(gel(x,2));
3966 : }
3967 238 : return op_ReIm(greal,x);
3968 : }
3969 : GEN
3970 29344 : gimag(GEN x)
3971 : {
3972 29344 : switch(typ(x))
3973 : {
3974 525 : case t_INT: case t_REAL: case t_FRAC:
3975 525 : return gen_0;
3976 :
3977 28203 : case t_COMPLEX:
3978 28203 : return gcopy(gel(x,2));
3979 :
3980 7 : case t_QUAD:
3981 7 : return gcopy(gel(x,3));
3982 : }
3983 609 : return op_ReIm(gimag,x);
3984 : }
3985 :
3986 : /* return Re(x * y), x and y scalars */
3987 : GEN
3988 15684708 : mulreal(GEN x, GEN y)
3989 : {
3990 15684708 : if (typ(x) == t_COMPLEX)
3991 : {
3992 15541047 : if (typ(y) == t_COMPLEX)
3993 : {
3994 14322546 : pari_sp av = avma;
3995 14322546 : GEN a = gmul(gel(x,1), gel(y,1));
3996 14322541 : GEN b = gmul(gel(x,2), gel(y,2));
3997 14322545 : return gerepileupto(av, gsub(a, b));
3998 : }
3999 1218501 : x = gel(x,1);
4000 : }
4001 : else
4002 143661 : if (typ(y) == t_COMPLEX) y = gel(y,1);
4003 1362162 : return gmul(x,y);
4004 : }
4005 : /* Compute Re(x * y), x and y matrices of compatible dimensions
4006 : * assume scalar entries */
4007 : GEN
4008 0 : RgM_mulreal(GEN x, GEN y)
4009 : {
4010 0 : long i, j, k, l, lx = lg(x), ly = lg(y);
4011 0 : GEN z = cgetg(ly,t_MAT);
4012 0 : l = (lx == 1)? 1: lgcols(x);
4013 0 : for (j=1; j<ly; j++)
4014 : {
4015 0 : GEN zj = cgetg(l,t_COL), yj = gel(y,j);
4016 0 : gel(z,j) = zj;
4017 0 : for (i=1; i<l; i++)
4018 : {
4019 0 : pari_sp av = avma;
4020 0 : GEN c = mulreal(gcoeff(x,i,1),gel(yj,1));
4021 0 : for (k=2; k<lx; k++) c = gadd(c, mulreal(gcoeff(x,i,k),gel(yj,k)));
4022 0 : gel(zj, i) = gerepileupto(av, c);
4023 : }
4024 : }
4025 0 : return z;
4026 : }
4027 :
4028 : /* Compute Re(x * y), symmetric result, x and y vectors of compatible
4029 : * dimensions; assume scalar entries */
4030 : GEN
4031 21315 : RgC_RgV_mulrealsym(GEN x, GEN y)
4032 : {
4033 21315 : long i, j, l = lg(x);
4034 21315 : GEN q = cgetg(l, t_MAT);
4035 85260 : for (j = 1; j < l; j++)
4036 : {
4037 63945 : gel(q,j) = cgetg(l,t_COL);
4038 191835 : for (i = 1; i <= j; i++)
4039 127890 : gcoeff(q,i,j) = gcoeff(q,j,i) = mulreal(gel(x,i), gel(y,j));
4040 : }
4041 21315 : return q;
4042 : }
4043 :
4044 : /*******************************************************************/
4045 : /* */
4046 : /* LOGICAL OPERATIONS */
4047 : /* */
4048 : /*******************************************************************/
4049 : static long
4050 108043681 : _egal_i(GEN x, GEN y)
4051 : {
4052 108043681 : x = simplify_shallow(x);
4053 108043684 : y = simplify_shallow(y);
4054 108043683 : if (typ(y) == t_INT)
4055 : {
4056 107056501 : if (equali1(y)) return gequal1(x);
4057 61892603 : if (equalim1(y)) return gequalm1(x);
4058 : }
4059 987182 : else if (typ(x) == t_INT)
4060 : {
4061 133 : if (equali1(x)) return gequal1(y);
4062 84 : if (equalim1(x)) return gequalm1(y);
4063 : }
4064 62748164 : return gequal(x, y);
4065 : }
4066 : static long
4067 108043681 : _egal(GEN x, GEN y)
4068 108043681 : { pari_sp av = avma; return gc_long(av, _egal_i(x, y)); }
4069 :
4070 : GEN
4071 6328186 : glt(GEN x, GEN y) { return gcmp(x,y)<0? gen_1: gen_0; }
4072 :
4073 : GEN
4074 7628237 : gle(GEN x, GEN y) { return gcmp(x,y)<=0? gen_1: gen_0; }
4075 :
4076 : GEN
4077 248429 : gge(GEN x, GEN y) { return gcmp(x,y)>=0? gen_1: gen_0; }
4078 :
4079 : GEN
4080 2374836 : ggt(GEN x, GEN y) { return gcmp(x,y)>0? gen_1: gen_0; }
4081 :
4082 : GEN
4083 47185640 : geq(GEN x, GEN y) { return _egal(x,y)? gen_1: gen_0; }
4084 :
4085 : GEN
4086 60858041 : gne(GEN x, GEN y) { return _egal(x,y)? gen_0: gen_1; }
4087 :
4088 : GEN
4089 603995 : gnot(GEN x) { return gequal0(x)? gen_1: gen_0; }
4090 :
4091 : /*******************************************************************/
4092 : /* */
4093 : /* FORMAL SIMPLIFICATIONS */
4094 : /* */
4095 : /*******************************************************************/
4096 :
4097 : GEN
4098 11020 : geval_gp(GEN x, GEN t)
4099 : {
4100 11020 : long lx, i, tx = typ(x);
4101 : pari_sp av;
4102 : GEN y, z;
4103 :
4104 11020 : if (is_const_t(tx) || tx==t_VECSMALL) return gcopy(x);
4105 10999 : switch(tx)
4106 : {
4107 10992 : case t_STR:
4108 10992 : return localvars_read_str(GSTR(x),t);
4109 :
4110 0 : case t_POLMOD:
4111 0 : av = avma;
4112 0 : return gerepileupto(av, gmodulo(geval_gp(gel(x,2),t),
4113 0 : geval_gp(gel(x,1),t)));
4114 :
4115 7 : case t_POL:
4116 7 : lx=lg(x); if (lx==2) return gen_0;
4117 7 : z = fetch_var_value(varn(x),t);
4118 7 : if (!z) return RgX_copy(x);
4119 7 : av = avma; y = geval_gp(gel(x,lx-1),t);
4120 14 : for (i=lx-2; i>1; i--)
4121 7 : y = gadd(geval_gp(gel(x,i),t), gmul(z,y));
4122 7 : return gerepileupto(av, y);
4123 :
4124 0 : case t_SER:
4125 0 : pari_err_IMPL( "evaluation of a power series");
4126 :
4127 0 : case t_RFRAC:
4128 0 : av = avma;
4129 0 : return gerepileupto(av, gdiv(geval_gp(gel(x,1),t), geval_gp(gel(x,2),t)));
4130 :
4131 0 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
4132 0 : pari_APPLY_same(geval_gp(gel(x,i),t));
4133 :
4134 0 : case t_CLOSURE:
4135 0 : if (closure_arity(x) || closure_is_variadic(x))
4136 0 : pari_err_IMPL("eval on functions with parameters");
4137 0 : return closure_evalres(x);
4138 : }
4139 0 : pari_err_TYPE("geval",x);
4140 : return NULL; /* LCOV_EXCL_LINE */
4141 : }
4142 : GEN
4143 0 : geval(GEN x) { return geval_gp(x,NULL); }
4144 :
4145 : GEN
4146 527279240 : simplify_shallow(GEN x)
4147 : {
4148 : long v, lx;
4149 : GEN y, z;
4150 527279240 : if (!x) pari_err_BUG("simplify, NULL input");
4151 :
4152 527279240 : switch(typ(x))
4153 : {
4154 444343785 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
4155 : case t_PADIC: case t_QFB: case t_LIST: case t_STR: case t_VECSMALL:
4156 : case t_CLOSURE: case t_ERROR: case t_INFINITY:
4157 444343785 : return x;
4158 735534 : case t_COMPLEX: return isintzero(gel(x,2))? gel(x,1): x;
4159 805 : case t_QUAD: return isintzero(gel(x,3))? gel(x,2): x;
4160 :
4161 197222 : case t_POLMOD: y = cgetg(3,t_POLMOD);
4162 197222 : z = gel(x,1); v = varn(z); z = simplify_shallow(z);
4163 197222 : if (typ(z) != t_POL || varn(z) != v) z = scalarpol_shallow(z, v);
4164 197222 : gel(y,1) = z;
4165 197222 : gel(y,2) = simplify_shallow(gel(x,2)); return y;
4166 :
4167 68997603 : case t_POL:
4168 68997603 : lx = lg(x);
4169 68997603 : if (lx==2) return gen_0;
4170 61412405 : if (lx==3) return simplify_shallow(gel(x,2));
4171 198146054 : pari_APPLY_pol(simplify_shallow(gel(x,i)));
4172 :
4173 3633 : case t_SER:
4174 1065337 : pari_APPLY_ser(simplify_shallow(gel(x,i)));
4175 :
4176 651691 : case t_RFRAC: y = cgetg(3,t_RFRAC);
4177 651691 : gel(y,1) = simplify_shallow(gel(x,1));
4178 651691 : z = simplify_shallow(gel(x,2));
4179 651691 : if (typ(z) != t_POL) return gdiv(gel(y,1), z);
4180 651691 : gel(y,2) = z; return y;
4181 :
4182 12348967 : case t_VEC: case t_COL: case t_MAT:
4183 68400592 : pari_APPLY_same(simplify_shallow(gel(x,i)));
4184 : }
4185 0 : pari_err_BUG("simplify_shallow, type unknown");
4186 : return NULL; /* LCOV_EXCL_LINE */
4187 : }
4188 :
4189 : GEN
4190 10810060 : simplify(GEN x)
4191 : {
4192 10810060 : pari_sp av = avma;
4193 10810060 : GEN y = simplify_shallow(x);
4194 10810060 : return av == avma ? gcopy(y): gerepilecopy(av, y);
4195 : }
4196 :
4197 : /*******************************************************************/
4198 : /* */
4199 : /* EVALUATION OF SOME SIMPLE OBJECTS */
4200 : /* */
4201 : /*******************************************************************/
4202 : /* q is a real symmetric matrix, x a RgV. Horner-type evaluation of q(x)
4203 : * using (n^2+3n-2)/2 mul */
4204 : GEN
4205 16848 : qfeval(GEN q, GEN x)
4206 : {
4207 16848 : pari_sp av = avma;
4208 16848 : long i, j, l = lg(q);
4209 : GEN z;
4210 16848 : if (lg(x) != l) pari_err_DIM("qfeval");
4211 16841 : if (l==1) return gen_0;
4212 16841 : if (lgcols(q) != l) pari_err_DIM("qfeval");
4213 : /* l = lg(x) = lg(q) > 1 */
4214 16834 : z = gmul(gcoeff(q,1,1), gsqr(gel(x,1)));
4215 73490 : for (i=2; i<l; i++)
4216 : {
4217 56656 : GEN c = gel(q,i), s;
4218 56656 : if (isintzero(gel(x,i))) continue;
4219 43006 : s = gmul(gel(c,1), gel(x,1));
4220 131285 : for (j=2; j<i; j++) s = gadd(s, gmul(gel(c,j),gel(x,j)));
4221 43006 : s = gadd(gshift(s,1), gmul(gel(c,i),gel(x,i)));
4222 43006 : z = gadd(z, gmul(gel(x,i), s));
4223 : }
4224 16834 : return gerepileupto(av,z);
4225 : }
4226 :
4227 : static GEN
4228 347816 : qfbeval(GEN q, GEN z)
4229 : {
4230 347816 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3), x = gel(z,1), y = gel(z,2);
4231 347816 : pari_sp av = avma;
4232 347816 : A = gadd(gmul(x, gadd(gmul(a,x), gmul(b,y))), gmul(c, gsqr(y)));
4233 347816 : return gerepileupto(av, A);
4234 : }
4235 : static GEN
4236 7 : qfbevalb(GEN q, GEN z, GEN z2)
4237 : {
4238 7 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3);
4239 7 : GEN x = gel(z,1), y = gel(z,2);
4240 7 : GEN X = gel(z2,1), Y = gel(z2,2);
4241 7 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4242 7 : pari_sp av = avma;
4243 : /* a2 x X + b (x Y + X y) + c2 y Y */
4244 7 : A = gadd(gmul(x, gadd(gmul(a2,X), gmul(b,Y))),
4245 : gmul(y, gadd(gmul(c2,Y), gmul(b,X))));
4246 7 : return gerepileupto(av, gmul2n(A, -1));
4247 : }
4248 : GEN
4249 61920 : qfb_ZM_apply(GEN q, GEN M)
4250 : {
4251 61920 : pari_sp av = avma;
4252 61920 : GEN A, B, C, a = gel(q,1), b = gel(q,2), c = gel(q,3);
4253 61920 : GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
4254 61920 : GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
4255 61920 : GEN by = mulii(b,y), bt = mulii(b,t), bz = mulii(b,z);
4256 61920 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4257 :
4258 61920 : A = addii(mulii(x, addii(mulii(a,x), by)), mulii(c, sqri(y)));
4259 61920 : B = addii(mulii(x, addii(mulii(a2,z), bt)),
4260 : mulii(y, addii(mulii(c2,t), bz)));
4261 61920 : C = addii(mulii(z, addii(mulii(a,z), bt)), mulii(c, sqri(t)));
4262 61920 : q = leafcopy(q);
4263 61920 : gel(q,1) = A; gel(q,2) = B; gel(q,3) = C;
4264 61920 : return gerepilecopy(av, q);
4265 : }
4266 :
4267 : static GEN
4268 348012 : qfnorm0(GEN q, GEN x)
4269 : {
4270 348012 : if (!q) switch(typ(x))
4271 : {
4272 7 : case t_VEC: case t_COL: return RgV_dotsquare(x);
4273 7 : case t_MAT: return gram_matrix(x);
4274 7 : default: pari_err_TYPE("qfeval",x);
4275 : }
4276 347991 : switch(typ(q))
4277 : {
4278 161 : case t_MAT: break;
4279 347823 : case t_QFB:
4280 347823 : if (lg(x) == 3) switch(typ(x))
4281 : {
4282 347816 : case t_VEC:
4283 347816 : case t_COL: return qfbeval(q,x);
4284 7 : case t_MAT: if (RgM_is_ZM(x)) return qfb_ZM_apply(q,x);
4285 0 : default: pari_err_TYPE("qfeval",x);
4286 : }
4287 7 : default: pari_err_TYPE("qfeval",q);
4288 : }
4289 161 : switch(typ(x))
4290 : {
4291 154 : case t_VEC: case t_COL: break;
4292 7 : case t_MAT: return qf_RgM_apply(q, x);
4293 0 : default: pari_err_TYPE("qfeval",x);
4294 : }
4295 154 : return qfeval(q,x);
4296 : }
4297 : /* obsolete, use qfeval0 */
4298 : GEN
4299 0 : qfnorm(GEN x, GEN q) { return qfnorm0(q,x); }
4300 :
4301 : /* assume q is square, x~ * q * y using n^2+n mul */
4302 : GEN
4303 21 : qfevalb(GEN q, GEN x, GEN y)
4304 : {
4305 21 : pari_sp av = avma;
4306 21 : long l = lg(q);
4307 21 : if (lg(x) != l || lg(y) != l) pari_err_DIM("qfevalb");
4308 14 : return gerepileupto(av, RgV_dotproduct(RgV_RgM_mul(x,q), y));
4309 : }
4310 :
4311 : /* obsolete, use qfeval0 */
4312 : GEN
4313 0 : qfbil(GEN x, GEN y, GEN q)
4314 : {
4315 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfbil",x);
4316 0 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfbil",y);
4317 0 : if (!q) {
4318 0 : if (lg(x) != lg(y)) pari_err_DIM("qfbil");
4319 0 : return RgV_dotproduct(x,y);
4320 : }
4321 0 : if (typ(q) != t_MAT) pari_err_TYPE("qfbil",q);
4322 0 : return qfevalb(q,x,y);
4323 : }
4324 : GEN
4325 348068 : qfeval0(GEN q, GEN x, GEN y)
4326 : {
4327 348068 : if (!y) return qfnorm0(q, x);
4328 56 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfeval",x);
4329 42 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfeval",y);
4330 42 : if (!q) {
4331 14 : if (lg(x) != lg(y)) pari_err_DIM("qfeval");
4332 7 : return RgV_dotproduct(x,y);
4333 : }
4334 28 : switch(typ(q))
4335 : {
4336 21 : case t_MAT: break;
4337 7 : case t_QFB:
4338 7 : if (lg(x) == 3 && lg(y) == 3) return qfbevalb(q,x,y);
4339 0 : default: pari_err_TYPE("qfeval",q);
4340 : }
4341 21 : return qfevalb(q,x,y);
4342 : }
4343 :
4344 : /* q a hermitian complex matrix, x a RgV */
4345 : GEN
4346 0 : hqfeval(GEN q, GEN x)
4347 : {
4348 0 : pari_sp av = avma;
4349 0 : long i, j, l = lg(q);
4350 : GEN z, xc;
4351 :
4352 0 : if (lg(x) != l) pari_err_DIM("hqfeval");
4353 0 : if (l==1) return gen_0;
4354 0 : if (lgcols(q) != l) pari_err_DIM("hqfeval");
4355 0 : if (l == 2) return gerepileupto(av, gmul(gcoeff(q,1,1), gnorm(gel(x,1))));
4356 : /* l = lg(x) = lg(q) > 2 */
4357 0 : xc = conj_i(x);
4358 0 : z = mulreal(gcoeff(q,2,1), gmul(gel(x,2),gel(xc,1)));
4359 0 : for (i=3;i<l;i++)
4360 0 : for (j=1;j<i;j++)
4361 0 : z = gadd(z, mulreal(gcoeff(q,i,j), gmul(gel(x,i),gel(xc,j))));
4362 0 : z = gshift(z,1);
4363 0 : for (i=1;i<l;i++) z = gadd(z, gmul(gcoeff(q,i,i), gnorm(gel(x,i))));
4364 0 : return gerepileupto(av,z);
4365 : }
4366 :
4367 : static void
4368 500579 : init_qf_apply(GEN q, GEN M, long *l)
4369 : {
4370 500579 : long k = lg(M);
4371 500579 : *l = lg(q);
4372 500579 : if (*l == 1) { if (k == 1) return; }
4373 500572 : else { if (k != 1 && lgcols(M) == *l) return; }
4374 0 : pari_err_DIM("qf_RgM_apply");
4375 : }
4376 : /* Return X = M'.q.M, assuming q is a symmetric matrix and M is a
4377 : * matrix of compatible dimensions. X_ij are X_ji identical, not copies */
4378 : GEN
4379 7570 : qf_RgM_apply(GEN q, GEN M)
4380 : {
4381 7570 : pari_sp av = avma;
4382 7570 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4383 7570 : return gerepileupto(av, RgM_transmultosym(M, RgM_mul(q, M)));
4384 : }
4385 : GEN
4386 493009 : qf_ZM_apply(GEN q, GEN M)
4387 : {
4388 493009 : pari_sp av = avma;
4389 493009 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4390 : /* FIXME: ZM_transmultosym is asymptotically slow, so choose some random
4391 : * threshold defaulting to default implementation: this is suboptimal but
4392 : * has the right complexity in the dimension. Should implement M'*q*M as an
4393 : * atomic operation with the right complexity, see ZM_mul_i. */
4394 493002 : if (l > 20)
4395 161 : M = ZM_mul(shallowtrans(M), ZM_mul(q, M));
4396 : else
4397 492841 : M = ZM_transmultosym(M, ZM_mul(q, M));
4398 493002 : return gerepileupto(av, M);
4399 : }
4400 :
4401 : GEN
4402 2907068 : poleval(GEN x, GEN y)
4403 : {
4404 2907068 : long i, j, imin, tx = typ(x);
4405 2907068 : pari_sp av0 = avma, av;
4406 : GEN t, t2, r, s;
4407 :
4408 2907068 : if (is_scalar_t(tx)) return gcopy(x);
4409 2725770 : switch(tx)
4410 : {
4411 2723845 : case t_POL:
4412 2723845 : i = lg(x)-1; imin = 2; break;
4413 :
4414 1568 : case t_RFRAC:
4415 1568 : return gerepileupto(av0, gdiv(poleval(gel(x,1),y),
4416 1568 : poleval(gel(x,2),y)));
4417 :
4418 357 : case t_VEC: case t_COL:
4419 357 : i = lg(x)-1; imin = 1; break;
4420 0 : default: pari_err_TYPE("poleval",x);
4421 : return NULL; /* LCOV_EXCL_LINE */
4422 : }
4423 2724202 : if (i<=imin) return (i==imin)? gmul(gel(x,imin),Rg_get_1(y)): Rg_get_0(y);
4424 2576933 : if (isintzero(y)) return gcopy(gel(x,imin));
4425 :
4426 2569785 : t = gel(x,i); i--;
4427 2569785 : if (typ(y)!=t_COMPLEX)
4428 : {
4429 : #if 0 /* standard Horner's rule */
4430 : for ( ; i>=imin; i--)
4431 : t = gadd(gmul(t,y),gel(x,i));
4432 : #endif
4433 : /* specific attention to sparse polynomials */
4434 18646043 : for ( ; i>=imin; i=j-1)
4435 : {
4436 18859640 : for (j=i; isexactzero(gel(x,j)); j--)
4437 2685760 : if (j==imin)
4438 : {
4439 970573 : if (i!=j) y = gpowgs(y, i-j+1);
4440 970573 : return gerepileupto(av0, gmul(t,y));
4441 : }
4442 16173884 : r = (i==j)? y: gpowgs(y, i-j+1);
4443 16173884 : t = gadd(gmul(t,r), gel(x,j));
4444 16173844 : if (gc_needed(av0,2))
4445 : {
4446 103 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4447 103 : t = gerepileupto(av0, t);
4448 : }
4449 : }
4450 1501590 : return gerepileupto(av0, t);
4451 : }
4452 :
4453 97586 : t2 = gel(x,i); i--; r = gtrace(y); s = gneg_i(gnorm(y));
4454 97586 : av = avma;
4455 4957695 : for ( ; i>=imin; i--)
4456 : {
4457 4860109 : GEN p = gadd(t2, gmul(r, t));
4458 4860109 : t2 = gadd(gel(x,i), gmul(s, t)); t = p;
4459 4860109 : if (gc_needed(av0,2))
4460 : {
4461 0 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4462 0 : gerepileall(av, 2, &t, &t2);
4463 : }
4464 : }
4465 97586 : return gerepileupto(av0, gadd(t2, gmul(y, t)));
4466 : }
4467 :
4468 : /* Evaluate a polynomial using Horner. Unstable!
4469 : * If ui != NULL, ui = 1/u, evaluate P(1/u)*u^(deg P): useful for |u|>1 */
4470 : GEN
4471 2398268 : RgX_cxeval(GEN T, GEN u, GEN ui)
4472 : {
4473 2398268 : pari_sp ltop = avma;
4474 : GEN S;
4475 2398268 : long n, lim = lg(T)-1;
4476 2398268 : if (lim == 1) return gen_0;
4477 2398268 : if (lim == 2) return gcopy(gel(T,2));
4478 2397110 : if (!ui)
4479 : {
4480 1653590 : n = lim; S = gel(T,n);
4481 14877297 : for (n--; n >= 2; n--) S = gadd(gmul(u,S), gel(T,n));
4482 : }
4483 : else
4484 : {
4485 743520 : n = 2; S = gel(T,2);
4486 4468545 : for (n++; n <= lim; n++) S = gadd(gmul(ui,S), gel(T,n));
4487 743530 : S = gmul(gpowgs(u, lim-2), S);
4488 : }
4489 2396948 : return gerepileupto(ltop, S);
4490 : }
4491 :
4492 : GEN
4493 63 : RgXY_cxevalx(GEN x, GEN u, GEN ui)
4494 : {
4495 427 : pari_APPLY_pol(typ(gel(x,i))==t_POL? RgX_cxeval(gel(x,i), u, ui): gel(x,i));
4496 : }
|