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 26558 : recvar(hashtable *h, GEN x)
31 : {
32 26558 : long i = 1, lx = lg(x);
33 : void *v;
34 26558 : switch(typ(x))
35 : {
36 6083 : case t_POL: case t_SER:
37 6083 : v = (void*)varn(x);
38 6083 : if (!hash_search(h, v)) hash_insert(h, v, NULL);
39 6083 : 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 19418 : default:
46 19418 : return;
47 : }
48 32655 : for (; i < lx; i++) recvar(h, gel(x,i));
49 : }
50 :
51 : GEN
52 1043 : variables_vecsmall(GEN x)
53 : {
54 1043 : hashtable *h = hash_create_ulong(100, 1);
55 1043 : recvar(h, x);
56 1043 : return vars_sort_inplace(hash_keys(h));
57 : }
58 :
59 : GEN
60 28 : variables_vec(GEN x)
61 28 : { return x? vars_to_RgXV(variables_vecsmall(x)): gpolvar(NULL); }
62 :
63 : long
64 126731760 : gvar(GEN x)
65 : {
66 : long i, v, w, lx;
67 126731760 : switch(typ(x))
68 : {
69 48340175 : case t_POL: case t_SER: return varn(x);
70 169540 : case t_POLMOD: return varn(gel(x,1));
71 14556723 : case t_RFRAC: return varn(gel(x,2));
72 3482736 : case t_VEC: case t_COL: case t_MAT:
73 3482736 : lx = lg(x); break;
74 14 : case t_LIST:
75 14 : x = list_data(x);
76 14 : lx = x? lg(x): 1; break;
77 60182572 : default:
78 60182572 : return NO_VARIABLE;
79 : }
80 3482750 : v = NO_VARIABLE;
81 32346513 : for (i=1; i < lx; i++) { w = gvar(gel(x,i)); if (varncmp(w,v) < 0) v = w; }
82 3482755 : 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 54327 : gvar9(GEN x)
103 54327 : { 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 20637892 : gvar2(GEN x)
108 : {
109 : long i, v, w;
110 20637892 : switch(typ(x))
111 : {
112 56 : case t_POLMOD:
113 56 : return var2_polmod(x);
114 16744 : case t_POL: case t_SER:
115 16744 : v = NO_VARIABLE;
116 70049 : for (i=2; i < lg(x); i++) {
117 53305 : w = gvar9(gel(x,i));
118 53305 : if (varncmp(w,v) < 0) v=w;
119 : }
120 16744 : 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 20614008 : return NO_VARIABLE;
132 : }
133 :
134 : /*******************************************************************/
135 : /* */
136 : /* PRECISION OF SCALAR OBJECTS */
137 : /* */
138 : /*******************************************************************/
139 : static long
140 9916519 : prec0(long e) { return (e < 0)? nbits2prec(-e): 3; }
141 : static long
142 668156655 : 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 600779 : precrealexact(GEN x, GEN y)
146 : {
147 600779 : long lx, ey = gexpo(y), ex, e;
148 600816 : if (ey == -(long)HIGHEXPOBIT) return precREAL(x);
149 185840 : ex = expo(x);
150 185840 : e = ey - ex;
151 185840 : if (!signe(x)) return prec0((e >= 0)? -e: ex);
152 185791 : lx = realprec(x);
153 185791 : return (e > 0)? lx + nbits2extraprec(e): lx;
154 : }
155 : static long
156 15630618 : precCOMPLEX(GEN z)
157 : { /* ~ precision(|x| + |y|) */
158 15630618 : GEN x = gel(z,1), y = gel(z,2);
159 : long e, ex, ey, lz, lx, ly;
160 15630618 : if (typ(x) != t_REAL) {
161 1383275 : if (typ(y) != t_REAL) return 0;
162 588528 : return precrealexact(y, x);
163 : }
164 14247343 : if (typ(y) != t_REAL) return precrealexact(x, y);
165 : /* x, y are t_REALs, cf addrr_sign */
166 14235090 : ex = expo(x);
167 14235090 : ey = expo(y);
168 14235090 : e = ey - ex;
169 14235090 : if (!signe(x)) {
170 498005 : if (!signe(y)) return prec0( minss(ex,ey) );
171 497956 : if (e <= 0) return prec0(ex);
172 497746 : lz = nbits2prec(e);
173 497746 : ly = realprec(y); if (lz > ly) lz = ly;
174 497746 : return lz;
175 : }
176 13737085 : if (!signe(y)) {
177 80221 : if (e >= 0) return prec0(ey);
178 80221 : lz = nbits2prec(-e);
179 80221 : lx = realprec(x); if (lz > lx) lz = lx;
180 80221 : return lz;
181 : }
182 13656864 : if (e < 0) { swap(x, y); e = -e; }
183 13656864 : lx = realprec(x);
184 13656864 : ly = realprec(y);
185 13656864 : if (e) {
186 10956995 : long d = nbits2extraprec(e), l = ly-d;
187 10957001 : return (l > lx)? lx + d: ly;
188 : }
189 2699869 : return minss(lx, ly);
190 : }
191 : long
192 677124998 : precision(GEN z)
193 : {
194 677124998 : switch(typ(z))
195 : {
196 665545921 : case t_REAL: return precREAL(z);
197 11549040 : case t_COMPLEX: return precCOMPLEX(z);
198 : }
199 30037 : return 0;
200 : }
201 :
202 : long
203 9718643 : gprecision(GEN x)
204 : {
205 : long i, k, l;
206 :
207 9718643 : switch(typ(x))
208 : {
209 2196159 : case t_REAL: return precREAL(x);
210 4081606 : case t_COMPLEX: return precCOMPLEX(x);
211 872316 : case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT:
212 : case t_PADIC: case t_QUAD: case t_POLMOD:
213 872316 : return 0;
214 :
215 182 : case t_POL: case t_SER:
216 182 : k = LONG_MAX;
217 826 : for (i=lg(x)-1; i>1; i--)
218 : {
219 644 : l = gprecision(gel(x,i));
220 644 : if (l && l<k) k = l;
221 : }
222 182 : return (k==LONG_MAX)? 0: k;
223 2568441 : case t_VEC: case t_COL: case t_MAT:
224 2568441 : k = LONG_MAX;
225 10291784 : for (i=lg(x)-1; i>0; i--)
226 : {
227 7723338 : l = gprecision(gel(x,i));
228 7723343 : if (l && l<k) k = l;
229 : }
230 2568446 : 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 + valp(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 101930 : poldegree(GEN x, long v)
370 : {
371 101930 : const long DEGREE0 = -LONG_MAX;
372 101930 : long tx = typ(x), lx,w,i,d;
373 :
374 101930 : if (is_scalar_t(tx)) return gequal0(x)? DEGREE0: 0;
375 101576 : switch(tx)
376 : {
377 101499 : case t_POL:
378 101499 : if (!signe(x)) return DEGREE0;
379 101492 : w = varn(x);
380 101492 : 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 77 : case t_RFRAC:
391 : {
392 77 : GEN a = gel(x,1), b = gel(x,2);
393 77 : if (gequal0(a)) return DEGREE0;
394 70 : if (v < 0)
395 : {
396 70 : v = varn(b); d = -degpol(b);
397 70 : if (typ(a) == t_POL && varn(a) == v) d += degpol(a);
398 70 : 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 30811 : gppoldegree(GEN x, long v)
408 : {
409 30811 : long d = poldegree(x,v);
410 30811 : 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 484696 : RgX_degree(GEN x, long v)
416 : {
417 484696 : long tx = typ(x), lx, w, i, d;
418 :
419 484696 : if (is_scalar_t(tx)) return gequal0(x)? -1: 0;
420 285900 : switch(tx)
421 : {
422 285900 : case t_POL:
423 285900 : if (!signe(x)) return -1;
424 285879 : w = varn(x);
425 285879 : if (v == w) return degpol(x);
426 106609 : if (varncmp(v, w) < 0) return 0;
427 106609 : lx = lg(x); d = -1;
428 481189 : for (i=2; i<lx; i++)
429 : {
430 374580 : long e = RgX_degree(gel(x,i), v);
431 374580 : if (e > d) d = e;
432 : }
433 106609 : 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 259 : pollead(GEN x, long v)
454 : {
455 259 : long tx = typ(x), w;
456 : pari_sp av;
457 :
458 259 : if (is_scalar_t(tx)) return gcopy(x);
459 259 : w = varn(x);
460 259 : switch(tx)
461 : {
462 224 : case t_POL:
463 224 : if (v < 0 || v == w)
464 : {
465 189 : long l = lg(x);
466 189 : 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, valp(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 14049 : isinexactreal(GEN x)
489 : {
490 : long i;
491 14049 : switch(typ(x))
492 : {
493 973 : 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 9891 : case t_INT: case t_INTMOD: case t_FRAC:
497 : case t_FFELT: case t_PADIC: case t_QUAD:
498 9891 : 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 1126876 : isrealappr(GEN x, long e)
518 : {
519 : long i;
520 1126876 : switch(typ(x))
521 : {
522 251815 : case t_INT: case t_REAL: case t_FRAC:
523 251815 : return 1;
524 875065 : case t_COMPLEX:
525 875065 : 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 133243813 : isinexact(GEN x)
547 : {
548 : long lx, i;
549 :
550 133243813 : switch(typ(x))
551 : {
552 653049 : case t_REAL: case t_PADIC: case t_SER:
553 653049 : return 1;
554 90029646 : case t_INT: case t_INTMOD: case t_FFELT: case t_FRAC:
555 : case t_QFB:
556 90029646 : return 0;
557 2410045 : case t_COMPLEX: case t_QUAD: case t_RFRAC: case t_POLMOD:
558 2410045 : return isinexact(gel(x,1)) || isinexact(gel(x,2));
559 40151073 : case t_POL:
560 125938368 : for (i=lg(x)-1; i>1; i--)
561 86065867 : if (isinexact(gel(x,i))) return 1;
562 39872501 : return 0;
563 0 : case t_VEC: case t_COL: case t_MAT:
564 0 : for (i=lg(x)-1; i>0; i--)
565 0 : if (isinexact(gel(x,i))) return 1;
566 0 : 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 176636 : is_realext(GEN x)
613 176636 : { long t = typ(x);
614 176636 : 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)) { set_avma(av); return stoi(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 3436373 : modRr_safe(GEN x, GEN y)
732 : {
733 : GEN q, f;
734 : long e;
735 3436373 : if (isintzero(x)) return gen_0;
736 3436384 : q = gdiv(x,y); /* t_REAL */
737 :
738 3436252 : e = expo(q);
739 3436252 : if (e >= 0 && nbits2prec(e+1) > realprec(q)) return NULL;
740 3436248 : f = floorr(q);
741 3436151 : if (signe(y) < 0 && signe(subri(q,f))) f = addiu(f, 1);
742 3436319 : return signe(f)? gsub(x, mulir(f,y)): x;
743 : }
744 :
745 : GEN
746 35879670 : gmod(GEN x, GEN y)
747 : {
748 : pari_sp av;
749 : long ty, tx;
750 : GEN z;
751 :
752 35879670 : tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gmodsg(itos(x),y);
753 1145417 : ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gmodgs(x,itos(y));
754 1745195 : if (is_matvec_t(tx)) pari_APPLY_same(gmod(gel(x,i), y));
755 807527 : if (tx == t_POL || ty == t_POL) return grem(x,y);
756 511505 : if (!is_scalar_t(tx) || !is_scalar_t(ty)) pari_err_TYPE2("%",x,y);
757 511442 : switch(ty)
758 : {
759 510938 : case t_INT:
760 : switch(tx)
761 : {
762 507634 : case t_INT: return modii(x,y);
763 7 : case t_INTMOD: z=cgetg(3, t_INTMOD);
764 7 : gel(z,1) = gcdii(gel(x,1),y);
765 7 : gel(z,2) = modii(gel(x,2),gel(z,1)); return z;
766 491 : case t_FRAC: return Fp_div(gel(x,1),gel(x,2),y);
767 2771 : case t_PADIC: return padic_to_Fp(x, y);
768 14 : case t_QUAD: if (!is_realquad(x)) break;
769 : case t_REAL:
770 14 : av = avma; /* NB: conflicting semantic with lift(x * Mod(1,y)). */
771 14 : return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
772 : }
773 21 : break;
774 126 : case t_QUAD:
775 126 : if (!is_realquad(y)) break;
776 : case t_REAL: case t_FRAC:
777 189 : if (!is_realext(x)) break;
778 84 : av = avma;
779 84 : return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
780 : }
781 441 : pari_err_TYPE2("%",x,y);
782 : return NULL; /* LCOV_EXCL_LINE */
783 : }
784 :
785 : GEN
786 22027554 : gmodgs(GEN x, long y)
787 : {
788 : ulong u;
789 22027554 : long i, tx = typ(x);
790 : GEN z;
791 43821506 : if (is_matvec_t(tx)) pari_APPLY_same(gmodgs(gel(x,i), y));
792 19649334 : if (!y) pari_err_INV("gmodgs",gen_0);
793 19649334 : switch(tx)
794 : {
795 19566419 : case t_INT: return modis(x,y);
796 14 : case t_REAL: return modrs(x,y);
797 :
798 21 : case t_INTMOD: z=cgetg(3, t_INTMOD);
799 21 : u = (ulong)labs(y);
800 21 : i = ugcdiu(gel(x,1), u);
801 21 : gel(z,1) = utoi(i);
802 21 : gel(z,2) = modis(gel(x,2), i); return z;
803 :
804 81360 : case t_FRAC:
805 81360 : u = (ulong)labs(y);
806 81360 : return utoi( Fl_div(umodiu(gel(x,1), u),
807 81360 : umodiu(gel(x,2), u), u) );
808 28 : case t_QUAD:
809 : {
810 28 : pari_sp av = avma;
811 28 : if (!is_realquad(x)) break;
812 14 : return gerepileupto(av, gsub(x, gmulgs(_quotqs(x,y),y)));
813 : }
814 1436 : case t_PADIC: return padic_to_Fp(x, stoi(y));
815 14 : case t_POL: return scalarpol(Rg_get_0(x), varn(x));
816 14 : case t_POLMOD: return gmul(gen_0,x);
817 : }
818 42 : pari_err_TYPE2("%",x,stoi(y));
819 : return NULL; /* LCOV_EXCL_LINE */
820 : }
821 : GEN
822 34734253 : gmodsg(long x, GEN y)
823 : {
824 34734253 : switch(typ(y))
825 : {
826 34733889 : case t_INT: return modsi(x,y);
827 35 : case t_REAL: return modsr(x,y);
828 35 : case t_FRAC: return modsf(x,y);
829 63 : case t_QUAD:
830 : {
831 63 : pari_sp av = avma;
832 63 : if (!is_realquad(y)) break;
833 35 : return gerepileupto(av, gsubsg(x, gmul(_quotsq(x,y),y)));
834 : }
835 119 : case t_POL:
836 119 : if (!signe(y)) pari_err_INV("gmodsg",y);
837 119 : return degpol(y)? gmulsg(x, Rg_get_1(y)): Rg_get_0(y);
838 : }
839 140 : pari_err_TYPE2("%",stoi(x),y);
840 : return NULL; /* LCOV_EXCL_LINE */
841 : }
842 : /* divisibility: return 1 if y | x, 0 otherwise */
843 : int
844 15539 : gdvd(GEN x, GEN y)
845 : {
846 15539 : pari_sp av = avma;
847 15539 : return gc_bool(av, gequal0( gmod(x,y) ));
848 : }
849 :
850 : GEN
851 689981 : gmodulss(long x, long y)
852 : {
853 689981 : if (!y) pari_err_INV("%",gen_0);
854 689974 : y = labs(y);
855 689974 : retmkintmod(utoi(umodsu(x, y)), utoipos(y));
856 : }
857 : GEN
858 1000005 : gmodulsg(long x, GEN y)
859 : {
860 1000005 : switch(typ(y))
861 : {
862 752379 : case t_INT:
863 752379 : if (!is_bigint(y)) return gmodulss(x,itos(y));
864 62411 : retmkintmod(modsi(x,y), absi(y));
865 247619 : case t_POL:
866 247619 : if (!signe(y)) pari_err_INV("%", y);
867 247612 : retmkpolmod(degpol(y)? stoi(x): gen_0,RgX_copy(y));
868 : }
869 7 : pari_err_TYPE2("%",stoi(x),y);
870 : return NULL; /* LCOV_EXCL_LINE */
871 : }
872 : GEN
873 1516489 : gmodulo(GEN x,GEN y)
874 : {
875 1516489 : long tx = typ(x), vx, vy;
876 1516489 : if (tx == t_INT && !is_bigint(x)) return gmodulsg(itos(x), y);
877 1062819 : if (is_matvec_t(tx)) pari_APPLY_same(gmodulo(gel(x,i), y));
878 511499 : switch(typ(y))
879 : {
880 98191 : case t_INT:
881 98191 : if (!is_const_t(tx)) return gmul(x, gmodulsg(1,y));
882 98149 : if (tx == t_INTMOD) return gmod(x,y);
883 98142 : retmkintmod(Rg_to_Fp(x,y), absi(y));
884 413308 : case t_POL:
885 413308 : vx = gvar(x); vy = varn(y);
886 413308 : if (varncmp(vy, vx) > 0) return gmul(x, gmodulsg(1,y));
887 409556 : if (vx == vy && tx == t_POLMOD) return grem(x,y);
888 396242 : retmkpolmod(grem(x,y), RgX_copy(y));
889 : }
890 0 : pari_err_TYPE2("%",x,y);
891 : return NULL; /* LCOV_EXCL_LINE */
892 : }
893 :
894 : /*******************************************************************/
895 : /* */
896 : /* GENERIC EUCLIDEAN DIVISION */
897 : /* */
898 : /*******************************************************************/
899 : GEN
900 6199466 : gdivent(GEN x, GEN y)
901 : {
902 : long tx, ty;
903 6199466 : tx = typ(x); if (tx == t_INT && !is_bigint(x)) return gdiventsg(itos(x),y);
904 2122 : ty = typ(y); if (ty == t_INT && !is_bigint(y)) return gdiventgs(x,itos(y));
905 1960 : if (is_matvec_t(tx)) pari_APPLY_same(gdivent(gel(x,i), y));
906 1610 : if (tx == t_POL || ty == t_POL) return gdeuc(x,y);
907 1148 : switch(ty)
908 : {
909 112 : case t_INT:
910 : switch(tx)
911 : {
912 7 : case t_INT: return truedivii(x,y);
913 7 : case t_REAL: return quotri(x,y);
914 7 : case t_FRAC: return quotfi(x,y);
915 21 : case t_QUAD:
916 21 : if (!is_realquad(x)) break;
917 7 : return quot(x,y);
918 : }
919 84 : break;
920 252 : case t_QUAD:
921 252 : if (!is_realext(x) || !is_realquad(y)) break;
922 : case t_REAL: case t_FRAC:
923 252 : return quot(x,y);
924 : }
925 868 : pari_err_TYPE2("\\",x,y);
926 : return NULL; /* LCOV_EXCL_LINE */
927 : }
928 :
929 : GEN
930 284053 : gdiventgs(GEN x, long y)
931 : {
932 284053 : switch(typ(x))
933 : {
934 235663 : case t_INT: return truedivis(x,y);
935 14 : case t_REAL: return quotrs(x,y);
936 301 : case t_FRAC: return quotfs(x,y);
937 42 : case t_QUAD: if (!is_realquad(x)) break;
938 14 : return quotqs(x,y);
939 28 : case t_POL: return gdivgs(x,y);
940 275043 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gdiventgs(gel(x,i),y));
941 : }
942 168 : pari_err_TYPE2("\\",x,stoi(y));
943 : return NULL; /* LCOV_EXCL_LINE */
944 : }
945 : GEN
946 6197344 : gdiventsg(long x, GEN y)
947 : {
948 6197344 : switch(typ(y))
949 : {
950 6196889 : case t_INT: return truedivsi(x,y);
951 35 : case t_REAL: return quotsr(x,y);
952 35 : case t_FRAC: return quotsf(x,y);
953 91 : case t_QUAD: if (!is_realquad(y)) break;
954 35 : return quotsq(x,y);
955 70 : case t_POL:
956 70 : if (!signe(y)) pari_err_INV("gdiventsg",y);
957 70 : return degpol(y)? Rg_get_0(y): gdivsg(x,gel(y,2));
958 : }
959 280 : pari_err_TYPE2("\\",stoi(x),y);
960 : return NULL; /* LCOV_EXCL_LINE */
961 : }
962 :
963 : /* with remainder */
964 : static GEN
965 518 : quotrem(GEN x, GEN y, GEN *r)
966 : {
967 518 : GEN q = quot(x,y);
968 448 : pari_sp av = avma;
969 448 : *r = gerepileupto(av, gsub(x, gmul(q,y)));
970 448 : return q;
971 : }
972 :
973 : GEN
974 1064 : gdiventres(GEN x, GEN y)
975 : {
976 1064 : long tx = typ(x), ty = typ(y);
977 : GEN z;
978 :
979 1078 : if (is_matvec_t(tx)) pari_APPLY_same(gdiventres(gel(x,i), y));
980 1057 : z = cgetg(3,t_COL);
981 1057 : if (tx == t_POL || ty == t_POL)
982 : {
983 182 : gel(z,1) = poldivrem(x,y,(GEN*)(z+2));
984 168 : return z;
985 : }
986 875 : switch(ty)
987 : {
988 252 : case t_INT:
989 : switch(tx)
990 : { /* equal to, but more efficient than next case */
991 84 : case t_INT:
992 84 : gel(z,1) = truedvmdii(x,y,(GEN*)(z+2));
993 84 : return z;
994 42 : case t_QUAD:
995 42 : if (!is_realquad(x)) break;
996 : case t_REAL: case t_FRAC:
997 63 : gel(z,1) = quotrem(x,y,&gel(z,2));
998 63 : return z;
999 : }
1000 105 : break;
1001 154 : case t_QUAD:
1002 154 : if (!is_realext(x) || !is_realquad(y)) break;
1003 : case t_REAL: case t_FRAC:
1004 196 : gel(z,1) = quotrem(x,y,&gel(z,2));
1005 126 : return z;
1006 : }
1007 532 : pari_err_TYPE2("\\",x,y);
1008 : return NULL; /* LCOV_EXCL_LINE */
1009 : }
1010 :
1011 : GEN
1012 1057 : divrem(GEN x, GEN y, long v)
1013 : {
1014 1057 : pari_sp av = avma;
1015 : long vx, vy;
1016 : GEN q, r;
1017 1057 : if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
1018 7 : vx = varn(x); if (vx != v) x = swap_vars(x,v);
1019 7 : vy = varn(y); if (vy != v) y = swap_vars(y,v);
1020 7 : q = RgX_divrem(x,y, &r);
1021 7 : if (v && (vx != v || vy != v))
1022 : {
1023 7 : GEN X = pol_x(v);
1024 7 : q = gsubst(q, v, X); /* poleval broken for t_RFRAC, subst is safe */
1025 7 : r = gsubst(r, v, X);
1026 : }
1027 7 : return gerepilecopy(av, mkcol2(q, r));
1028 : }
1029 :
1030 : GEN
1031 58270151 : diviiround(GEN x, GEN y)
1032 : {
1033 58270151 : pari_sp av1, av = avma;
1034 : GEN q,r;
1035 : int fl;
1036 :
1037 58270151 : q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
1038 58243883 : if (r==gen_0) return q;
1039 31749344 : av1 = avma;
1040 31749344 : fl = abscmpii(shifti(r,1),y);
1041 31773722 : set_avma(av1); cgiv(r);
1042 31776824 : if (fl >= 0) /* If 2*|r| >= |y| */
1043 : {
1044 16275894 : long sz = signe(x)*signe(y);
1045 16275894 : if (fl || sz > 0) q = gerepileuptoint(av, addis(q,sz));
1046 : }
1047 31781400 : return q;
1048 : }
1049 :
1050 : static GEN
1051 518 : _abs(GEN x)
1052 : {
1053 518 : if (typ(x) == t_QUAD) return (gsigne(x) < 0)? gneg(x): x;
1054 364 : return R_abs_shallow(x);
1055 : }
1056 :
1057 : /* If x and y are not both scalars, same as gdivent.
1058 : * Otherwise, compute the quotient x/y, rounded to the nearest integer
1059 : * (towards +oo in case of tie). */
1060 : GEN
1061 1462594 : gdivround(GEN x, GEN y)
1062 : {
1063 : pari_sp av;
1064 1462594 : long tx = typ(x), ty = typ(y);
1065 : GEN q, r;
1066 :
1067 1462594 : if (tx == t_INT && ty == t_INT) return diviiround(x,y);
1068 175433 : av = avma;
1069 175433 : if (is_realext(x) && is_realext(y))
1070 : { /* same as diviiround, less efficient */
1071 : pari_sp av1;
1072 : int fl;
1073 259 : q = quotrem(x,y,&r); av1 = avma;
1074 259 : fl = gcmp(gmul2n(_abs(r),1), _abs(y));
1075 259 : set_avma(av1); cgiv(r);
1076 259 : if (fl >= 0) /* If 2*|r| >= |y| */
1077 : {
1078 84 : long sz = gsigne(y);
1079 84 : if (fl || sz > 0) q = gerepileupto(av, gaddgs(q, sz));
1080 : }
1081 259 : return q;
1082 : }
1083 1584771 : if (is_matvec_t(tx)) pari_APPLY_same(gdivround(gel(x,i),y));
1084 931 : return gdivent(x,y);
1085 : }
1086 :
1087 : GEN
1088 0 : gdivmod(GEN x, GEN y, GEN *pr)
1089 : {
1090 0 : switch(typ(x))
1091 : {
1092 0 : case t_INT:
1093 0 : switch(typ(y))
1094 : {
1095 0 : case t_INT: return dvmdii(x,y,pr);
1096 0 : case t_POL: *pr=icopy(x); return gen_0;
1097 : }
1098 0 : break;
1099 0 : case t_POL: return poldivrem(x,y,pr);
1100 : }
1101 0 : pari_err_TYPE2("gdivmod",x,y);
1102 : return NULL;/*LCOV_EXCL_LINE*/
1103 : }
1104 :
1105 : /*******************************************************************/
1106 : /* */
1107 : /* SHIFT */
1108 : /* */
1109 : /*******************************************************************/
1110 :
1111 : /* Shift tronque si n<0 (multiplication tronquee par 2^n) */
1112 :
1113 : GEN
1114 43132122 : gshift(GEN x, long n)
1115 : {
1116 43132122 : switch(typ(x))
1117 : {
1118 36429822 : case t_INT: return shifti(x,n);
1119 6552351 : case t_REAL:return shiftr(x,n);
1120 35035 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gshift(gel(x,i),n));
1121 : }
1122 139526 : return gmul2n(x,n);
1123 : }
1124 :
1125 : /*******************************************************************/
1126 : /* */
1127 : /* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
1128 : /* */
1129 : /*******************************************************************/
1130 :
1131 : /* Convert t_SER --> t_POL, ignoring valp. INTERNAL ! */
1132 : GEN
1133 9619978 : ser2pol_i(GEN x, long lx)
1134 : {
1135 9619978 : long i = lx-1;
1136 : GEN y;
1137 13162020 : while (i > 1 && isrationalzero(gel(x,i))) i--;
1138 9619978 : if (!signe(x))
1139 : { /* danger */
1140 105 : if (i == 1) return zeropol(varn(x));
1141 105 : y = cgetg(3,t_POL); y[1] = x[1] & ~VALPBITS;
1142 105 : gel(y,2) = gel(x,2); return y;
1143 : }
1144 9619873 : y = cgetg(i+1, t_POL); y[1] = x[1] & ~VALPBITS;
1145 40603278 : for ( ; i > 1; i--) gel(y,i) = gel(x,i);
1146 9619873 : return y;
1147 : }
1148 :
1149 : GEN
1150 731589 : ser2pol_approx(GEN x, long l, long *v)
1151 : {
1152 731589 : long i = 2, j = l-1, k;
1153 : GEN y;
1154 731624 : while (i < l && gequal0(gel(x,i))) i++;
1155 731589 : if (i != 2) pari_warn(warner,"normalizing a series with 0 leading term");
1156 731589 : *v = i - 2; if (i == l) return zeropol(varn(x));
1157 971352 : while (j > i && gequal0(gel(x,j))) j--;
1158 731575 : l = j - *v + 1;
1159 731575 : y = cgetg(l, t_POL); y[1] = x[1] & ~VALPBITS;
1160 3757394 : k = l; while (k > 2) gel(y, --k) = gel(x,j--);
1161 731575 : return y;
1162 : }
1163 :
1164 : GEN
1165 37849 : ser_inv(GEN b)
1166 : {
1167 37849 : pari_sp av = avma;
1168 37849 : long e, l = lg(b);
1169 : GEN x, y;
1170 37849 : y = ser2pol_approx(b, l, &e);
1171 37849 : y = RgXn_inv_i(y, l - 2 - e);
1172 37842 : x = RgX_to_ser(y, l); setvalp(x, - valp(b) - e);
1173 37842 : return gerepilecopy(av, x);
1174 : }
1175 :
1176 : /* T t_POL in var v, mod out by T components of x which are t_POL/t_RFRAC in v.
1177 : * Recursively. Make sure that resulting polynomials of degree 0 in v are
1178 : * simplified (map K[X]_0 to K) */
1179 : static GEN
1180 196 : mod_r(GEN x, long v, GEN T)
1181 : {
1182 196 : long w, tx = typ(x);
1183 : GEN y;
1184 :
1185 196 : if (is_const_t(tx)) return x;
1186 175 : switch(tx)
1187 : {
1188 7 : case t_POLMOD:
1189 7 : w = varn(gel(x,1));
1190 7 : if (w == v) pari_err_PRIORITY("subst", gel(x,1), "=", v);
1191 7 : if (varncmp(v, w) < 0) return x;
1192 7 : return gmodulo(mod_r(gel(x,2),v,T), mod_r(gel(x,1),v,T));
1193 7 : case t_SER:
1194 7 : w = varn(x);
1195 7 : if (w == v) break; /* fail */
1196 7 : if (varncmp(v, w) < 0 || ser_isexactzero(x)) return x;
1197 21 : pari_APPLY_ser(mod_r(gel(x,i),v,T));
1198 133 : case t_POL:
1199 133 : w = varn(x);
1200 133 : if (w == v)
1201 : {
1202 105 : x = RgX_rem(x, T);
1203 105 : if (!degpol(x)) x = gel(x,2);
1204 105 : return x;
1205 : }
1206 28 : if (varncmp(v, w) < 0) return x;
1207 98 : pari_APPLY_pol(mod_r(gel(x,i),v,T));
1208 14 : case t_RFRAC:
1209 14 : x = gdiv(mod_r(gel(x,1),v,T), mod_r(gel(x,2),v,T));
1210 14 : if (typ(x) == t_POL && varn(x) == v && lg(x) == 3) x = gel(x,2);
1211 14 : return x;
1212 7 : case t_VEC: case t_COL: case t_MAT:
1213 21 : pari_APPLY_same(mod_r(gel(x,i),v,T));
1214 7 : case t_LIST:
1215 7 : y = mklist();
1216 7 : list_data(y) = list_data(x)? mod_r(list_data(x),v,T): NULL;
1217 7 : return y;
1218 : }
1219 0 : pari_err_TYPE("substpol",x);
1220 : return NULL;/*LCOV_EXCL_LINE*/
1221 : }
1222 : GEN
1223 2037 : gsubstpol(GEN x, GEN T, GEN y)
1224 : {
1225 2037 : pari_sp av = avma;
1226 : long v;
1227 : GEN z;
1228 2037 : if (typ(T) == t_POL && RgX_is_monomial(T) && gequal1(leading_coeff(T)))
1229 : { /* T = t^d */
1230 2016 : long d = degpol(T);
1231 2016 : v = varn(T); z = (d==1)? x: gdeflate(x, v, d);
1232 2002 : if (z) return gerepileupto(av, gsubst(z, v, y));
1233 : }
1234 49 : v = fetch_var(); T = simplify_shallow(T);
1235 49 : if (typ(T) == t_RFRAC)
1236 21 : z = gsub(gel(T,1), gmul(pol_x(v), gel(T,2)));
1237 : else
1238 28 : z = gsub(T, pol_x(v));
1239 49 : z = mod_r(x, gvar(T), z);
1240 49 : z = gsubst(z, v, y); (void)delete_var();
1241 49 : return gerepileupto(av, z);
1242 : }
1243 :
1244 : long
1245 295205 : RgX_deflate_order(GEN x)
1246 : {
1247 295205 : ulong d = 0, i, lx = (ulong)lg(x);
1248 1340621 : for (i=3; i<lx; i++)
1249 1228287 : if (!gequal0(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1250 112334 : return d? (long)d: 1;
1251 : }
1252 : long
1253 508481 : ZX_deflate_order(GEN x)
1254 : {
1255 508481 : ulong d = 0, i, lx = (ulong)lg(x);
1256 1610814 : for (i=3; i<lx; i++)
1257 1423591 : if (signe(gel(x,i))) { d = ugcd(d,i-2); if (d == 1) return 1; }
1258 187223 : return d? (long)d: 1;
1259 : }
1260 :
1261 : /* deflate (non-leaf) x recursively */
1262 : static GEN
1263 63 : vdeflate(GEN x, long v, long d)
1264 : {
1265 63 : long i = lontyp[typ(x)], lx;
1266 63 : GEN z = cgetg_copy(x, &lx);
1267 63 : if (i == 2) z[1] = x[1];
1268 154 : for (; i<lx; i++)
1269 : {
1270 133 : gel(z,i) = gdeflate(gel(x,i),v,d);
1271 133 : if (!z[i]) return NULL;
1272 : }
1273 21 : return z;
1274 : }
1275 :
1276 : /* don't return NULL if substitution fails (fallback won't be able to handle
1277 : * t_SER anyway), fail with a meaningful message */
1278 : static GEN
1279 2793 : serdeflate(GEN x, long v, long d)
1280 : {
1281 2793 : long V, dy, lx, vx = varn(x);
1282 : pari_sp av;
1283 : GEN y;
1284 2793 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1285 2786 : if (varncmp(vx, v) > 0) return gcopy(x);
1286 2786 : av = avma;
1287 2786 : V = valp(x);
1288 2786 : lx = lg(x);
1289 2786 : if (lx == 2) return zeroser(v, V / d);
1290 2786 : y = ser2pol_i(x, lx);
1291 2786 : dy = degpol(y);
1292 2786 : if (V % d != 0 || (dy > 0 && RgX_deflate_order(y) % d != 0))
1293 : {
1294 14 : const char *s = stack_sprintf("valuation(x) %% %ld", d);
1295 14 : pari_err_DOMAIN("gdeflate", s, "!=", gen_0,x);
1296 : }
1297 2772 : if (dy > 0) y = RgX_deflate(y, d);
1298 2772 : y = RgX_to_ser(y, 3 + (lx-3)/d);
1299 2772 : setvalp(y, V/d); return gerepilecopy(av, y);
1300 : }
1301 : static GEN
1302 2051 : poldeflate(GEN x, long v, long d)
1303 : {
1304 2051 : long vx = varn(x);
1305 : pari_sp av;
1306 2051 : if (varncmp(vx, v) < 0) return vdeflate(x,v,d);
1307 2023 : if (varncmp(vx, v) > 0 || degpol(x) <= 0) return gcopy(x);
1308 1988 : av = avma;
1309 : /* x nonconstant */
1310 1988 : if (RgX_deflate_order(x) % d != 0) return NULL;
1311 1960 : return gerepilecopy(av, RgX_deflate(x,d));
1312 : }
1313 : static GEN
1314 21 : listdeflate(GEN x, long v, long d)
1315 : {
1316 21 : GEN y = NULL, z = mklist();
1317 21 : if (list_data(x))
1318 : {
1319 14 : y = vdeflate(list_data(x),v,d);
1320 14 : if (!y) return NULL;
1321 : }
1322 14 : list_data(z) = y; return z;
1323 : }
1324 : /* return NULL if substitution fails */
1325 : GEN
1326 4907 : gdeflate(GEN x, long v, long d)
1327 : {
1328 4907 : if (d <= 0) pari_err_DOMAIN("gdeflate", "degree", "<=", gen_0,stoi(d));
1329 4907 : switch(typ(x))
1330 : {
1331 28 : case t_INT:
1332 : case t_REAL:
1333 : case t_INTMOD:
1334 : case t_FRAC:
1335 : case t_FFELT:
1336 : case t_COMPLEX:
1337 : case t_PADIC:
1338 28 : case t_QUAD: return gcopy(x);
1339 2051 : case t_POL: return poldeflate(x,v,d);
1340 2793 : case t_SER: return serdeflate(x,v,d);
1341 7 : case t_POLMOD:
1342 7 : if (varncmp(varn(gel(x,1)), v) >= 0) return gcopy(x);
1343 : /* fall through */
1344 : case t_RFRAC:
1345 : case t_VEC:
1346 : case t_COL:
1347 14 : case t_MAT: return vdeflate(x,v,d);
1348 21 : case t_LIST: return listdeflate(x,v,d);
1349 : }
1350 0 : pari_err_TYPE("gdeflate",x);
1351 : return NULL; /* LCOV_EXCL_LINE */
1352 : }
1353 :
1354 : /* set *m to the largest d such that x0 = A(X^d); return A */
1355 : GEN
1356 291173 : RgX_deflate_max(GEN x, long *m)
1357 : {
1358 291173 : *m = RgX_deflate_order(x);
1359 291171 : return RgX_deflate(x, *m);
1360 : }
1361 : GEN
1362 305161 : ZX_deflate_max(GEN x, long *m)
1363 : {
1364 305161 : *m = ZX_deflate_order(x);
1365 305161 : return RgX_deflate(x, *m);
1366 : }
1367 :
1368 : static int
1369 20181 : serequalXk(GEN x)
1370 : {
1371 20181 : long i, l = lg(x);
1372 20181 : if (l == 2 || !isint1(gel(x,2))) return 0;
1373 9870 : for (i = 3; i < l; i++)
1374 7917 : if (!isintzero(gel(x,i))) return 0;
1375 1953 : return 1;
1376 : }
1377 :
1378 : static GEN
1379 84 : gsubst_v(GEN e, long v, GEN x)
1380 245 : { pari_APPLY_same(gsubst(e, v, gel(x,i))); }
1381 :
1382 : static GEN
1383 14 : constmat(GEN z, long n)
1384 : {
1385 14 : GEN y = cgetg(n+1, t_MAT), c = const_col(n, gcopy(z));
1386 : long i;
1387 35 : for (i = 1; i <= n; i++) gel(y, i) = c;
1388 14 : return y;
1389 : }
1390 : static GEN
1391 56 : scalarmat2(GEN o, GEN z, long n)
1392 : {
1393 : GEN y;
1394 : long i;
1395 56 : if (n == 0) return cgetg(1, t_MAT);
1396 56 : if (n == 1) retmkmat(mkcol(gcopy(o)));
1397 35 : y = cgetg(n+1, t_MAT); z = gcopy(z); o = gcopy(o);
1398 105 : for (i = 1; i <= n; i++) { gel(y, i) = const_col(n, z); gcoeff(y,i,i) = o; }
1399 35 : return y;
1400 : }
1401 : /* x * y^0, n = dim(y) if t_MAT, else -1 */
1402 : static GEN
1403 745746 : subst_higher(GEN x, GEN y, long n)
1404 : {
1405 745746 : GEN o = Rg_get_1(y);
1406 745746 : if (o == gen_1) return n < 0? gcopy(x): scalarmat(x,n);
1407 63 : x = gmul(x,o); return n < 0? x: scalarmat2(x, Rg_get_0(y), n);
1408 : }
1409 :
1410 : /* x t_POLMOD, v strictly lower priority than var(x) */
1411 : static GEN
1412 10948 : subst_polmod(GEN x, long v, GEN y)
1413 : {
1414 : long l, i;
1415 10948 : GEN a = gsubst(gel(x,2),v,y), b = gsubst(gel(x,1),v,y), z;
1416 :
1417 10948 : if (typ(b) != t_POL) pari_err_TYPE2("substitution",x,y);
1418 10948 : if (typ(a) != t_POL || varncmp(varn(a), varn(b)) >= 0) return gmodulo(a, b);
1419 511 : l = lg(a); z = cgetg(l,t_POL); z[1] = a[1];
1420 3948 : for (i = 2; i < l; i++) gel(z,i) = gmodulo(gel(a,i),b);
1421 511 : return normalizepol_lg(z, l);
1422 : }
1423 : /* Trunc to n terms; x + O(t^(n + v(x))). FIXME: export ? */
1424 : static GEN
1425 70 : sertrunc(GEN x, long n)
1426 : {
1427 70 : long i, l = n + 2;
1428 : GEN y;
1429 70 : if (l >= lg(x)) return x;
1430 14 : if (n <= 0) return zeroser(varn(x), n + valp(x));
1431 14 : y = cgetg(l, t_SER);
1432 28 : for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
1433 14 : y[1] = x[1]; return y;
1434 : }
1435 : /* FIXME: export ? */
1436 : static GEN
1437 1960 : sertrunc_copy(GEN x, long n)
1438 : {
1439 1960 : long i, l = minss(n + 2, lg(x));
1440 1960 : GEN y = cgetg(l, t_SER);
1441 13433 : for (i = 2; i < l; i++) gel(y,i) = gcopy(gel(x,i));
1442 1960 : y[1] = x[1]; return y;
1443 : }
1444 :
1445 : GEN
1446 2721468 : gsubst(GEN x, long v, GEN y)
1447 : {
1448 2721468 : long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
1449 : long l, vx, vy, ex, ey, i, j, k, jb, matn;
1450 : pari_sp av, av2;
1451 : GEN X, t, z;
1452 :
1453 2721468 : switch(ty)
1454 : {
1455 84 : case t_VEC: case t_COL:
1456 84 : return gsubst_v(x, v, y);
1457 175 : case t_MAT:
1458 175 : if (ly==1) return cgetg(1,t_MAT);
1459 168 : if (ly == lgcols(y)) { matn = ly - 1; break; }
1460 : /* fall through */
1461 : case t_QFB:
1462 7 : pari_err_TYPE2("substitution",x,y);
1463 2721209 : default: matn = -1;
1464 : }
1465 2721370 : if (is_scalar_t(tx))
1466 : {
1467 408493 : if (tx == t_POLMOD && varncmp(v, varn(gel(x,1))) > 0)
1468 : {
1469 10948 : av = avma;
1470 10948 : return gerepileupto(av, subst_polmod(x, v, y));
1471 : }
1472 397545 : return subst_higher(x, y, matn);
1473 : }
1474 :
1475 2312877 : switch(tx)
1476 : {
1477 2074760 : case t_POL:
1478 2074760 : vx = varn(x);
1479 2074760 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1480 2072667 : if (varncmp(vx, v) < 0)
1481 : {
1482 175167 : av = avma; z = cgetg(lx, t_POL); z[1] = x[1];
1483 175167 : if (lx == 2) return z;
1484 814380 : for (i = 2; i < lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
1485 174726 : z = normalizepol_lg(z, lx); lx = lg(z);
1486 174726 : if (lx == 2) { set_avma(av); return zeropol(vx); }
1487 174712 : if (lx == 3) return gerepileupto(av, gmul(pol_1(vx), gel(z,2)));
1488 152410 : return gerepileupto(av, poleval(z, pol_x(vx)));
1489 : }
1490 : /* v = vx */
1491 1897500 : if (lx == 2)
1492 : {
1493 27860 : GEN z = Rg_get_0(y);
1494 27860 : return matn >= 0? constmat(z, matn): z;
1495 : }
1496 1869640 : if (lx == 3)
1497 : {
1498 346101 : x = subst_higher(gel(x,2), y, matn);
1499 346101 : if (matn >= 0) return x;
1500 346087 : vy = gvar(y);
1501 346087 : return (vy == NO_VARIABLE)? x: gmul(x, pol_1(vy));
1502 : }
1503 1523539 : return matn >= 0? RgX_RgM_eval(x, y): poleval(x,y);
1504 :
1505 24549 : case t_SER:
1506 24549 : vx = varn(x);
1507 24549 : if (varncmp(vx, v) > 0) return subst_higher(x, y, matn);
1508 24542 : ex = valp(x);
1509 24542 : if (varncmp(vx, v) < 0)
1510 : {
1511 49 : if (lx == 2) return matn >= 0? scalarmat(x, matn): gcopy(x);
1512 49 : av = avma; X = pol_x(vx);
1513 49 : av2 = avma;
1514 49 : z = gadd(gsubst(gel(x,lx-1),v,y), zeroser(vx,1));
1515 189 : for (i = lx-2; i>=2; i--)
1516 : {
1517 140 : z = gadd(gmul(z,X), gsubst(gel(x,i),v,y));
1518 140 : if (gc_needed(av2,1))
1519 : {
1520 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1521 0 : z = gerepileupto(av2, z);
1522 : }
1523 : }
1524 49 : if (ex) z = gmul(z, pol_xnall(ex,vx));
1525 49 : return gerepileupto(av, z);
1526 : }
1527 24493 : switch(ty) /* here vx == v */
1528 : {
1529 20293 : case t_SER:
1530 20293 : vy = varn(y); ey = valp(y);
1531 20293 : if (ey < 1 || lx == 2) return zeroser(vy, ey*(ex+lx-2));
1532 20293 : if (ey == 1 && serequalXk(y)
1533 1953 : && (varncmp(vx,vy) >= 0 || varncmp(gvar2(x), vy) >= 0))
1534 : { /* y = t + O(t^N) */
1535 1953 : if (lx > ly)
1536 : { /* correct number of significant terms */
1537 1624 : l = ly;
1538 1624 : if (!ex)
1539 1603 : for (i = 3; i < lx; i++)
1540 1603 : if (++l >= lx || !gequal0(gel(x,i))) break;
1541 1624 : lx = l;
1542 : }
1543 1953 : z = sertrunc_copy(x, lx - 2); if (vx != vy) setvarn(z,vy);
1544 1953 : return z;
1545 : }
1546 18340 : if (vy != vx)
1547 : {
1548 28 : long nx = lx - 2, n = minss(ey * nx, ly - 2);
1549 28 : av = avma; z = gel(x, nx+1);
1550 91 : for (i = nx; i > 1; i--)
1551 : {
1552 63 : z = gadd(gmul(y,z), gel(x,i));
1553 63 : if (gc_needed(av,1))
1554 : {
1555 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst (i = %ld)", i);
1556 0 : z = gerepileupto(av, z);
1557 : }
1558 : }
1559 28 : if (ex)
1560 : {
1561 21 : if (ex < 0) { y = ginv(y); ex = -ex; }
1562 21 : z = gmul(z, gpowgs(sertrunc(y, n), ex));
1563 : }
1564 28 : if (lg(z)-2 > n) z = sertrunc_copy(z, n);
1565 28 : return gerepileupto(av,z);
1566 : }
1567 18312 : l = (lx-2)*ey + 2;
1568 18312 : if (ex) { if (l>ly) l = ly; }
1569 18263 : else if (lx != 3)
1570 : {
1571 18277 : for (i = 3; i < lx; i++)
1572 18277 : if (!isexactzero(gel(x,i))) break;
1573 18263 : l = minss(l, (i-2)*ey + (gequal0(y)? 2 : ly));
1574 : }
1575 18312 : av = avma; t = leafcopy(y);
1576 18312 : if (l < ly) setlg(t, l);
1577 18312 : z = scalarser(gen_1, varn(y), l-2);
1578 18312 : gel(z,2) = gel(x,2); /* ensure lg(z) = l even if x[2] = 0 */
1579 71260 : for (i = 3, jb = ey; jb <= l-2; i++,jb += ey)
1580 : {
1581 52955 : if (i < lx) {
1582 114793 : for (j = jb+2; j < minss(l, jb+ly); j++)
1583 61901 : gel(z,j) = gadd(gel(z,j), gmul(gel(x,i),gel(t,j-jb)));
1584 : }
1585 79772 : for (j = minss(ly-1, l-1-jb-ey); j > 1; j--)
1586 : {
1587 26824 : GEN a = gmul(gel(t,2), gel(y,j));
1588 63805 : for (k=2; k<j; k++) a = gadd(a, gmul(gel(t,j-k+2), gel(y,k)));
1589 26824 : gel(t,j) = a;
1590 : }
1591 52948 : if (gc_needed(av,1))
1592 : {
1593 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gsubst");
1594 0 : gerepileall(av,2, &z,&t);
1595 : }
1596 : }
1597 18305 : if (!ex) return gerepilecopy(av,z);
1598 49 : if (ex < 0) { ex = -ex; y = ginv(y); }
1599 49 : return gerepileupto(av, gmul(z, gpowgs(sertrunc(y, l-2), ex)));
1600 :
1601 4158 : case t_POL: case t_RFRAC:
1602 : {
1603 4158 : long N, n = lx-2;
1604 4158 : vy = gvar(y); ey = gval(y,vy);
1605 4158 : if (ey == LONG_MAX)
1606 : { /* y = 0 */
1607 49 : if (ex < 0) pari_err_INV("gsubst",y);
1608 35 : if (!n) return gcopy(x);
1609 28 : if (ex > 0) return Rg_get_0(ty == t_RFRAC? gel(y,2): y);
1610 14 : y = Rg_get_1(ty == t_RFRAC? gel(y,2): y);
1611 14 : return gmul(y, gel(x,2));
1612 : }
1613 4109 : if (ey < 1 || n == 0) return zeroser(vy, ey*(ex+n));
1614 4102 : av = avma;
1615 4102 : n *= ey;
1616 4102 : N = ex? n: maxss(n-ey,1);
1617 4102 : y = (ty == t_RFRAC)? rfrac_to_ser_i(y, N+2): RgX_to_ser(y, N+2);
1618 4102 : if (lg(y)-2 > n) setlg(y, n+2);
1619 4102 : x = ser2pol_i(x, lx);
1620 4102 : if (varncmp(vy,vx) > 0)
1621 42 : z = gadd(poleval(x, y), zeroser(vy,n));
1622 : else
1623 : {
1624 4060 : z = RgXn_eval(x, ser2rfrac_i(y), n);
1625 4060 : if (varn(z) == vy) z = RgX_to_ser(z, n+2);
1626 : }
1627 4102 : switch(typ(z))
1628 : {
1629 4102 : case t_SER:
1630 : case t_POL:
1631 4102 : if (varncmp(varn(z),vy) <= 0) break;
1632 0 : default: z = scalarser(z, vy, n);
1633 : }
1634 4102 : if (!ex) return gerepilecopy(av, z);
1635 4039 : return gerepileupto(av, gmul(z, gpowgs(y,ex)));
1636 : }
1637 :
1638 42 : default:
1639 42 : if (isexactzero(y))
1640 : {
1641 35 : if (ex < 0) pari_err_INV("gsubst",y);
1642 14 : if (ex > 0) return gcopy(y);
1643 7 : if (lx > 2) return gadd(gel(x,2), y); /*add maps to correct ring*/
1644 : }
1645 7 : pari_err_TYPE2("substitution",x,y);
1646 : }
1647 0 : break;
1648 :
1649 1246 : case t_RFRAC:
1650 : {
1651 1246 : GEN a = gel(x,1), b = gel(x,2);
1652 1246 : av = avma;
1653 1246 : a = gsubst(a, v, y);
1654 1246 : b = gsubst(b, v, y); return gerepileupto(av, gdiv(a, b));
1655 : }
1656 :
1657 659849 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gsubst(gel(x,i),v,y));
1658 56 : case t_LIST:
1659 56 : z = mklist();
1660 56 : list_data(z) = list_data(x)? gsubst(list_data(x),v,y): NULL;
1661 56 : return z;
1662 : }
1663 0 : return gcopy(x);
1664 : }
1665 :
1666 : /* Return P(x * h), not memory clean */
1667 : GEN
1668 4193 : ser_unscale(GEN P, GEN h)
1669 : {
1670 4193 : long l = lg(P);
1671 4193 : GEN Q = cgetg(l,t_SER);
1672 4193 : Q[1] = P[1];
1673 4193 : if (l != 2)
1674 : {
1675 4193 : long i = 2;
1676 4193 : GEN hi = gpowgs(h, valp(P));
1677 4193 : gel(Q,i) = gmul(gel(P,i), hi);
1678 200508 : for (i++; i<l; i++) { hi = gmul(hi,h); gel(Q,i) = gmul(gel(P,i), hi); }
1679 : }
1680 4193 : return Q;
1681 : }
1682 :
1683 : GEN
1684 959 : gsubstvec(GEN e, GEN v, GEN r)
1685 : {
1686 959 : pari_sp av = avma;
1687 959 : long i, j, k, l = lg(v);
1688 : GEN w, z, R;
1689 959 : if ( !is_vec_t(typ(v)) ) pari_err_TYPE("substvec",v);
1690 959 : if ( !is_vec_t(typ(r)) ) pari_err_TYPE("substvec",r);
1691 959 : if (lg(r)!=l) pari_err_DIM("substvec");
1692 959 : w = cgetg(l, t_VECSMALL);
1693 959 : z = cgetg(l, t_VECSMALL);
1694 959 : R = cgetg(l, t_VEC); k = 0;
1695 4235 : for(i = j = 1; i < l; i++)
1696 : {
1697 3276 : GEN T = gel(v,i), ri = gel(r,i);
1698 3276 : if (!gequalX(T)) pari_err_TYPE("substvec [not a variable]", T);
1699 3276 : if (gvar(ri) == NO_VARIABLE) /* no need to take precautions */
1700 : {
1701 1855 : e = gsubst(e, varn(T), ri);
1702 1855 : if (is_vec_t(typ(ri)) && k++) e = shallowconcat1(e);
1703 : }
1704 : else
1705 : {
1706 1421 : w[j] = varn(T);
1707 1421 : z[j] = fetch_var();
1708 1421 : gel(R,j) = ri; j++;
1709 : }
1710 : }
1711 2380 : for(i = 1; i < j; i++) e = gsubst(e,w[i],pol_x(z[i]));
1712 2380 : for(i = 1; i < j; i++)
1713 : {
1714 1421 : e = gsubst(e,z[i],gel(R,i));
1715 1421 : if (is_vec_t(typ(gel(R,i))) && k++) e = shallowconcat1(e);
1716 : }
1717 2380 : for(i = 1; i < j; i++) (void)delete_var();
1718 959 : return k > 1? gerepilecopy(av, e): gerepileupto(av, e);
1719 : }
1720 :
1721 : /*******************************************************************/
1722 : /* */
1723 : /* SERIE RECIPROQUE D'UNE SERIE */
1724 : /* */
1725 : /*******************************************************************/
1726 :
1727 : GEN
1728 98 : serreverse(GEN x)
1729 : {
1730 98 : long v=varn(x), lx = lg(x), i, mi;
1731 98 : pari_sp av0 = avma, av;
1732 : GEN a, y, u;
1733 :
1734 98 : if (typ(x)!=t_SER) pari_err_TYPE("serreverse",x);
1735 98 : if (valp(x)!=1) pari_err_DOMAIN("serreverse", "valuation", "!=", gen_1,x);
1736 91 : if (lx < 3) pari_err_DOMAIN("serreverse", "x", "=", gen_0,x);
1737 91 : y = ser_normalize(x);
1738 91 : if (y == x) a = NULL; else { a = gel(x,2); x = y; }
1739 91 : av = avma;
1740 252 : mi = lx-1; while (mi>=3 && gequal0(gel(x,mi))) mi--;
1741 91 : u = cgetg(lx,t_SER);
1742 91 : y = cgetg(lx,t_SER);
1743 91 : u[1] = y[1] = evalsigne(1) | _evalvalp(1) | evalvarn(v);
1744 91 : gel(u,2) = gel(y,2) = gen_1;
1745 91 : if (lx > 3)
1746 : {
1747 84 : gel(u,3) = gmulsg(-2,gel(x,3));
1748 84 : gel(y,3) = gneg(gel(x,3));
1749 : }
1750 1113 : for (i=3; i<lx-1; )
1751 : {
1752 : pari_sp av2;
1753 : GEN p1;
1754 1022 : long j, k, K = minss(i,mi);
1755 8456 : for (j=3; j<i+1; j++)
1756 : {
1757 7434 : av2 = avma; p1 = gel(x,j);
1758 39291 : for (k = maxss(3,j+2-mi); k < j; k++)
1759 31857 : p1 = gadd(p1, gmul(gel(u,k),gel(x,j-k+2)));
1760 7434 : p1 = gneg(p1);
1761 7434 : gel(u,j) = gerepileupto(av2, gadd(gel(u,j), p1));
1762 : }
1763 1022 : av2 = avma;
1764 1022 : p1 = gmulsg(i,gel(x,i+1));
1765 8309 : for (k = 2; k < K; k++)
1766 : {
1767 7287 : GEN p2 = gmul(gel(x,k+1),gel(u,i-k+2));
1768 7287 : p1 = gadd(p1, gmulsg(k,p2));
1769 : }
1770 1022 : i++;
1771 1022 : gel(u,i) = gerepileupto(av2, gneg(p1));
1772 1022 : gel(y,i) = gdivgu(gel(u,i), i-1);
1773 1022 : if (gc_needed(av,2))
1774 : {
1775 0 : GEN dummy = cgetg(1,t_VEC);
1776 0 : if(DEBUGMEM>1) pari_warn(warnmem,"serreverse");
1777 0 : for(k=i+1; k<lx; k++) gel(u,k) = gel(y,k) = dummy;
1778 0 : gerepileall(av,2, &u,&y);
1779 : }
1780 : }
1781 91 : if (a) y = ser_unscale(y, ginv(a));
1782 91 : return gerepilecopy(av0,y);
1783 : }
1784 :
1785 : /*******************************************************************/
1786 : /* */
1787 : /* DERIVATION ET INTEGRATION */
1788 : /* */
1789 : /*******************************************************************/
1790 : GEN
1791 23499 : derivser(GEN x)
1792 : {
1793 23499 : long i, vx = varn(x), e = valp(x), lx = lg(x);
1794 : GEN y;
1795 23499 : if (ser_isexactzero(x))
1796 : {
1797 7 : x = gcopy(x);
1798 7 : if (e) setvalp(x,e-1);
1799 7 : return x;
1800 : }
1801 23492 : if (e)
1802 : {
1803 602 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|evalvalp(e-1) | evalvarn(vx);
1804 22960 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i+e-2,gel(x,i));
1805 : } else {
1806 22890 : if (lx == 3) return zeroser(vx, 0);
1807 19964 : lx--;
1808 19964 : y = cgetg(lx,t_SER); y[1] = evalsigne(1)|_evalvalp(0) | evalvarn(vx);
1809 62461 : for (i=2; i<lx; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1810 : }
1811 20566 : return normalizeser(y);
1812 : }
1813 :
1814 : static GEN
1815 56 : rfrac_deriv(GEN x, long v)
1816 : {
1817 56 : pari_sp av = avma;
1818 56 : GEN y = cgetg(3,t_RFRAC), a = gel(x,1), b = gel(x,2), bp, b0, t, T;
1819 56 : long vx = varn(b);
1820 :
1821 56 : bp = deriv(b, v);
1822 56 : t = simplify_shallow(RgX_gcd(bp, b));
1823 56 : if (typ(t) != t_POL || varn(t) != vx)
1824 : {
1825 35 : if (gequal1(t)) b0 = b;
1826 : else
1827 : {
1828 0 : b0 = RgX_Rg_div(b, t);
1829 0 : bp = RgX_Rg_div(bp, t);
1830 : }
1831 35 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1832 35 : if (isexactzero(a)) return gerepileupto(av, a);
1833 35 : if (b0 == b)
1834 : {
1835 35 : gel(y,1) = gerepileupto((pari_sp)y, a);
1836 35 : gel(y,2) = RgX_sqr(b);
1837 : }
1838 : else
1839 : {
1840 0 : gel(y,1) = a;
1841 0 : gel(y,2) = RgX_Rg_mul(RgX_sqr(b0), t);
1842 0 : y = gerepilecopy(av, y);
1843 : }
1844 35 : return y;
1845 : }
1846 21 : b0 = gdivexact(b, t);
1847 21 : bp = gdivexact(bp,t);
1848 21 : a = gsub(gmul(b0, deriv(a,v)), gmul(a, bp));
1849 21 : if (isexactzero(a)) return gerepileupto(av, a);
1850 14 : T = RgX_gcd(a, t);
1851 14 : if (typ(T) != t_POL || varn(T) != vx)
1852 : {
1853 0 : a = gdiv(a, T);
1854 0 : t = gdiv(t, T);
1855 : }
1856 14 : else if (!gequal1(T))
1857 : {
1858 0 : a = gdivexact(a, T);
1859 0 : t = gdivexact(t, T);
1860 : }
1861 14 : gel(y,1) = a;
1862 14 : gel(y,2) = gmul(RgX_sqr(b0), t);
1863 14 : return gerepilecopy(av, y);
1864 : }
1865 :
1866 : GEN
1867 110054 : deriv(GEN x, long v)
1868 : {
1869 110054 : long tx = typ(x);
1870 110054 : if (is_const_t(tx))
1871 39067 : switch(tx)
1872 : {
1873 14 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
1874 14 : case t_FFELT: return FF_zero(x);
1875 39039 : default: return gen_0;
1876 : }
1877 70987 : if (v < 0)
1878 : {
1879 49 : if (tx == t_CLOSURE) return closure_deriv(x);
1880 49 : v = gvar9(x);
1881 : }
1882 70987 : switch(tx)
1883 : {
1884 14 : case t_POLMOD:
1885 : {
1886 14 : GEN a = gel(x,2), b = gel(x,1);
1887 14 : if (v == varn(b)) return Rg_get_0(b);
1888 7 : retmkpolmod(deriv(a,v), RgX_copy(b));
1889 : }
1890 70728 : case t_POL:
1891 70728 : switch(varncmp(varn(x), v))
1892 : {
1893 0 : case 1: return Rg_get_0(x);
1894 63315 : case 0: return RgX_deriv(x);
1895 : }
1896 109431 : pari_APPLY_pol(deriv(gel(x,i),v));
1897 :
1898 147 : case t_SER:
1899 147 : switch(varncmp(varn(x), v))
1900 : {
1901 0 : case 1: return Rg_get_0(x);
1902 133 : case 0: return derivser(x);
1903 : }
1904 14 : if (ser_isexactzero(x)) return gcopy(x);
1905 28 : pari_APPLY_ser(deriv(gel(x,i),v));
1906 :
1907 56 : case t_RFRAC:
1908 56 : return rfrac_deriv(x,v);
1909 :
1910 42 : case t_VEC: case t_COL: case t_MAT:
1911 84 : pari_APPLY_same(deriv(gel(x,i),v));
1912 : }
1913 0 : pari_err_TYPE("deriv",x);
1914 : return NULL; /* LCOV_EXCL_LINE */
1915 : }
1916 :
1917 : /* n-th derivative of t_SER x, n > 0 */
1918 : static GEN
1919 189 : derivnser(GEN x, long n)
1920 : {
1921 189 : long i, vx = varn(x), e = valp(x), lx = lg(x);
1922 : GEN y;
1923 189 : if (ser_isexactzero(x))
1924 : {
1925 7 : x = gcopy(x);
1926 7 : if (e) setvalp(x,e-n);
1927 7 : return x;
1928 : }
1929 336 : if (e < 0 || e >= n)
1930 : {
1931 154 : y = cgetg(lx,t_SER);
1932 154 : y[1] = evalsigne(1)| evalvalp(e-n) | evalvarn(vx);
1933 742 : for (i=0; i<lx-2; i++)
1934 588 : gel(y,i+2) = gmul(muls_interval(i+e-n+1,i+e), gel(x,i+2));
1935 : } else {
1936 28 : if (lx <= n+2) return zeroser(vx, 0);
1937 28 : lx -= n;
1938 28 : y = cgetg(lx,t_SER);
1939 28 : y[1] = evalsigne(1)|_evalvalp(0) | evalvarn(vx);
1940 105 : for (i=0; i<lx-2; i++)
1941 77 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n-e));
1942 : }
1943 182 : return normalizeser(y);
1944 : }
1945 :
1946 : /* n-th derivative of t_POL x, n > 0 */
1947 : static GEN
1948 826 : RgX_derivn(GEN x, long n)
1949 : {
1950 826 : long i, vx = varn(x), lx = lg(x);
1951 : GEN y;
1952 826 : if (lx <= n+2) return pol_0(vx);
1953 742 : lx -= n;
1954 742 : y = cgetg(lx,t_POL);
1955 742 : y[1] = evalsigne(1)| evalvarn(vx);
1956 36883 : for (i=0; i<lx-2; i++)
1957 36141 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n));
1958 742 : return normalizepol_lg(y, lx);
1959 : }
1960 :
1961 : static GEN
1962 42 : rfrac_derivn(GEN x, long n, long vs)
1963 : {
1964 42 : pari_sp av = avma;
1965 42 : GEN u = gel(x,1), v = gel(x,2);
1966 42 : GEN dv = deriv(v, vs);
1967 : long i;
1968 112 : for (i=1; i<=n; i++)
1969 : {
1970 70 : GEN du = deriv(u, vs);
1971 70 : u = gadd(gmul(du, v), gmulsg (-i, gmul(dv, u)));
1972 : }
1973 42 : v = gpowgs(v, n+1);
1974 42 : return gerepileupto(av, gdiv(u, v));
1975 : }
1976 :
1977 : GEN
1978 1344 : derivn(GEN x, long n, long v)
1979 : {
1980 : long tx;
1981 1344 : if (n < 0) pari_err_DOMAIN("derivn","n","<", gen_0, stoi(n));
1982 1337 : if (n == 0) return gcopy(x);
1983 1337 : tx = typ(x);
1984 1337 : if (is_const_t(tx))
1985 49 : switch(tx)
1986 : {
1987 21 : case t_INTMOD: retmkintmod(gen_0, icopy(gel(x,1)));
1988 21 : case t_FFELT: return FF_zero(x);
1989 7 : default: return gen_0;
1990 : }
1991 1288 : if (v < 0)
1992 : {
1993 1050 : if (tx == t_CLOSURE) return closure_derivn(x, n);
1994 945 : v = gvar9(x);
1995 : }
1996 1183 : switch(tx)
1997 : {
1998 21 : case t_POLMOD:
1999 : {
2000 21 : GEN a = gel(x,2), b = gel(x,1);
2001 21 : if (v == varn(b)) return Rg_get_0(b);
2002 14 : retmkpolmod(derivn(a,n,v), RgX_copy(b));
2003 : }
2004 854 : case t_POL:
2005 854 : switch(varncmp(varn(x), v))
2006 : {
2007 0 : case 1: return Rg_get_0(x);
2008 826 : case 0: return RgX_derivn(x,n);
2009 : }
2010 84 : pari_APPLY_pol(derivn(gel(x,i),n,v));
2011 :
2012 196 : case t_SER:
2013 196 : switch(varncmp(varn(x), v))
2014 : {
2015 0 : case 1: return Rg_get_0(x);
2016 189 : case 0: return derivnser(x, n);
2017 : }
2018 7 : if (ser_isexactzero(x)) return gcopy(x);
2019 28 : pari_APPLY_ser(derivn(gel(x,i),n,v));
2020 :
2021 42 : case t_RFRAC:
2022 42 : return rfrac_derivn(x, n, v);
2023 :
2024 63 : case t_VEC: case t_COL: case t_MAT:
2025 126 : pari_APPLY_same(derivn(gel(x,i),n,v));
2026 : }
2027 7 : pari_err_TYPE("derivn",x);
2028 : return NULL; /* LCOV_EXCL_LINE */
2029 : }
2030 :
2031 : static long
2032 833 : lookup(GEN v, long vx)
2033 : {
2034 833 : long i ,l = lg(v);
2035 1491 : for(i=1; i<l; i++)
2036 1253 : if (varn(gel(v,i)) == vx) return i;
2037 238 : return 0;
2038 : }
2039 :
2040 : GEN
2041 3535 : diffop(GEN x, GEN v, GEN dv)
2042 : {
2043 : pari_sp av;
2044 3535 : long i, idx, lx, tx = typ(x), vx;
2045 : GEN y;
2046 3535 : if (!is_vec_t(typ(v))) pari_err_TYPE("diffop",v);
2047 3535 : if (!is_vec_t(typ(dv))) pari_err_TYPE("diffop",dv);
2048 3535 : if (lg(v)!=lg(dv)) pari_err_DIM("diffop");
2049 3535 : if (is_const_t(tx)) return gen_0;
2050 1148 : switch(tx)
2051 : {
2052 84 : case t_POLMOD:
2053 84 : av = avma;
2054 84 : vx = varn(gel(x,1)); idx = lookup(v,vx);
2055 84 : if (idx) /*Assume the users now what they are doing */
2056 0 : y = gmodulo(diffop(gel(x,2),v,dv), gel(x,1));
2057 : else
2058 : {
2059 84 : GEN m = gel(x,1), pol=gel(x,2);
2060 84 : GEN u = gneg(gdiv(diffop(m,v,dv),RgX_deriv(m)));
2061 84 : y = diffop(pol,v,dv);
2062 84 : if (typ(pol)==t_POL && varn(pol)==varn(m))
2063 70 : y = gadd(y, gmul(u,RgX_deriv(pol)));
2064 84 : y = gmodulo(y, gel(x,1));
2065 : }
2066 84 : return gerepileupto(av, y);
2067 952 : case t_POL:
2068 952 : if (signe(x)==0) return gen_0;
2069 742 : vx = varn(x); idx = lookup(v,vx);
2070 742 : av = avma; lx = lg(x);
2071 742 : y = diffop(gel(x,lx-1),v,dv);
2072 2842 : for (i=lx-2; i>=2; i--) y = gadd(gmul(y,pol_x(vx)),diffop(gel(x,i),v,dv));
2073 742 : if (idx) y = gadd(y, gmul(gel(dv,idx),RgX_deriv(x)));
2074 742 : return gerepileupto(av, y);
2075 :
2076 7 : case t_SER:
2077 7 : if (signe(x)==0) return gen_0;
2078 7 : vx = varn(x); idx = lookup(v,vx);
2079 7 : if (!idx) return gen_0;
2080 7 : av = avma;
2081 7 : if (ser_isexactzero(x)) y = x;
2082 : else
2083 : {
2084 7 : y = cgetg_copy(x, &lx); y[1] = x[1];
2085 119 : for (i=2; i<lx; i++) gel(y,i) = diffop(gel(x,i),v,dv);
2086 7 : y = normalizeser(y); /* y is probably invalid */
2087 7 : y = gsubst(y, vx, pol_x(vx)); /* Fix that */
2088 : }
2089 7 : y = gadd(y, gmul(gel(dv,idx),derivser(x)));
2090 7 : return gerepileupto(av, y);
2091 :
2092 105 : case t_RFRAC: {
2093 105 : GEN a = gel(x,1), b = gel(x,2), ap, bp;
2094 105 : av = avma;
2095 105 : ap = diffop(a, v, dv); bp = diffop(b, v, dv);
2096 105 : y = gsub(gdiv(ap,b),gdiv(gmul(a,bp),gsqr(b)));
2097 105 : return gerepileupto(av, y);
2098 : }
2099 :
2100 0 : case t_VEC: case t_COL: case t_MAT:
2101 0 : pari_APPLY_same(diffop(gel(x,i),v,dv));
2102 : }
2103 0 : pari_err_TYPE("diffop",x);
2104 : return NULL; /* LCOV_EXCL_LINE */
2105 : }
2106 :
2107 : GEN
2108 42 : diffop0(GEN x, GEN v, GEN dv, long n)
2109 : {
2110 42 : pari_sp av=avma;
2111 : long i;
2112 245 : for(i=1; i<=n; i++)
2113 203 : x = gerepileupto(av, diffop(x,v,dv));
2114 42 : return x;
2115 : }
2116 :
2117 : /********************************************************************/
2118 : /** **/
2119 : /** TAYLOR SERIES **/
2120 : /** **/
2121 : /********************************************************************/
2122 : /* swap vars (vx,v) in x (assume vx < v, vx main variable in x), then call
2123 : * act(data, v, x). FIXME: use in other places */
2124 : static GEN
2125 21 : swapvar_act(GEN x, long vx, long v, GEN (*act)(void*, long, GEN), void *data)
2126 : {
2127 21 : long v0 = fetch_var();
2128 21 : GEN y = act(data, v, gsubst(x,vx,pol_x(v0)));
2129 14 : y = gsubst(y,v0,pol_x(vx));
2130 14 : (void)delete_var(); return y;
2131 : }
2132 : /* x + O(v^data) */
2133 : static GEN
2134 7 : tayl_act(void *data, long v, GEN x) { return gadd(zeroser(v, (long)data), x); }
2135 : static GEN
2136 14 : integ_act(void *data, long v, GEN x) { (void)data; return integ(x,v); }
2137 :
2138 : GEN
2139 7 : tayl(GEN x, long v, long precS)
2140 : {
2141 7 : long vx = gvar9(x);
2142 : pari_sp av;
2143 :
2144 7 : if (varncmp(v, vx) <= 0) return gadd(zeroser(v,precS), x);
2145 7 : av = avma;
2146 7 : return gerepileupto(av, swapvar_act(x, vx, v, tayl_act, (void*)precS));
2147 : }
2148 :
2149 : GEN
2150 6426 : ggrando(GEN x, long n)
2151 : {
2152 : long m, v;
2153 :
2154 6426 : switch(typ(x))
2155 : {
2156 3633 : case t_INT:/* bug 3 + O(1) */
2157 3633 : if (signe(x) <= 0) pari_err_DOMAIN("O", "x", "<=", gen_0, x);
2158 3633 : if (!is_pm1(x)) return zeropadic(x,n);
2159 : /* +/-1 = x^0 */
2160 91 : v = m = 0; break;
2161 2786 : case t_POL:
2162 2786 : if (!signe(x)) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2163 2786 : v = varn(x);
2164 2786 : m = n * RgX_val(x); break;
2165 7 : case t_RFRAC:
2166 7 : if (gequal0(gel(x,1))) pari_err_DOMAIN("O", "x", "=", gen_0, x);
2167 7 : v = gvar(x);
2168 7 : m = n * gval(x,v); break;
2169 0 : default: pari_err_TYPE("O", x);
2170 : v = m = 0; /* LCOV_EXCL_LINE */
2171 : }
2172 2884 : return zeroser(v,m);
2173 : }
2174 :
2175 : /*******************************************************************/
2176 : /* */
2177 : /* FORMAL INTEGRATION */
2178 : /* */
2179 : /*******************************************************************/
2180 : GEN
2181 105 : RgX_integ(GEN x)
2182 : {
2183 105 : long i, lx = lg(x);
2184 : GEN y;
2185 105 : if (lx == 2) return RgX_copy(x);
2186 91 : y = cgetg(lx+1, t_POL); y[1] = x[1]; gel(y,2) = gen_0;
2187 273 : for (i=3; i<=lx; i++) gel(y,i) = gdivgu(gel(x,i-1),i-2);
2188 91 : return y;
2189 : }
2190 :
2191 : static void
2192 35 : err_intformal(GEN x)
2193 35 : { pari_err_DOMAIN("intformal", "residue(series, pole)", "!=", gen_0, x); }
2194 :
2195 : GEN
2196 24122 : integser(GEN x)
2197 : {
2198 24122 : long i, lx = lg(x), vx = varn(x), e = valp(x);
2199 : GEN y;
2200 24122 : if (lx == 2) return zeroser(vx, e+1);
2201 21175 : y = cgetg(lx, t_SER);
2202 91434 : for (i=2; i<lx; i++)
2203 : {
2204 70266 : long j = i+e-1;
2205 70266 : GEN c = gel(x,i);
2206 70266 : if (j)
2207 69951 : c = gdivgs(c, j);
2208 : else
2209 : { /* should be isexactzero, but try to avoid error */
2210 315 : if (!gequal0(c)) err_intformal(x);
2211 308 : c = gen_0;
2212 : }
2213 70259 : gel(y,i) = c;
2214 : }
2215 21168 : y[1] = evalsigne(1) | evalvarn(vx) | evalvalp(e+1); return y;
2216 : }
2217 :
2218 : GEN
2219 350 : integ(GEN x, long v)
2220 : {
2221 350 : long tx = typ(x), vx, n;
2222 350 : pari_sp av = avma;
2223 : GEN y, p1;
2224 :
2225 350 : if (v < 0) { v = gvar9(x); if (v == NO_VARIABLE) v = 0; }
2226 350 : if (is_scalar_t(tx))
2227 : {
2228 63 : if (tx == t_POLMOD)
2229 : {
2230 14 : GEN a = gel(x,2), b = gel(x,1);
2231 14 : vx = varn(b);
2232 14 : if (varncmp(v, vx) > 0) retmkpolmod(integ(a,v), RgX_copy(b));
2233 7 : if (v == vx) pari_err_PRIORITY("intformal",x,"=",v);
2234 : }
2235 49 : return deg1pol(x, gen_0, v);
2236 : }
2237 :
2238 287 : switch(tx)
2239 : {
2240 112 : case t_POL:
2241 112 : vx = varn(x);
2242 112 : if (v == vx) return RgX_integ(x);
2243 42 : if (lg(x) == 2) {
2244 14 : if (varncmp(vx, v) < 0) v = vx;
2245 14 : return zeropol(v);
2246 : }
2247 28 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2248 84 : pari_APPLY_pol(integ(gel(x,i),v));
2249 :
2250 77 : case t_SER:
2251 77 : vx = varn(x);
2252 77 : if (v == vx) return integser(x);
2253 21 : if (lg(x) == 2) {
2254 14 : if (varncmp(vx, v) < 0) v = vx;
2255 14 : return zeroser(v, valp(x));
2256 : }
2257 7 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2258 28 : pari_APPLY_ser(integ(gel(x,i),v));
2259 :
2260 56 : case t_RFRAC:
2261 : {
2262 56 : GEN a = gel(x,1), b = gel(x,2), c, d, s;
2263 56 : vx = varn(b);
2264 56 : if (varncmp(vx, v) > 0) return deg1pol(x, gen_0, v);
2265 49 : if (varncmp(vx, v) < 0)
2266 14 : return gerepileupto(av, swapvar_act(x, vx, v, integ_act, NULL));
2267 :
2268 35 : n = degpol(b);
2269 35 : if (typ(a) == t_POL && varn(a) == vx) n += degpol(a);
2270 35 : y = integ(gadd(x, zeroser(v,n + 2)), v);
2271 35 : y = gdiv(gtrunc(gmul(b, y)), b);
2272 35 : if (typ(y) != t_RFRAC) pari_err_BUG("intformal(t_RFRAC)");
2273 35 : c = gel(y,1); d = gel(y,2);
2274 35 : s = gsub(gmul(deriv(c,v),d), gmul(c,deriv(d,v)));
2275 : /* (c'd-cd')/d^2 = y' = x = a/b ? */
2276 35 : if (!gequal(gmul(s,b), gmul(a,gsqr(d)))) err_intformal(x);
2277 7 : if (typ(y)==t_RFRAC && lg(gel(y,1)) == lg(gel(y,2)))
2278 : {
2279 7 : GEN p2 = leading_coeff(gel(y,2));
2280 7 : p1 = gel(y,1);
2281 7 : if (typ(p1) == t_POL && varn(p1) == vx) p1 = leading_coeff(p1);
2282 7 : y = gsub(y, gdiv(p1,p2));
2283 : }
2284 7 : return gerepileupto(av,y);
2285 : }
2286 :
2287 42 : case t_VEC: case t_COL: case t_MAT:
2288 84 : pari_APPLY_same(integ(gel(x,i),v));
2289 : }
2290 0 : pari_err_TYPE("integ",x);
2291 : return NULL; /* LCOV_EXCL_LINE */
2292 : }
2293 :
2294 : /*******************************************************************/
2295 : /* */
2296 : /* FLOOR */
2297 : /* */
2298 : /*******************************************************************/
2299 : static GEN
2300 518 : quad_floor(GEN x)
2301 : {
2302 518 : GEN Q = gel(x,1), D = quad_disc(x), u, v, b, d, z;
2303 518 : if (signe(D) < 0) return NULL;
2304 490 : x = Q_remove_denom(x, &d);
2305 490 : u = gel(x,2);
2306 490 : v = gel(x,3); b = gel(Q,3);
2307 490 : if (typ(u) != t_INT || typ(v) != t_INT) return NULL;
2308 : /* x0 = (2u + v*(-b + sqrt(D))) / (2d) */
2309 483 : z = sqrti(mulii(D, sqri(v)));
2310 483 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2311 : /* z = floor(v * sqrt(D)) */
2312 483 : z = addii(subii(shifti(u,1), mulii(v,b)), z);
2313 483 : return truedivii(z, d? shifti(d,1): gen_2);
2314 : }
2315 : GEN
2316 5285274 : gfloor(GEN x)
2317 : {
2318 5285274 : switch(typ(x))
2319 : {
2320 5262173 : case t_INT: return icopy(x);
2321 35 : case t_POL: return RgX_copy(x);
2322 10648 : case t_REAL: return floorr(x);
2323 11592 : case t_FRAC: return truedivii(gel(x,1),gel(x,2));
2324 511 : case t_QUAD:
2325 : {
2326 511 : pari_sp av = avma;
2327 : GEN y;
2328 511 : if (!(y = quad_floor(x))) break;
2329 476 : return gerepileuptoint(av, y);
2330 : }
2331 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2332 98 : case t_VEC: case t_COL: case t_MAT:
2333 1533 : pari_APPLY_same(gfloor(gel(x,i)));
2334 : }
2335 238 : pari_err_TYPE("gfloor",x);
2336 : return NULL; /* LCOV_EXCL_LINE */
2337 : }
2338 :
2339 : GEN
2340 231408 : gfrac(GEN x)
2341 : {
2342 : pari_sp av;
2343 : GEN y;
2344 231408 : switch(typ(x))
2345 : {
2346 21101 : case t_INT: return gen_0;
2347 7 : case t_POL: return pol_0(varn(x));
2348 177175 : case t_REAL: av = avma; return gerepileuptoleaf(av, subri(x, floorr(x)));
2349 29668 : case t_FRAC: retmkfrac(modii(gel(x,1),gel(x,2)), icopy(gel(x,2)));
2350 7 : case t_QUAD:
2351 7 : av = avma; if (!(y = quad_floor(x))) break;
2352 7 : return gerepileupto(av, gsub(x, y));
2353 7 : case t_RFRAC: retmkrfrac(grem(gel(x,1),gel(x,2)), gcopy(gel(x,2)));
2354 3415 : case t_VEC: case t_COL: case t_MAT:
2355 32869 : pari_APPLY_same(gfrac(gel(x,i)));
2356 : }
2357 28 : pari_err_TYPE("gfrac",x);
2358 : return NULL; /* LCOV_EXCL_LINE */
2359 : }
2360 :
2361 : /* assume x t_REAL */
2362 : GEN
2363 3339 : ceilr(GEN x)
2364 : {
2365 3339 : pari_sp av = avma;
2366 3339 : GEN y = floorr(x);
2367 3339 : if (cmpri(x, y)) return gerepileuptoint(av, addui(1,y));
2368 28 : return y;
2369 : }
2370 :
2371 : GEN
2372 31206 : gceil(GEN x)
2373 : {
2374 : pari_sp av;
2375 : GEN y;
2376 31206 : switch(typ(x))
2377 : {
2378 17481 : case t_INT: return icopy(x);
2379 21 : case t_POL: return RgX_copy(x);
2380 3261 : case t_REAL: return ceilr(x);
2381 10317 : case t_FRAC:
2382 10317 : av = avma; y = divii(gel(x,1),gel(x,2));
2383 10317 : if (signe(gel(x,1)) > 0) y = gerepileuptoint(av, addui(1,y));
2384 10317 : return y;
2385 49 : case t_QUAD:
2386 49 : if (!is_realquad(x)) break;
2387 42 : if (gequal0(gel(x,3))) return gceil(gel(x,2));
2388 35 : av = avma; return gerepileupto(av, addiu(gfloor(x), 1));
2389 7 : case t_RFRAC:
2390 7 : return gdeuc(gel(x,1),gel(x,2));
2391 35 : case t_VEC: case t_COL: case t_MAT:
2392 105 : pari_APPLY_same(gceil(gel(x,i)));
2393 : }
2394 42 : pari_err_TYPE("gceil",x);
2395 : return NULL; /* LCOV_EXCL_LINE */
2396 : }
2397 :
2398 : GEN
2399 5439 : round0(GEN x, GEN *pte)
2400 : {
2401 5439 : if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
2402 5432 : return ground(x);
2403 : }
2404 :
2405 : /* x t_REAL, return q=floor(x+1/2), set e = expo(x-q) */
2406 : static GEN
2407 43984229 : round_i(GEN x, long *pe)
2408 : {
2409 : long e;
2410 43984229 : GEN B, q,r, m = mantissa_real(x, &e); /* x = m/2^e */
2411 43981380 : if (e <= 0)
2412 : {
2413 3153314 : if (e) m = shifti(m,-e);
2414 3153278 : if (pe) *pe = -e;
2415 3153278 : return m;
2416 : }
2417 40828066 : B = int2n(e-1);
2418 40828152 : m = addii(m, B);
2419 40828156 : q = shifti(m, -e);
2420 40827796 : r = remi2n(m, e);
2421 : /* 2^e (x+1/2) = m = 2^e q + r, sgn(r)=sgn(m), |r|<2^e */
2422 40829758 : if (!signe(r))
2423 8195 : { if (pe) *pe = -1; }
2424 : else
2425 : {
2426 40821563 : if (signe(m) < 0)
2427 : {
2428 19652114 : q = subiu(q,1);
2429 19652048 : r = addii(r, B);
2430 : }
2431 : else
2432 21169449 : r = subii(r, B);
2433 : /* |x - q| = |r| / 2^e */
2434 40821305 : if (pe) *pe = signe(r)? expi(r) - e: -e;
2435 40821798 : cgiv(r);
2436 : }
2437 40830383 : return q;
2438 : }
2439 : /* assume x a t_REAL */
2440 : GEN
2441 3082044 : roundr(GEN x)
2442 : {
2443 3082044 : long ex, s = signe(x);
2444 : pari_sp av;
2445 3082044 : if (!s || (ex=expo(x)) < -1) return gen_0;
2446 2381008 : if (ex == -1) return s>0? gen_1:
2447 211685 : absrnz_equal2n(x)? gen_0: gen_m1;
2448 1740610 : av = avma; x = round_i(x, &ex);
2449 1740674 : if (ex >= 0) pari_err_PREC( "roundr (precision loss in truncation)");
2450 1740674 : return gerepileuptoint(av, x);
2451 : }
2452 : GEN
2453 332013 : roundr_safe(GEN x)
2454 : {
2455 332013 : long ex, s = signe(x);
2456 : pari_sp av;
2457 :
2458 332013 : if (!s || (ex = expo(x)) < -1) return gen_0;
2459 331969 : if (ex == -1) return s>0? gen_1:
2460 0 : absrnz_equal2n(x)? gen_0: gen_m1;
2461 331942 : av = avma; x = round_i(x, NULL);
2462 331942 : return gerepileuptoint(av, x);
2463 : }
2464 :
2465 : GEN
2466 2622316 : ground(GEN x)
2467 : {
2468 : pari_sp av;
2469 : GEN y;
2470 :
2471 2622316 : switch(typ(x))
2472 : {
2473 549860 : case t_INT: return icopy(x);
2474 14 : case t_INTMOD: return gcopy(x);
2475 1521566 : case t_REAL: return roundr(x);
2476 56637 : case t_FRAC: return diviiround(gel(x,1), gel(x,2));
2477 49 : case t_QUAD:
2478 : {
2479 49 : GEN Q = gel(x,1), u, v, b, d, z;
2480 49 : av = avma;
2481 49 : if (is_realquad(x)) return gerepileupto(av, gfloor(gadd(x, ghalf)));
2482 7 : u = gel(x,2);
2483 7 : v = gel(x,3); b = gel(Q,3);
2484 7 : u = ground(gsub(u, gmul2n(gmul(v,b),-1)));
2485 7 : v = Q_remove_denom(v, &d);
2486 7 : if (!d) d = gen_1;
2487 : /* Im x = v sqrt(|D|) / (2d),
2488 : * Im(round(x)) = floor((d + v sqrt(|D|)) / (2d))
2489 : * = floor(floor(d + v sqrt(|D|)) / (2d)) */
2490 7 : z = sqrti(mulii(sqri(v), quad_disc(x)));
2491 7 : if (signe(v) < 0) { z = addiu(z,1); togglesign(z); }
2492 : /* z = floor(v * sqrt(|D|)) */
2493 7 : v = truedivii(addii(z, d), shifti(d,1));
2494 7 : return gerepilecopy(av, signe(v)? mkcomplex(u,v): u);
2495 : }
2496 14 : case t_POLMOD:
2497 14 : retmkpolmod(ground(gel(x,2)), RgX_copy(gel(x,1)));
2498 :
2499 4333 : case t_COMPLEX:
2500 4333 : av = avma; y = cgetg(3, t_COMPLEX);
2501 4333 : gel(y,2) = ground(gel(x,2));
2502 4333 : if (!signe(gel(y,2))) { set_avma(av); return ground(gel(x,1)); }
2503 490 : gel(y,1) = ground(gel(x,1)); return y;
2504 :
2505 91 : case t_POL:
2506 4011 : pari_APPLY_pol(ground(gel(x,i)));
2507 182 : case t_SER:
2508 182 : if (ser_isexactzero(x)) return gcopy(x);
2509 42 : pari_APPLY_ser(ground(gel(x,i)));
2510 21 : case t_RFRAC:
2511 21 : av = avma;
2512 21 : return gerepileupto(av, gdiv(ground(gel(x,1)), ground(gel(x,2))));
2513 489578 : case t_VEC: case t_COL: case t_MAT:
2514 2355117 : pari_APPLY_same(ground(gel(x,i)));
2515 : }
2516 6 : pari_err_TYPE("ground",x);
2517 : return NULL; /* LCOV_EXCL_LINE */
2518 : }
2519 :
2520 : /* e = number of error bits on integral part */
2521 : GEN
2522 50871010 : grndtoi(GEN x, long *e)
2523 : {
2524 : GEN y;
2525 : long i, lx, e1;
2526 : pari_sp av;
2527 :
2528 50871010 : if (e) *e = -(long)HIGHEXPOBIT;
2529 50871010 : switch(typ(x))
2530 : {
2531 2165776 : case t_INT: return icopy(x);
2532 43574181 : case t_REAL: {
2533 43574181 : long ex = expo(x);
2534 43574181 : if (!signe(x) || ex < -1)
2535 : {
2536 1662420 : if (e) *e = ex;
2537 1662420 : return gen_0;
2538 : }
2539 41911761 : av = avma; x = round_i(x, e);
2540 41910835 : return gerepileuptoint(av, x);
2541 : }
2542 29084 : case t_FRAC:
2543 29084 : y = diviiround(gel(x,1), gel(x,2)); if (!e) return y;
2544 28090 : av = avma; *e = gexpo(gsub(x, y)); set_avma(av);
2545 28090 : return y;
2546 7 : case t_INTMOD: return gcopy(x);
2547 7 : case t_QUAD:
2548 7 : y = ground(x); av = avma; if (!e) return y;
2549 7 : *e = gexpo(gsub(x,y)); set_avma(avma); return y;
2550 1632587 : case t_COMPLEX:
2551 1632587 : av = avma; y = cgetg(3, t_COMPLEX);
2552 1632587 : gel(y,2) = grndtoi(gel(x,2), e);
2553 1632587 : if (!signe(gel(y,2))) {
2554 243977 : set_avma(av);
2555 243977 : y = grndtoi(gel(x,1), e? &e1: NULL);
2556 : }
2557 : else
2558 1388610 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL);
2559 1632587 : if (e && e1 > *e) *e = e1;
2560 1632587 : return y;
2561 :
2562 7 : case t_POLMOD:
2563 7 : retmkpolmod(grndtoi(gel(x,2), e), RgX_copy(gel(x,1)));
2564 :
2565 314459 : case t_POL:
2566 314459 : y = cgetg_copy(x, &lx); y[1] = x[1];
2567 3307492 : for (i=2; i<lx; i++)
2568 : {
2569 2993049 : gel(y,i) = grndtoi(gel(x,i), &e1);
2570 2993033 : if (e && e1 > *e) *e = e1;
2571 : }
2572 314443 : return normalizepol_lg(y, lx);
2573 168 : case t_SER:
2574 168 : if (ser_isexactzero(x)) return gcopy(x);
2575 154 : y = cgetg_copy(x, &lx); y[1] = x[1];
2576 462 : for (i=2; i<lx; i++)
2577 : {
2578 308 : gel(y,i) = grndtoi(gel(x,i), &e1);
2579 308 : if (e && e1 > *e) *e = e1;
2580 : }
2581 154 : return normalizeser(y);
2582 7 : case t_RFRAC:
2583 7 : y = cgetg(3,t_RFRAC);
2584 7 : gel(y,1) = grndtoi(gel(x,1), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2585 7 : gel(y,2) = grndtoi(gel(x,2), e? &e1: NULL); if (e && e1 > *e) *e = e1;
2586 7 : return y;
2587 3155060 : case t_VEC: case t_COL: case t_MAT:
2588 3155060 : y = cgetg_copy(x, &lx);
2589 13675384 : for (i=1; i<lx; i++)
2590 : {
2591 10520465 : gel(y,i) = grndtoi(gel(x,i), e? &e1: NULL);
2592 10520322 : if (e && e1 > *e) *e = e1;
2593 : }
2594 3154919 : return y;
2595 : }
2596 6 : pari_err_TYPE("grndtoi",x);
2597 : return NULL; /* LCOV_EXCL_LINE */
2598 : }
2599 :
2600 : /* trunc(x * 2^s). lindep() sanity checks rely on this function to return a
2601 : * t_INT or fail when fed a non-t_COMPLEX input; so do not make this one
2602 : * recursive [ or change the lindep call ] */
2603 : GEN
2604 90290827 : gtrunc2n(GEN x, long s)
2605 : {
2606 : GEN z;
2607 90290827 : switch(typ(x))
2608 : {
2609 27768157 : case t_INT: return shifti(x, s);
2610 46206383 : case t_REAL: return trunc2nr(x, s);
2611 483 : case t_FRAC: {
2612 : pari_sp av;
2613 483 : GEN a = gel(x,1), b = gel(x,2), q;
2614 483 : if (s == 0) return divii(a, b);
2615 483 : av = avma;
2616 483 : if (s < 0) q = divii(shifti(a, s), b);
2617 : else {
2618 : GEN r;
2619 483 : q = dvmdii(a, b, &r);
2620 483 : q = addii(shifti(q,s), divii(shifti(r,s), b));
2621 : }
2622 483 : return gerepileuptoint(av, q);
2623 : }
2624 16403948 : case t_COMPLEX:
2625 16403948 : z = cgetg(3, t_COMPLEX);
2626 16407520 : gel(z,2) = gtrunc2n(gel(x,2), s);
2627 16402977 : if (!signe(gel(z,2))) {
2628 1550922 : set_avma((pari_sp)(z + 3));
2629 1550924 : return gtrunc2n(gel(x,1), s);
2630 : }
2631 14852055 : gel(z,1) = gtrunc2n(gel(x,1), s);
2632 14850831 : return z;
2633 0 : default: pari_err_TYPE("gtrunc2n",x);
2634 : return NULL; /* LCOV_EXCL_LINE */
2635 : }
2636 : }
2637 :
2638 : /* e = number of error bits on integral part */
2639 : GEN
2640 871850 : gcvtoi(GEN x, long *e)
2641 : {
2642 871850 : long tx = typ(x), lx, e1;
2643 : GEN y;
2644 :
2645 871850 : if (tx == t_REAL)
2646 : {
2647 871718 : long ex = expo(x); if (ex < 0) { *e = ex; return gen_0; }
2648 871612 : e1 = ex - bit_prec(x) + 1;
2649 871616 : y = mantissa2nr(x, e1);
2650 871607 : if (e1 <= 0) { pari_sp av = avma; e1 = expo(subri(x,y)); set_avma(av); }
2651 871594 : *e = e1; return y;
2652 : }
2653 132 : *e = -(long)HIGHEXPOBIT;
2654 132 : if (is_matvec_t(tx))
2655 : {
2656 : long i;
2657 28 : y = cgetg_copy(x, &lx);
2658 84 : for (i=1; i<lx; i++)
2659 : {
2660 56 : gel(y,i) = gcvtoi(gel(x,i),&e1);
2661 56 : if (e1 > *e) *e = e1;
2662 : }
2663 28 : return y;
2664 : }
2665 105 : return gtrunc(x);
2666 : }
2667 :
2668 : int
2669 438193 : isint(GEN n, GEN *ptk)
2670 : {
2671 438193 : switch(typ(n))
2672 : {
2673 390649 : case t_INT: *ptk = n; return 1;
2674 1260 : case t_REAL: {
2675 1260 : pari_sp av0 = avma;
2676 1260 : GEN z = floorr(n);
2677 1260 : pari_sp av = avma;
2678 1260 : long s = signe(subri(n, z));
2679 1260 : if (s) return gc_bool(av0,0);
2680 21 : *ptk = z; return gc_bool(av,1);
2681 : }
2682 29134 : case t_FRAC: return 0;
2683 16961 : case t_COMPLEX: return gequal0(gel(n,2)) && isint(gel(n,1),ptk);
2684 0 : case t_QUAD: return gequal0(gel(n,3)) && isint(gel(n,2),ptk);
2685 189 : default: pari_err_TYPE("isint",n);
2686 : return 0; /* LCOV_EXCL_LINE */
2687 : }
2688 : }
2689 :
2690 : int
2691 321846 : issmall(GEN n, long *ptk)
2692 : {
2693 321846 : pari_sp av = avma;
2694 : GEN z;
2695 : long k;
2696 321846 : if (!isint(n, &z)) return 0;
2697 320278 : k = itos_or_0(z); set_avma(av);
2698 320278 : if (k || lgefint(z) == 2) { *ptk = k; return 1; }
2699 0 : return 0;
2700 : }
2701 :
2702 : /* smallest integer greater than any incarnations of the real x
2703 : * Avoid "precision loss in truncation" */
2704 : GEN
2705 799492 : ceil_safe(GEN x)
2706 : {
2707 799492 : pari_sp av = avma;
2708 799492 : long e, tx = typ(x);
2709 : GEN y;
2710 :
2711 799492 : if (is_rational_t(tx)) return gceil(x);
2712 799236 : y = gcvtoi(x,&e);
2713 799222 : if (gsigne(x) >= 0)
2714 : {
2715 798723 : if (e < 0) e = 0;
2716 798723 : y = addii(y, int2n(e));
2717 : }
2718 799182 : return gerepileuptoint(av, y);
2719 : }
2720 : /* largest integer smaller than any incarnations of the real x
2721 : * Avoid "precision loss in truncation" */
2722 : GEN
2723 61504 : floor_safe(GEN x)
2724 : {
2725 61504 : pari_sp av = avma;
2726 61504 : long e, tx = typ(x);
2727 : GEN y;
2728 :
2729 61504 : if (is_rational_t(tx)) return gfloor(x);
2730 61505 : y = gcvtoi(x,&e);
2731 61504 : if (gsigne(x) <= 0)
2732 : {
2733 21 : if (e < 0) e = 0;
2734 21 : y = subii(y, int2n(e));
2735 : }
2736 61504 : return gerepileuptoint(av, y);
2737 : }
2738 :
2739 : GEN
2740 13258 : ser2rfrac_i(GEN x)
2741 : {
2742 13258 : long e = valp(x);
2743 13258 : GEN a = ser2pol_i(x, lg(x));
2744 13258 : if (e) {
2745 5222 : if (e > 0) a = RgX_shift_shallow(a, e);
2746 0 : else a = gred_rfrac_simple(a, pol_xn(-e, varn(a)));
2747 : }
2748 13258 : return a;
2749 : }
2750 :
2751 : static GEN
2752 490 : ser2rfrac(GEN x)
2753 : {
2754 490 : pari_sp av = avma;
2755 490 : return gerepilecopy(av, ser2rfrac_i(x));
2756 : }
2757 :
2758 : /* x t_PADIC, truncate to rational (t_INT/t_FRAC) */
2759 : GEN
2760 269709 : padic_to_Q(GEN x)
2761 : {
2762 269709 : GEN u = gel(x,4), p;
2763 : long v;
2764 269709 : if (!signe(u)) return gen_0;
2765 266615 : v = valp(x);
2766 266615 : if (!v) return icopy(u);
2767 167806 : p = gel(x,2);
2768 167806 : if (v>0)
2769 : {
2770 167692 : pari_sp av = avma;
2771 167692 : return gerepileuptoint(av, mulii(u, powiu(p,v)));
2772 : }
2773 114 : retmkfrac(icopy(u), powiu(p,-v));
2774 : }
2775 : GEN
2776 14 : padic_to_Q_shallow(GEN x)
2777 : {
2778 14 : GEN u = gel(x,4), p, q, q2;
2779 : long v;
2780 14 : if (!signe(u)) return gen_0;
2781 14 : q = gel(x,3); q2 = shifti(q,-1);
2782 14 : v = valp(x);
2783 14 : u = Fp_center_i(u, q, q2);
2784 14 : if (!v) return u;
2785 7 : p = gel(x,2);
2786 7 : if (v>0) return mulii(powiu(p,v), u);
2787 7 : return mkfrac(u, powiu(p,-v));
2788 : }
2789 : GEN
2790 196 : QpV_to_QV(GEN v)
2791 : {
2792 : long i, l;
2793 196 : GEN w = cgetg_copy(v, &l);
2794 931 : for (i = 1; i < l; i++)
2795 : {
2796 735 : GEN c = gel(v,i);
2797 735 : switch(typ(c))
2798 : {
2799 721 : case t_INT:
2800 721 : case t_FRAC: break;
2801 14 : case t_PADIC: c = padic_to_Q_shallow(c); break;
2802 0 : default: pari_err_TYPE("padic_to_Q", v);
2803 : }
2804 735 : gel(w,i) = c;
2805 : }
2806 196 : return w;
2807 : }
2808 :
2809 : GEN
2810 185867 : gtrunc(GEN x)
2811 : {
2812 185867 : switch(typ(x))
2813 : {
2814 82523 : case t_INT: return icopy(x);
2815 48454 : case t_REAL: return truncr(x);
2816 56 : case t_FRAC: return divii(gel(x,1),gel(x,2));
2817 23915 : case t_PADIC: return padic_to_Q(x);
2818 42 : case t_POL: return RgX_copy(x);
2819 14 : case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
2820 462 : case t_SER: return ser2rfrac(x);
2821 183925 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(gtrunc(gel(x,i)));
2822 : }
2823 56 : pari_err_TYPE("gtrunc",x);
2824 : return NULL; /* LCOV_EXCL_LINE */
2825 : }
2826 :
2827 : GEN
2828 217 : trunc0(GEN x, GEN *pte)
2829 : {
2830 217 : if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
2831 189 : return gtrunc(x);
2832 : }
2833 : /*******************************************************************/
2834 : /* */
2835 : /* CONVERSIONS --> INT, POL & SER */
2836 : /* */
2837 : /*******************************************************************/
2838 :
2839 : /* return a_(n-1) B^(n-1) + ... + a_0, where B = 2^32.
2840 : * The a_i are 32bits integers */
2841 : GEN
2842 15470 : mkintn(long n, ...)
2843 : {
2844 : va_list ap;
2845 : GEN x, y;
2846 : long i;
2847 : #ifdef LONG_IS_64BIT
2848 13260 : long e = (n&1);
2849 13260 : n = (n+1) >> 1;
2850 : #endif
2851 15470 : va_start(ap,n);
2852 15470 : x = cgetipos(n+2);
2853 15470 : y = int_MSW(x);
2854 54740 : for (i=0; i <n; i++)
2855 : {
2856 : #ifdef LONG_IS_64BIT
2857 30600 : ulong a = (e && !i)? 0: (ulong) va_arg(ap, unsigned int);
2858 30600 : ulong b = (ulong) va_arg(ap, unsigned int);
2859 30600 : *y = (a << 32) | b;
2860 : #else
2861 8670 : *y = (ulong) va_arg(ap, unsigned int);
2862 : #endif
2863 39270 : y = int_precW(y);
2864 : }
2865 15470 : va_end(ap);
2866 15470 : return int_normalize(x, 0);
2867 : }
2868 :
2869 : /* 2^32 a + b */
2870 : GEN
2871 227654 : uu32toi(ulong a, ulong b)
2872 : {
2873 : #ifdef LONG_IS_64BIT
2874 183197 : return utoi((a<<32) | b);
2875 : #else
2876 44457 : return uutoi(a, b);
2877 : #endif
2878 : }
2879 : /* - (2^32 a + b), assume a or b != 0 */
2880 : GEN
2881 72276 : uu32toineg(ulong a, ulong b)
2882 : {
2883 : #ifdef LONG_IS_64BIT
2884 61909 : return utoineg((a<<32) | b);
2885 : #else
2886 10367 : return uutoineg(a, b);
2887 : #endif
2888 : }
2889 :
2890 : /* return a_(n-1) x^(n-1) + ... + a_0 */
2891 : GEN
2892 3468297 : mkpoln(long n, ...)
2893 : {
2894 : va_list ap;
2895 : GEN x, y;
2896 3468297 : long i, l = n+2;
2897 3468297 : va_start(ap,n);
2898 3468297 : x = cgetg(l, t_POL); y = x + 2;
2899 3473820 : x[1] = evalvarn(0);
2900 15847587 : for (i=n-1; i >= 0; i--) gel(y,i) = va_arg(ap, GEN);
2901 3475133 : va_end(ap); return normalizepol_lg(x, l);
2902 : }
2903 :
2904 : /* return [a_1, ..., a_n] */
2905 : GEN
2906 1067392 : mkvecn(long n, ...)
2907 : {
2908 : va_list ap;
2909 : GEN x;
2910 : long i;
2911 1067392 : va_start(ap,n);
2912 1067392 : x = cgetg(n+1, t_VEC);
2913 6785128 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2914 1067397 : va_end(ap); return x;
2915 : }
2916 :
2917 : GEN
2918 1379 : mkcoln(long n, ...)
2919 : {
2920 : va_list ap;
2921 : GEN x;
2922 : long i;
2923 1379 : va_start(ap,n);
2924 1379 : x = cgetg(n+1, t_COL);
2925 11032 : for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
2926 1379 : va_end(ap); return x;
2927 : }
2928 :
2929 : GEN
2930 41555 : mkvecsmalln(long n, ...)
2931 : {
2932 : va_list ap;
2933 : GEN x;
2934 : long i;
2935 41555 : va_start(ap,n);
2936 41555 : x = cgetg(n+1, t_VECSMALL);
2937 268779 : for (i=1; i <= n; i++) x[i] = va_arg(ap, long);
2938 41555 : va_end(ap); return x;
2939 : }
2940 :
2941 : GEN
2942 16343928 : scalarpol(GEN x, long v)
2943 : {
2944 : GEN y;
2945 16343928 : if (isrationalzero(x)) return zeropol(v);
2946 4593611 : y = cgetg(3,t_POL);
2947 4593668 : y[1] = gequal0(x)? evalvarn(v)
2948 4593667 : : evalvarn(v) | evalsigne(1);
2949 4593667 : gel(y,2) = gcopy(x); return y;
2950 : }
2951 : GEN
2952 2779468 : scalarpol_shallow(GEN x, long v)
2953 : {
2954 : GEN y;
2955 2779468 : if (isrationalzero(x)) return zeropol(v);
2956 2745040 : y = cgetg(3,t_POL);
2957 2745044 : y[1] = gequal0(x)? evalvarn(v)
2958 2745044 : : evalvarn(v) | evalsigne(1);
2959 2745044 : gel(y,2) = x; return y;
2960 : }
2961 :
2962 : /* x0 + x1*T, do not assume x1 != 0 */
2963 : GEN
2964 437978 : deg1pol(GEN x1, GEN x0,long v)
2965 : {
2966 437978 : GEN x = cgetg(4,t_POL);
2967 437979 : x[1] = evalsigne(1) | evalvarn(v);
2968 437979 : gel(x,2) = x0 == gen_0? x0: gcopy(x0); /* gen_0 frequent */
2969 437980 : gel(x,3) = gcopy(x1); return normalizepol_lg(x,4);
2970 : }
2971 :
2972 : /* same, no copy */
2973 : GEN
2974 7229184 : deg1pol_shallow(GEN x1, GEN x0,long v)
2975 : {
2976 7229184 : GEN x = cgetg(4,t_POL);
2977 7229361 : x[1] = evalsigne(1) | evalvarn(v);
2978 7229361 : gel(x,2) = x0;
2979 7229361 : gel(x,3) = x1; return normalizepol_lg(x,4);
2980 : }
2981 :
2982 : GEN
2983 309395 : deg2pol_shallow(GEN x2, GEN x1, GEN x0, long v)
2984 : {
2985 309395 : GEN x = cgetg(5,t_POL);
2986 309395 : x[1] = evalsigne(1) | evalvarn(v);
2987 309395 : gel(x,2) = x0;
2988 309395 : gel(x,3) = x1;
2989 309395 : gel(x,4) = x2;
2990 309395 : return normalizepol_lg(x,5);
2991 : }
2992 :
2993 : static GEN
2994 3472383 : _gtopoly(GEN x, long v, int reverse)
2995 : {
2996 3472383 : long tx = typ(x);
2997 : GEN y;
2998 :
2999 3472383 : if (v<0) v = 0;
3000 3472383 : switch(tx)
3001 : {
3002 21 : case t_POL:
3003 21 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3004 21 : y = RgX_copy(x); break;
3005 28 : case t_SER:
3006 28 : if (varncmp(varn(x), v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3007 28 : y = ser2rfrac(x);
3008 28 : if (typ(y) != t_POL)
3009 0 : pari_err_DOMAIN("gtopoly", "valuation", "<", gen_0, x);
3010 28 : break;
3011 42 : case t_RFRAC:
3012 : {
3013 42 : GEN a = gel(x,1), b = gel(x,2);
3014 42 : long vb = varn(b);
3015 42 : if (varncmp(vb, v) < 0) pari_err_PRIORITY("gtopoly", x, "<", v);
3016 42 : if (typ(a) != t_POL || varn(a) != vb) return zeropol(v);
3017 21 : y = RgX_div(a,b); break;
3018 : }
3019 337043 : case t_VECSMALL: x = zv_to_ZV(x); /* fall through */
3020 3471746 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
3021 : {
3022 3471746 : long j, k, lx = lg(x);
3023 : GEN z;
3024 3471746 : if (tx == t_QFB) lx--;
3025 3471746 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("gtopoly", x, "<=", v);
3026 3471746 : y = cgetg(lx+1, t_POL);
3027 3471742 : y[1] = evalvarn(v);
3028 3471742 : if (reverse) {
3029 2433795 : x--;
3030 26223351 : for (j=2; j<=lx; j++) gel(y,j) = gel(x,j);
3031 : } else {
3032 5311218 : for (j=2, k=lx; k>=2; j++) gel(y,j) = gel(x,--k);
3033 : }
3034 3471742 : z = RgX_copy(normalizepol_lg(y,lx+1));
3035 3471751 : settyp(y, t_VECSMALL);/* left on stack */
3036 3471751 : return z;
3037 : }
3038 546 : default:
3039 546 : if (is_scalar_t(tx)) return scalarpol(x,v);
3040 7 : pari_err_TYPE("gtopoly",x);
3041 : return NULL; /* LCOV_EXCL_LINE */
3042 : }
3043 70 : setvarn(y,v); return y;
3044 : }
3045 :
3046 : GEN
3047 2433830 : gtopolyrev(GEN x, long v) { return _gtopoly(x,v,1); }
3048 :
3049 : GEN
3050 1038551 : gtopoly(GEN x, long v) { return _gtopoly(x,v,0); }
3051 :
3052 : static GEN
3053 1064 : gtovecpost(GEN x, long n)
3054 : {
3055 1064 : long i, imax, lx, tx = typ(x);
3056 1064 : GEN y = zerovec(n);
3057 :
3058 1064 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,1) = gcopy(x); return y; }
3059 315 : switch(tx)
3060 : {
3061 56 : case t_POL:
3062 56 : lx=lg(x); imax = minss(lx-2, n);
3063 224 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,lx-i));
3064 56 : return y;
3065 28 : case t_SER:
3066 28 : lx=lg(x); imax = minss(lx-2, n); x++;
3067 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3068 28 : return y;
3069 28 : case t_QFB: case t_VEC: case t_COL:
3070 28 : lx=lg(x); imax = minss(lx-1, n);
3071 84 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3072 28 : return y;
3073 63 : case t_LIST:
3074 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3075 56 : x = list_data(x); lx = x? lg(x): 1;
3076 56 : imax = minss(lx-1, n);
3077 252 : for (i=1; i<=imax; i++) gel(y,i) = gcopy(gel(x,i));
3078 56 : return y;
3079 56 : case t_STR:
3080 : {
3081 56 : char *s = GSTR(x);
3082 56 : imax = minss(strlen(s), n); s--;
3083 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3084 56 : return y;
3085 : }
3086 28 : case t_VECSMALL:
3087 28 : lx=lg(x);
3088 28 : imax = minss(lx-1, n);
3089 84 : for (i=1; i<=imax; i++) gel(y,i) = stoi(x[i]);
3090 28 : return y;
3091 56 : default: pari_err_TYPE("gtovec",x);
3092 : return NULL; /*LCOV_EXCL_LINE*/
3093 : }
3094 : }
3095 :
3096 : static GEN
3097 3535 : init_vectopre(long a, long n, GEN y, long *imax)
3098 : {
3099 3535 : *imax = minss(a, n);
3100 3535 : return (n == *imax)? y: y + n - a;
3101 : }
3102 : static GEN
3103 3591 : gtovecpre(GEN x, long n)
3104 : {
3105 3591 : long i, imax, lx, tx = typ(x);
3106 3591 : GEN y = zerovec(n), y0;
3107 :
3108 3591 : if (is_scalar_t(tx) || tx == t_RFRAC) { gel(y,n) = gcopy(x); return y; }
3109 3535 : switch(tx)
3110 : {
3111 56 : case t_POL:
3112 56 : lx=lg(x);
3113 56 : y0 = init_vectopre(lx-2, n, y, &imax);
3114 224 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,lx-i));
3115 56 : return y;
3116 3248 : case t_SER:
3117 3248 : lx=lg(x); x++;
3118 3248 : y0 = init_vectopre(lx-2, n, y, &imax);
3119 131425 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3120 3248 : return y;
3121 28 : case t_QFB: case t_VEC: case t_COL:
3122 28 : lx=lg(x);
3123 28 : y0 = init_vectopre(lx-1, n, y, &imax);
3124 84 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3125 28 : return y;
3126 63 : case t_LIST:
3127 63 : if (list_typ(x) == t_LIST_MAP) pari_err_TYPE("gtovec",x);
3128 56 : x = list_data(x); lx = x? lg(x): 1;
3129 56 : y0 = init_vectopre(lx-1, n, y, &imax);
3130 252 : for (i=1; i<=imax; i++) gel(y0,i) = gcopy(gel(x,i));
3131 56 : return y;
3132 56 : case t_STR:
3133 : {
3134 56 : char *s = GSTR(x);
3135 56 : y0 = init_vectopre(strlen(s), n, y, &imax); s--;
3136 224 : for (i=1; i<=imax; i++) gel(y,i) = chartoGENstr(s[i]);
3137 56 : return y;
3138 : }
3139 28 : case t_VECSMALL:
3140 28 : lx=lg(x);
3141 28 : y0 = init_vectopre(lx-1, n, y, &imax);
3142 84 : for (i=1; i<=imax; i++) gel(y0,i) = stoi(x[i]);
3143 28 : return y;
3144 56 : default: pari_err_TYPE("gtovec",x);
3145 : return NULL; /*LCOV_EXCL_LINE*/
3146 : }
3147 : }
3148 : GEN
3149 106000 : gtovec0(GEN x, long n)
3150 : {
3151 106000 : if (!n) return gtovec(x);
3152 4655 : if (n > 0) return gtovecpost(x, n);
3153 3591 : return gtovecpre(x, -n);
3154 : }
3155 :
3156 : GEN
3157 101660 : gtovec(GEN x)
3158 : {
3159 101660 : long i, lx, tx = typ(x);
3160 : GEN y;
3161 :
3162 101660 : if (is_scalar_t(tx)) return mkveccopy(x);
3163 101527 : switch(tx)
3164 : {
3165 15218 : case t_POL:
3166 15218 : lx=lg(x); y=cgetg(lx-1,t_VEC);
3167 1497342 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,lx-i));
3168 15218 : return y;
3169 385 : case t_SER:
3170 385 : lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
3171 12264 : for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
3172 385 : return y;
3173 28 : case t_RFRAC: return mkveccopy(x);
3174 70007 : case t_QFB:
3175 70007 : retmkvec3(icopy(gel(x,1)), icopy(gel(x,2)), icopy(gel(x,3)));
3176 14168 : case t_VEC: case t_COL: case t_MAT:
3177 14168 : lx=lg(x); y=cgetg(lx,t_VEC);
3178 1467200 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3179 14168 : return y;
3180 76 : case t_LIST:
3181 76 : if (list_typ(x) == t_LIST_MAP) return mapdomain(x);
3182 62 : x = list_data(x); lx = x? lg(x): 1;
3183 62 : y = cgetg(lx, t_VEC);
3184 304 : for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3185 62 : return y;
3186 105 : case t_STR:
3187 : {
3188 105 : char *s = GSTR(x);
3189 105 : lx = strlen(s)+1; y = cgetg(lx, t_VEC);
3190 105 : s--;
3191 340312 : for (i=1; i<lx; i++) gel(y,i) = chartoGENstr(s[i]);
3192 105 : return y;
3193 : }
3194 1498 : case t_VECSMALL:
3195 1498 : return vecsmall_to_vec(x);
3196 42 : case t_ERROR:
3197 42 : lx=lg(x); y = cgetg(lx,t_VEC);
3198 42 : gel(y,1) = errname(x);
3199 126 : for (i=2; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
3200 42 : return y;
3201 0 : default: pari_err_TYPE("gtovec",x);
3202 : return NULL; /*LCOV_EXCL_LINE*/
3203 : }
3204 : }
3205 :
3206 : GEN
3207 287 : gtovecrev0(GEN x, long n)
3208 : {
3209 287 : GEN y = gtovec0(x, -n);
3210 259 : vecreverse_inplace(y);
3211 259 : return y;
3212 : }
3213 : GEN
3214 0 : gtovecrev(GEN x) { return gtovecrev0(x, 0); }
3215 :
3216 : GEN
3217 3626 : gtocol0(GEN x, long n)
3218 : {
3219 : GEN y;
3220 3626 : if (!n) return gtocol(x);
3221 3430 : y = gtovec0(x, n);
3222 3374 : settyp(y, t_COL); return y;
3223 : }
3224 : GEN
3225 196 : gtocol(GEN x)
3226 : {
3227 : long lx, tx, i, j, h;
3228 : GEN y;
3229 196 : tx = typ(x);
3230 196 : if (tx != t_MAT) { y = gtovec(x); settyp(y, t_COL); return y; }
3231 14 : lx = lg(x); if (lx == 1) return cgetg(1, t_COL);
3232 14 : h = lgcols(x); y = cgetg(h, t_COL);
3233 42 : for (i = 1 ; i < h; i++) {
3234 28 : gel(y,i) = cgetg(lx, t_VEC);
3235 112 : for (j = 1; j < lx; j++) gmael(y,i,j) = gcopy(gcoeff(x,i,j));
3236 : }
3237 14 : return y;
3238 : }
3239 :
3240 : GEN
3241 273 : gtocolrev0(GEN x, long n)
3242 : {
3243 273 : GEN y = gtocol0(x, -n);
3244 245 : long ly = lg(y), lim = ly>>1, i;
3245 700 : for (i = 1; i <= lim; i++) swap(gel(y,i), gel(y,ly-i));
3246 245 : return y;
3247 : }
3248 : GEN
3249 0 : gtocolrev(GEN x) { return gtocolrev0(x, 0); }
3250 :
3251 : static long
3252 18669 : Itos(GEN x)
3253 : {
3254 18669 : if (typ(x) != t_INT) pari_err_TYPE("vectosmall",x);
3255 18669 : return itos(x);
3256 : }
3257 :
3258 : static GEN
3259 91 : gtovecsmallpost(GEN x, long n)
3260 : {
3261 : long i, imax, lx;
3262 91 : GEN y = zero_Flv(n);
3263 :
3264 91 : switch(typ(x))
3265 : {
3266 7 : case t_INT:
3267 7 : y[1] = itos(x); return y;
3268 14 : case t_POL:
3269 14 : lx=lg(x); imax = minss(lx-2, n);
3270 56 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,lx-i));
3271 14 : return y;
3272 7 : case t_SER:
3273 7 : lx=lg(x); imax = minss(lx-2, n); x++;
3274 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3275 7 : return y;
3276 7 : case t_VEC: case t_COL:
3277 7 : lx=lg(x); imax = minss(lx-1, n);
3278 21 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3279 7 : return y;
3280 14 : case t_LIST:
3281 14 : x = list_data(x); lx = x? lg(x): 1;
3282 14 : imax = minss(lx-1, n);
3283 63 : for (i=1; i<=imax; i++) y[i] = Itos(gel(x,i));
3284 14 : return y;
3285 7 : case t_VECSMALL:
3286 7 : lx=lg(x);
3287 7 : imax = minss(lx-1, n);
3288 21 : for (i=1; i<=imax; i++) y[i] = x[i];
3289 7 : return y;
3290 14 : case t_STR:
3291 : {
3292 14 : unsigned char *s = (unsigned char*)GSTR(x);
3293 14 : imax = minss(strlen((const char *)s), n); s--;
3294 56 : for (i=1; i<=imax; i++) y[i] = (long)s[i];
3295 14 : return y;
3296 : }
3297 21 : default: pari_err_TYPE("gtovecsmall",x);
3298 : return NULL; /*LCOV_EXCL_LINE*/
3299 : }
3300 : }
3301 : static GEN
3302 91 : gtovecsmallpre(GEN x, long n)
3303 : {
3304 : long i, imax, lx;
3305 91 : GEN y = zero_Flv(n), y0;
3306 :
3307 91 : switch(typ(x))
3308 : {
3309 7 : case t_INT:
3310 7 : y[n] = itos(x); return y;
3311 14 : case t_POL:
3312 14 : lx=lg(x);
3313 14 : y0 = init_vectopre(lx-2, n, y, &imax);
3314 56 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,lx-i));
3315 14 : return y;
3316 7 : case t_SER:
3317 7 : lx=lg(x); x++;
3318 7 : y0 = init_vectopre(lx-2, n, y, &imax);
3319 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3320 7 : return y;
3321 7 : case t_VEC: case t_COL:
3322 7 : lx=lg(x);
3323 7 : y0 = init_vectopre(lx-1, n, y, &imax);
3324 21 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3325 7 : return y;
3326 14 : case t_LIST:
3327 14 : x = list_data(x); lx = x? lg(x): 1;
3328 14 : y0 = init_vectopre(lx-1, n, y, &imax);
3329 63 : for (i=1; i<=imax; i++) y0[i] = Itos(gel(x,i));
3330 14 : return y;
3331 7 : case t_VECSMALL:
3332 7 : lx=lg(x);
3333 7 : y0 = init_vectopre(lx-1, n, y, &imax);
3334 21 : for (i=1; i<=imax; i++) y0[i] = x[i];
3335 7 : return y;
3336 14 : case t_STR:
3337 : {
3338 14 : unsigned char *s = (unsigned char*)GSTR(x);
3339 14 : y0 = init_vectopre(strlen((const char *)s), n, y, &imax); s--;
3340 56 : for (i=1; i<=imax; i++) y0[i] = (long)s[i];
3341 14 : return y;
3342 : }
3343 21 : default: pari_err_TYPE("gtovecsmall",x);
3344 : return NULL; /*LCOV_EXCL_LINE*/
3345 : }
3346 : }
3347 :
3348 : GEN
3349 7686 : gtovecsmall0(GEN x, long n)
3350 : {
3351 7686 : if (!n) return gtovecsmall(x);
3352 182 : if (n > 0) return gtovecsmallpost(x, n);
3353 91 : return gtovecsmallpre(x, -n);
3354 : }
3355 :
3356 : GEN
3357 18669 : gtovecsmall(GEN x)
3358 : {
3359 : GEN V;
3360 : long l, i;
3361 :
3362 18669 : switch(typ(x))
3363 : {
3364 112 : case t_INT: return mkvecsmall(itos(x));
3365 28 : case t_STR:
3366 : {
3367 28 : unsigned char *s = (unsigned char*)GSTR(x);
3368 28 : l = strlen((const char *)s);
3369 28 : V = cgetg(l+1, t_VECSMALL);
3370 28 : s--;
3371 1953 : for (i=1; i<=l; i++) V[i] = (long)s[i];
3372 28 : return V;
3373 : }
3374 13083 : case t_VECSMALL: return leafcopy(x);
3375 14 : case t_LIST:
3376 14 : x = list_data(x);
3377 14 : if (!x) return cgetg(1, t_VECSMALL);
3378 : /* fall through */
3379 : case t_VEC: case t_COL:
3380 5404 : l = lg(x); V = cgetg(l,t_VECSMALL);
3381 23772 : for(i=1; i<l; i++) V[i] = Itos(gel(x,i));
3382 5404 : return V;
3383 14 : case t_POL:
3384 14 : l = lg(x); V = cgetg(l-1,t_VECSMALL);
3385 63 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,l-i));
3386 14 : return V;
3387 7 : case t_SER:
3388 7 : l = lg(x); V = cgetg(l-1,t_VECSMALL); x++;
3389 21 : for (i=1; i<=l-2; i++) V[i] = Itos(gel(x,i));
3390 7 : return V;
3391 21 : default:
3392 21 : pari_err_TYPE("vectosmall",x);
3393 : return NULL; /* LCOV_EXCL_LINE */
3394 : }
3395 : }
3396 :
3397 : GEN
3398 313 : compo(GEN x, long n)
3399 : {
3400 313 : long tx = typ(x);
3401 313 : ulong l, lx = (ulong)lg(x);
3402 :
3403 313 : if (!is_recursive_t(tx))
3404 : {
3405 63 : if (tx == t_VECSMALL)
3406 : {
3407 21 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3408 21 : if ((ulong)n >= lx) pari_err_COMPONENT("", ">", utoi(lx-1), stoi(n));
3409 7 : return stoi(x[n]);
3410 : }
3411 42 : pari_err_TYPE("component [leaf]", x);
3412 : }
3413 250 : if (n < 1) pari_err_COMPONENT("", "<", gen_1, stoi(n));
3414 243 : if (tx == t_LIST) {
3415 28 : x = list_data(x); tx = t_VEC;
3416 28 : lx = (ulong)(x? lg(x): 1);
3417 : }
3418 243 : l = (ulong)lontyp[tx] + (ulong)n-1; /* beware overflow */
3419 243 : if (l >= lx) pari_err_COMPONENT("", ">", utoi(lx-lontyp[tx]), stoi(n));
3420 166 : return gcopy(gel(x,l));
3421 : }
3422 :
3423 : /* assume x a t_POL */
3424 : static GEN
3425 2124228 : _polcoef(GEN x, long n, long v)
3426 : {
3427 2124228 : long i, w, lx = lg(x), dx = lx-3;
3428 : GEN z;
3429 2124228 : if (dx < 0) return gen_0;
3430 1483014 : if (v < 0 || v == (w=varn(x)))
3431 1113152 : return (n < 0 || n > dx)? gen_0: gel(x,n+2);
3432 369862 : if (varncmp(w,v) > 0) return n? gen_0: x;
3433 : /* w < v */
3434 369022 : z = cgetg(lx, t_POL); z[1] = x[1];
3435 1779911 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3436 369022 : z = normalizepol_lg(z, lx);
3437 369020 : switch(lg(z))
3438 : {
3439 27805 : case 2: z = gen_0; break;
3440 123845 : case 3: z = gel(z,2); break;
3441 : }
3442 369020 : return z;
3443 : }
3444 :
3445 : /* assume x a t_SER */
3446 : static GEN
3447 110901 : _sercoef(GEN x, long n, long v)
3448 : {
3449 110901 : long i, w = varn(x), lx = lg(x), dx = lx-3, N;
3450 : GEN z;
3451 110901 : if (v < 0) v = w;
3452 110901 : N = v == w? n - valp(x): n;
3453 110901 : if (dx < 0)
3454 : {
3455 21 : if (N >= 0) pari_err_DOMAIN("polcoef", "t_SER", "=", x, x);
3456 14 : return gen_0;
3457 : }
3458 110880 : if (v == w)
3459 : {
3460 110838 : if (!dx && !signe(x) && !isinexact(gel(x,2))) dx = -1;
3461 110838 : if (N > dx)
3462 28 : pari_err_DOMAIN("polcoef", "degree", ">", stoi(dx+valp(x)), stoi(n));
3463 110810 : return (N < 0)? gen_0: gel(x,N+2);
3464 : }
3465 42 : if (varncmp(w,v) > 0) return N? gen_0: x;
3466 : /* w < v */
3467 28 : z = cgetg(lx, t_SER); z[1] = x[1];
3468 91 : for (i = 2; i < lx; i++) gel(z,i) = polcoef_i(gel(x,i), n, v);
3469 28 : return normalizeser(z);
3470 : }
3471 :
3472 : /* assume x a t_RFRAC(n) */
3473 : static GEN
3474 21 : _rfraccoef(GEN x, long n, long v)
3475 : {
3476 21 : GEN P,Q, p = gel(x,1), q = gel(x,2);
3477 21 : long vp = gvar(p), vq = gvar(q);
3478 21 : if (v < 0) v = varncmp(vp, vq) < 0? vp: vq;
3479 21 : P = (vp == v)? p: swap_vars(p, v);
3480 21 : Q = (vq == v)? q: swap_vars(q, v);
3481 21 : if (!RgX_is_monomial(Q)) pari_err_TYPE("polcoef", x);
3482 21 : n += degpol(Q);
3483 21 : return gdiv(_polcoef(P, n, v), leading_coeff(Q));
3484 : }
3485 :
3486 : GEN
3487 2817602 : polcoef_i(GEN x, long n, long v)
3488 : {
3489 2817602 : long tx = typ(x);
3490 2817602 : switch(tx)
3491 : {
3492 2124208 : case t_POL: return _polcoef(x,n,v);
3493 110901 : case t_SER: return _sercoef(x,n,v);
3494 21 : case t_RFRAC: return _rfraccoef(x,n,v);
3495 : }
3496 582472 : if (!is_scalar_t(tx)) pari_err_TYPE("polcoef", x);
3497 582306 : return n? gen_0: x;
3498 : }
3499 :
3500 : /* with respect to the main variable if v<0, with respect to the variable v
3501 : * otherwise. v ignored if x is not a polynomial/series. */
3502 : GEN
3503 727531 : polcoef(GEN x, long n, long v)
3504 : {
3505 727531 : pari_sp av = avma;
3506 727531 : x = polcoef_i(x,n,v);
3507 727321 : if (x == gen_0) return x;
3508 130697 : return (avma == av)? gcopy(x): gerepilecopy(av, x);
3509 : }
3510 :
3511 : static GEN
3512 111300 : vecdenom(GEN v, long imin, long imax)
3513 : {
3514 111300 : long i = imin;
3515 : GEN s;
3516 111300 : if (imin > imax) return gen_1;
3517 111300 : s = denom_i(gel(v,i));
3518 222257 : for (i++; i<=imax; i++)
3519 : {
3520 110957 : GEN t = denom_i(gel(v,i));
3521 110957 : if (t != gen_1) s = glcm(s,t);
3522 : }
3523 111300 : return s;
3524 : }
3525 : static GEN denompol(GEN x, long v);
3526 : static GEN
3527 14 : vecdenompol(GEN v, long imin, long imax, long vx)
3528 : {
3529 14 : long i = imin;
3530 : GEN s;
3531 14 : if (imin > imax) return gen_1;
3532 14 : s = denompol(gel(v,i), vx);
3533 14 : for (i++; i<=imax; i++)
3534 : {
3535 0 : GEN t = denompol(gel(v,i), vx);
3536 0 : if (t != gen_1) s = glcm(s,t);
3537 : }
3538 14 : return s;
3539 : }
3540 : GEN
3541 10357650 : denom_i(GEN x)
3542 : {
3543 10357650 : switch(typ(x))
3544 : {
3545 2901034 : case t_INT:
3546 : case t_REAL:
3547 : case t_INTMOD:
3548 : case t_FFELT:
3549 : case t_PADIC:
3550 : case t_SER:
3551 2901034 : case t_VECSMALL: return gen_1;
3552 80886 : case t_FRAC: return gel(x,2);
3553 665 : case t_COMPLEX: return vecdenom(x,1,2);
3554 69069 : case t_QUAD: return vecdenom(x,2,3);
3555 42 : case t_POLMOD: return denom_i(gel(x,2));
3556 7263408 : case t_RFRAC: return gel(x,2);
3557 966 : case t_POL: return pol_1(varn(x));
3558 41566 : case t_VEC: case t_COL: case t_MAT: return vecdenom(x, 1, lg(x)-1);
3559 : }
3560 14 : pari_err_TYPE("denom",x);
3561 : return NULL; /* LCOV_EXCL_LINE */
3562 : }
3563 : /* v has lower (or equal) priority as x's main variable */
3564 : static GEN
3565 112 : denompol(GEN x, long v)
3566 : {
3567 112 : long vx, tx = typ(x);
3568 112 : if (is_scalar_t(tx)) return gen_1;
3569 98 : switch(typ(x))
3570 : {
3571 14 : case t_SER:
3572 14 : if (varn(x) != v) return x;
3573 14 : vx = valp(x); return vx < 0? pol_xn(-vx, v): pol_1(v);
3574 56 : case t_RFRAC: x = gel(x,2); return varn(x) == v? x: pol_1(v);
3575 14 : case t_POL: return pol_1(v);
3576 14 : case t_VEC: case t_COL: case t_MAT: return vecdenompol(x, 1, lg(x)-1, v);
3577 : }
3578 0 : pari_err_TYPE("denom",x);
3579 : return NULL; /* LCOV_EXCL_LINE */
3580 : }
3581 : GEN
3582 224168 : denom(GEN x) { pari_sp av = avma; return gerepilecopy(av, denom_i(x)); }
3583 :
3584 : static GEN
3585 280 : denominator_v(GEN x, long v)
3586 : {
3587 280 : long v0 = gvar(x);
3588 : GEN d;
3589 280 : if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
3590 98 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3591 98 : d = denompol(x, v0);
3592 98 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3593 98 : return d;
3594 : }
3595 : GEN
3596 124215 : denominator(GEN x, GEN D)
3597 : {
3598 124215 : pari_sp av = avma;
3599 : GEN d;
3600 124215 : if (!D) return denom(x);
3601 280 : if (isint1(D))
3602 : {
3603 140 : d = Q_denom_safe(x);
3604 140 : if (!d) { set_avma(av); return gen_1; }
3605 105 : return gerepilecopy(av, d);
3606 : }
3607 140 : if (!gequalX(D)) pari_err_TYPE("denominator", D);
3608 140 : return gerepileupto(av, denominator_v(x, varn(D)));
3609 : }
3610 : GEN
3611 8890 : numerator(GEN x, GEN D)
3612 : {
3613 8890 : pari_sp av = avma;
3614 : long v;
3615 8890 : if (!D) return numer(x);
3616 280 : if (isint1(D)) return Q_remove_denom(x,NULL);
3617 140 : if (!gequalX(D)) pari_err_TYPE("numerator", D);
3618 140 : v = varn(D); /* optimization */
3619 140 : if (typ(x) == t_RFRAC && varn(gel(x,2)) == v) return gcopy(gel(x,2));
3620 140 : return gerepileupto(av, gmul(x, denominator_v(x,v)));
3621 : }
3622 : GEN
3623 130991 : content0(GEN x, GEN D)
3624 : {
3625 130991 : pari_sp av = avma;
3626 : long v, v0;
3627 : GEN d;
3628 130991 : if (!D) return content(x);
3629 280 : if (isint1(D))
3630 : {
3631 140 : d = Q_content_safe(x);
3632 140 : return d? d: gen_1;
3633 : }
3634 140 : if (!gequalX(D)) pari_err_TYPE("content", D);
3635 140 : v = varn(D);
3636 140 : v0 = gvar(x); if (v0 == NO_VARIABLE || varncmp(v0,v) > 0) return pol_1(v);
3637 49 : if (v0 != v) { v0 = fetch_var_higher(); x = gsubst(x, v, pol_x(v0)); }
3638 49 : d = content(x);
3639 : /* gsubst is needed because of content([x]) = x */
3640 49 : if (v0 != v) { d = gsubst(d, v0, pol_x(v)); (void)delete_var(); }
3641 49 : return gerepileupto(av, d);
3642 : }
3643 :
3644 : GEN
3645 8984045 : numer_i(GEN x)
3646 : {
3647 8984045 : switch(typ(x))
3648 : {
3649 1722310 : case t_INT:
3650 : case t_REAL:
3651 : case t_INTMOD:
3652 : case t_FFELT:
3653 : case t_PADIC:
3654 : case t_SER:
3655 : case t_VECSMALL:
3656 1722310 : case t_POL: return x;
3657 28 : case t_POLMOD: return mkpolmod(numer_i(gel(x,2)), gel(x,1));
3658 7261518 : case t_FRAC:
3659 7261518 : case t_RFRAC: return gel(x,1);
3660 175 : case t_COMPLEX:
3661 : case t_QUAD:
3662 : case t_VEC:
3663 : case t_COL:
3664 175 : case t_MAT: return gmul(denom_i(x),x);
3665 : }
3666 14 : pari_err_TYPE("numer",x);
3667 : return NULL; /* LCOV_EXCL_LINE */
3668 : }
3669 : GEN
3670 8610 : numer(GEN x) { pari_sp av = avma; return gerepilecopy(av, numer_i(x)); }
3671 :
3672 : /* Lift only intmods if v does not occur in x, lift with respect to main
3673 : * variable of x if v < 0, with respect to variable v otherwise */
3674 : GEN
3675 2019796 : lift0(GEN x, long v)
3676 : {
3677 : GEN y;
3678 :
3679 2019796 : switch(typ(x))
3680 : {
3681 30338 : case t_INT: return icopy(x);
3682 1886243 : case t_INTMOD: return v < 0? icopy(gel(x,2)): gcopy(x);
3683 91707 : case t_POLMOD:
3684 91707 : if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
3685 14 : y = cgetg(3, t_POLMOD);
3686 14 : gel(y,1) = lift0(gel(x,1),v);
3687 14 : gel(y,2) = lift0(gel(x,2),v); return y;
3688 665 : case t_PADIC: return v < 0? padic_to_Q(x): gcopy(x);
3689 8715 : case t_POL:
3690 41405 : pari_APPLY_pol(lift0(gel(x,i), v));
3691 56 : case t_SER:
3692 56 : if (ser_isexactzero(x))
3693 : {
3694 14 : if (lg(x) == 2) return gcopy(x);
3695 14 : y = scalarser(lift0(gel(x,2),v), varn(x), 1);
3696 14 : setvalp(y, valp(x)); return y;
3697 : }
3698 434 : pari_APPLY_ser(lift0(gel(x,i), v));
3699 1890 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3700 : case t_VEC: case t_COL: case t_MAT:
3701 41062 : pari_APPLY_same(lift0(gel(x,i), v));
3702 182 : default: return gcopy(x);
3703 : }
3704 : }
3705 : /* same as lift, shallow */
3706 : GEN
3707 575001 : lift_shallow(GEN x)
3708 : {
3709 : GEN y;
3710 575001 : switch(typ(x))
3711 : {
3712 194222 : case t_INTMOD: case t_POLMOD: return gel(x,2);
3713 476 : case t_PADIC: return padic_to_Q(x);
3714 0 : case t_SER:
3715 0 : if (ser_isexactzero(x))
3716 : {
3717 0 : if (lg(x) == 2) return x;
3718 0 : y = scalarser(lift_shallow(gel(x,2)), varn(x), 1);
3719 0 : setvalp(y, valp(x)); return y;
3720 : }
3721 0 : pari_APPLY_ser(lift_shallow(gel(x,i)));
3722 46312 : case t_POL:
3723 257530 : pari_APPLY_pol(lift_shallow(gel(x,i)));
3724 11179 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3725 : case t_VEC: case t_COL: case t_MAT:
3726 275618 : pari_APPLY_same(lift_shallow(gel(x,i)));
3727 322812 : default: return x;
3728 : }
3729 : }
3730 : GEN
3731 1868976 : lift(GEN x) { return lift0(x,-1); }
3732 :
3733 : GEN
3734 2003183 : liftall_shallow(GEN x)
3735 : {
3736 : GEN y;
3737 2003183 : switch(typ(x))
3738 : {
3739 533778 : case t_INTMOD: return gel(x,2);
3740 547505 : case t_POLMOD:
3741 547505 : return liftall_shallow(gel(x,2));
3742 581 : case t_PADIC: return padic_to_Q(x);
3743 555898 : case t_POL:
3744 1356222 : pari_APPLY_pol(liftall_shallow(gel(x,i)));
3745 7 : case t_SER:
3746 7 : if (ser_isexactzero(x))
3747 : {
3748 0 : if (lg(x) == 2) return x;
3749 0 : y = scalarser(liftall_shallow(gel(x,2)), varn(x), 1);
3750 0 : setvalp(y, valp(x)); return y;
3751 : }
3752 35 : pari_APPLY_ser(liftall_shallow(gel(x,i)));
3753 132762 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3754 : case t_VEC: case t_COL: case t_MAT:
3755 760515 : pari_APPLY_same(liftall_shallow(gel(x,i)));
3756 232652 : default: return x;
3757 : }
3758 : }
3759 : GEN
3760 26229 : liftall(GEN x)
3761 26229 : { pari_sp av = avma; return gerepilecopy(av, liftall_shallow(x)); }
3762 :
3763 : GEN
3764 546 : liftint_shallow(GEN x)
3765 : {
3766 : GEN y;
3767 546 : switch(typ(x))
3768 : {
3769 266 : case t_INTMOD: return gel(x,2);
3770 28 : case t_PADIC: return padic_to_Q(x);
3771 21 : case t_POL:
3772 70 : pari_APPLY_pol(liftint_shallow(gel(x,i)));
3773 14 : case t_SER:
3774 14 : if (ser_isexactzero(x))
3775 : {
3776 7 : if (lg(x) == 2) return x;
3777 7 : y = scalarser(liftint_shallow(gel(x,2)), varn(x), 1);
3778 7 : setvalp(y, valp(x)); return y;
3779 : }
3780 35 : pari_APPLY_ser(liftint_shallow(gel(x,i)));
3781 161 : case t_POLMOD: case t_COMPLEX: case t_QUAD: case t_RFRAC:
3782 : case t_VEC: case t_COL: case t_MAT:
3783 504 : pari_APPLY_same(liftint_shallow(gel(x,i)));
3784 56 : default: return x;
3785 : }
3786 : }
3787 : GEN
3788 119 : liftint(GEN x)
3789 119 : { pari_sp av = avma; return gerepilecopy(av, liftint_shallow(x)); }
3790 :
3791 : GEN
3792 21100693 : liftpol_shallow(GEN x)
3793 : {
3794 : GEN y;
3795 21100693 : switch(typ(x))
3796 : {
3797 2030243 : case t_POLMOD:
3798 2030243 : return liftpol_shallow(gel(x,2));
3799 2795599 : case t_POL:
3800 11623396 : pari_APPLY_pol(liftpol_shallow(gel(x,i)));
3801 7 : case t_SER:
3802 7 : if (ser_isexactzero(x))
3803 : {
3804 0 : if (lg(x) == 2) return x;
3805 0 : y = scalarser(liftpol(gel(x,2)), varn(x), 1);
3806 0 : setvalp(y, valp(x)); return y;
3807 : }
3808 35 : pari_APPLY_ser(liftpol_shallow(gel(x,i)));
3809 131768 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
3810 944909 : pari_APPLY_same(liftpol_shallow(gel(x,i)));
3811 16143076 : default: return x;
3812 : }
3813 : }
3814 : GEN
3815 5586 : liftpol(GEN x)
3816 5586 : { pari_sp av = avma; return gerepilecopy(av, liftpol_shallow(x)); }
3817 :
3818 : static GEN
3819 42518 : centerliftii(GEN x, GEN y)
3820 : {
3821 42518 : pari_sp av = avma;
3822 42518 : long i = cmpii(shifti(x,1), y);
3823 42518 : set_avma(av); return (i > 0)? subii(x,y): icopy(x);
3824 : }
3825 :
3826 : /* see lift0 */
3827 : GEN
3828 707 : centerlift0(GEN x, long v)
3829 707 : { return v < 0? centerlift(x): lift0(x,v); }
3830 : GEN
3831 60473 : centerlift(GEN x)
3832 : {
3833 : long v;
3834 : GEN y;
3835 60473 : switch(typ(x))
3836 : {
3837 784 : case t_INT: return icopy(x);
3838 784 : case t_INTMOD: return centerliftii(gel(x,2), gel(x,1));
3839 7 : case t_POLMOD: return gcopy(gel(x,2));
3840 1554 : case t_POL:
3841 9912 : pari_APPLY_pol(centerlift(gel(x,i)));
3842 7 : case t_SER:
3843 7 : if (ser_isexactzero(x)) return lift(x);
3844 35 : pari_APPLY_ser(centerlift(gel(x,i)));
3845 5551 : case t_COMPLEX: case t_QUAD: case t_RFRAC:
3846 : case t_VEC: case t_COL: case t_MAT:
3847 56210 : pari_APPLY_same(centerlift(gel(x,i)));
3848 51779 : case t_PADIC:
3849 51779 : if (!signe(gel(x,4))) return gen_0;
3850 41734 : v = valp(x);
3851 41734 : if (v>=0)
3852 : { /* here p^v is an integer */
3853 41727 : GEN z = centerliftii(gel(x,4), gel(x,3));
3854 : pari_sp av;
3855 41727 : if (!v) return z;
3856 27027 : av = avma; y = powiu(gel(x,2),v);
3857 27027 : return gerepileuptoint(av, mulii(y,z));
3858 : }
3859 7 : y = cgetg(3,t_FRAC);
3860 7 : gel(y,1) = centerliftii(gel(x,4), gel(x,3));
3861 7 : gel(y,2) = powiu(gel(x,2),-v);
3862 7 : return y;
3863 7 : default: return gcopy(x);
3864 : }
3865 : }
3866 :
3867 : /*******************************************************************/
3868 : /* */
3869 : /* REAL & IMAGINARY PARTS */
3870 : /* */
3871 : /*******************************************************************/
3872 :
3873 : static GEN
3874 122938354 : op_ReIm(GEN f(GEN), GEN x)
3875 : {
3876 122938354 : switch(typ(x))
3877 : {
3878 378811077 : case t_POL: pari_APPLY_pol(f(gel(x,i)));
3879 62190 : case t_SER: pari_APPLY_ser(f(gel(x,i)));
3880 42 : case t_RFRAC: {
3881 42 : pari_sp av = avma;
3882 42 : GEN n, d, dxb = conj_i(gel(x,2));
3883 42 : n = gmul(gel(x,1), dxb);
3884 42 : d = gmul(gel(x,2), dxb);
3885 42 : return gerepileupto(av, gdiv(f(n), d));
3886 : }
3887 32888071 : case t_VEC: case t_COL: case t_MAT: pari_APPLY_same(f(gel(x, i)));
3888 : }
3889 12 : pari_err_TYPE("greal/gimag",x);
3890 : return NULL; /* LCOV_EXCL_LINE */
3891 : }
3892 :
3893 : GEN
3894 222166154 : real_i(GEN x)
3895 : {
3896 222166154 : switch(typ(x))
3897 : {
3898 105264813 : case t_INT: case t_REAL: case t_FRAC:
3899 105264813 : return x;
3900 52167271 : case t_COMPLEX:
3901 52167271 : return gel(x,1);
3902 0 : case t_QUAD:
3903 0 : return gel(x,2);
3904 : }
3905 64734070 : return op_ReIm(real_i,x);
3906 : }
3907 : GEN
3908 191534295 : imag_i(GEN x)
3909 : {
3910 191534295 : switch(typ(x))
3911 : {
3912 102500518 : case t_INT: case t_REAL: case t_FRAC:
3913 102500518 : return gen_0;
3914 30923610 : case t_COMPLEX:
3915 30923610 : return gel(x,2);
3916 0 : case t_QUAD:
3917 0 : return gel(x,3);
3918 : }
3919 58110167 : return op_ReIm(imag_i,x);
3920 : }
3921 : GEN
3922 5537 : greal(GEN x)
3923 : {
3924 5537 : switch(typ(x))
3925 : {
3926 294 : case t_INT: case t_REAL:
3927 294 : return mpcopy(x);
3928 :
3929 7 : case t_FRAC:
3930 7 : return gcopy(x);
3931 :
3932 5040 : case t_COMPLEX:
3933 5040 : return gcopy(gel(x,1));
3934 :
3935 7 : case t_QUAD:
3936 7 : return gcopy(gel(x,2));
3937 : }
3938 189 : return op_ReIm(greal,x);
3939 : }
3940 : GEN
3941 29288 : gimag(GEN x)
3942 : {
3943 29288 : switch(typ(x))
3944 : {
3945 525 : case t_INT: case t_REAL: case t_FRAC:
3946 525 : return gen_0;
3947 :
3948 28168 : case t_COMPLEX:
3949 28168 : return gcopy(gel(x,2));
3950 :
3951 7 : case t_QUAD:
3952 7 : return gcopy(gel(x,3));
3953 : }
3954 588 : return op_ReIm(gimag,x);
3955 : }
3956 :
3957 : /* return Re(x * y), x and y scalars */
3958 : GEN
3959 10507381 : mulreal(GEN x, GEN y)
3960 : {
3961 10507381 : if (typ(x) == t_COMPLEX)
3962 : {
3963 10480347 : if (typ(y) == t_COMPLEX)
3964 : {
3965 9808064 : pari_sp av = avma;
3966 9808064 : GEN a = gmul(gel(x,1), gel(y,1));
3967 9808047 : GEN b = gmul(gel(x,2), gel(y,2));
3968 9808059 : return gerepileupto(av, gsub(a, b));
3969 : }
3970 672283 : x = gel(x,1);
3971 : }
3972 : else
3973 27034 : if (typ(y) == t_COMPLEX) y = gel(y,1);
3974 699317 : return gmul(x,y);
3975 : }
3976 : /* Compute Re(x * y), x and y matrices of compatible dimensions
3977 : * assume scalar entries */
3978 : GEN
3979 0 : RgM_mulreal(GEN x, GEN y)
3980 : {
3981 0 : long i, j, k, l, lx = lg(x), ly = lg(y);
3982 0 : GEN z = cgetg(ly,t_MAT);
3983 0 : l = (lx == 1)? 1: lgcols(x);
3984 0 : for (j=1; j<ly; j++)
3985 : {
3986 0 : GEN zj = cgetg(l,t_COL), yj = gel(y,j);
3987 0 : gel(z,j) = zj;
3988 0 : for (i=1; i<l; i++)
3989 : {
3990 0 : pari_sp av = avma;
3991 0 : GEN c = mulreal(gcoeff(x,i,1),gel(yj,1));
3992 0 : for (k=2; k<lx; k++) c = gadd(c, mulreal(gcoeff(x,i,k),gel(yj,k)));
3993 0 : gel(zj, i) = gerepileupto(av, c);
3994 : }
3995 : }
3996 0 : return z;
3997 : }
3998 :
3999 : /* Compute Re(x * y), symetric result, x and y vectors of compatible
4000 : * dimensions; assume scalar entries */
4001 : GEN
4002 20853 : RgC_RgV_mulrealsym(GEN x, GEN y)
4003 : {
4004 20853 : long i, j, l = lg(x);
4005 20853 : GEN q = cgetg(l, t_MAT);
4006 83412 : for (j = 1; j < l; j++)
4007 : {
4008 62559 : gel(q,j) = cgetg(l,t_COL);
4009 187677 : for (i = 1; i <= j; i++)
4010 125118 : gcoeff(q,i,j) = gcoeff(q,j,i) = mulreal(gel(x,i), gel(y,j));
4011 : }
4012 20853 : return q;
4013 : }
4014 :
4015 : /*******************************************************************/
4016 : /* */
4017 : /* LOGICAL OPERATIONS */
4018 : /* */
4019 : /*******************************************************************/
4020 : static long
4021 96877652 : _egal_i(GEN x, GEN y)
4022 : {
4023 96877652 : x = simplify_shallow(x);
4024 96877652 : y = simplify_shallow(y);
4025 96877652 : if (typ(y) == t_INT)
4026 : {
4027 95984921 : if (equali1(y)) return gequal1(x);
4028 60997643 : if (equalim1(y)) return gequalm1(x);
4029 : }
4030 892731 : else if (typ(x) == t_INT)
4031 : {
4032 112 : if (equali1(x)) return gequal1(y);
4033 77 : if (equalim1(x)) return gequalm1(y);
4034 : }
4035 61758872 : return gequal(x, y);
4036 : }
4037 : static long
4038 96877652 : _egal(GEN x, GEN y)
4039 96877652 : { pari_sp av = avma; return gc_long(av, _egal_i(x, y)); }
4040 :
4041 : GEN
4042 6103563 : glt(GEN x, GEN y) { return gcmp(x,y)<0? gen_1: gen_0; }
4043 :
4044 : GEN
4045 7628238 : gle(GEN x, GEN y) { return gcmp(x,y)<=0? gen_1: gen_0; }
4046 :
4047 : GEN
4048 136428 : gge(GEN x, GEN y) { return gcmp(x,y)>=0? gen_1: gen_0; }
4049 :
4050 : GEN
4051 323299 : ggt(GEN x, GEN y) { return gcmp(x,y)>0? gen_1: gen_0; }
4052 :
4053 : GEN
4054 36913062 : geq(GEN x, GEN y) { return _egal(x,y)? gen_1: gen_0; }
4055 :
4056 : GEN
4057 59964590 : gne(GEN x, GEN y) { return _egal(x,y)? gen_0: gen_1; }
4058 :
4059 : GEN
4060 497490 : gnot(GEN x) { return gequal0(x)? gen_1: gen_0; }
4061 :
4062 : /*******************************************************************/
4063 : /* */
4064 : /* FORMAL SIMPLIFICATIONS */
4065 : /* */
4066 : /*******************************************************************/
4067 :
4068 : GEN
4069 10817 : geval_gp(GEN x, GEN t)
4070 : {
4071 10817 : long lx, i, tx = typ(x);
4072 : pari_sp av;
4073 : GEN y, z;
4074 :
4075 10817 : if (is_const_t(tx) || tx==t_VECSMALL) return gcopy(x);
4076 10796 : switch(tx)
4077 : {
4078 10789 : case t_STR:
4079 10789 : return localvars_read_str(GSTR(x),t);
4080 :
4081 0 : case t_POLMOD:
4082 0 : av = avma;
4083 0 : return gerepileupto(av, gmodulo(geval_gp(gel(x,2),t),
4084 0 : geval_gp(gel(x,1),t)));
4085 :
4086 7 : case t_POL:
4087 7 : lx=lg(x); if (lx==2) return gen_0;
4088 7 : z = fetch_var_value(varn(x),t);
4089 7 : if (!z) return RgX_copy(x);
4090 7 : av = avma; y = geval_gp(gel(x,lx-1),t);
4091 14 : for (i=lx-2; i>1; i--)
4092 7 : y = gadd(geval_gp(gel(x,i),t), gmul(z,y));
4093 7 : return gerepileupto(av, y);
4094 :
4095 0 : case t_SER:
4096 0 : pari_err_IMPL( "evaluation of a power series");
4097 :
4098 0 : case t_RFRAC:
4099 0 : av = avma;
4100 0 : return gerepileupto(av, gdiv(geval_gp(gel(x,1),t), geval_gp(gel(x,2),t)));
4101 :
4102 0 : case t_QFB: case t_VEC: case t_COL: case t_MAT:
4103 0 : pari_APPLY_same(geval_gp(gel(x,i),t));
4104 :
4105 0 : case t_CLOSURE:
4106 0 : if (closure_arity(x) || closure_is_variadic(x))
4107 0 : pari_err_IMPL("eval on functions with parameters");
4108 0 : return closure_evalres(x);
4109 : }
4110 0 : pari_err_TYPE("geval",x);
4111 : return NULL; /* LCOV_EXCL_LINE */
4112 : }
4113 : GEN
4114 0 : geval(GEN x) { return geval_gp(x,NULL); }
4115 :
4116 : GEN
4117 495232998 : simplify_shallow(GEN x)
4118 : {
4119 : long v, lx;
4120 : GEN y, z;
4121 495232998 : if (!x) pari_err_BUG("simplify, NULL input");
4122 :
4123 495232998 : switch(typ(x))
4124 : {
4125 412828678 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
4126 : case t_PADIC: case t_QFB: case t_LIST: case t_STR: case t_VECSMALL:
4127 : case t_CLOSURE: case t_ERROR: case t_INFINITY:
4128 412828678 : return x;
4129 688597 : case t_COMPLEX: return isintzero(gel(x,2))? gel(x,1): x;
4130 805 : case t_QUAD: return isintzero(gel(x,3))? gel(x,2): x;
4131 :
4132 198152 : case t_POLMOD: y = cgetg(3,t_POLMOD);
4133 198152 : z = gel(x,1); v = varn(z); z = simplify_shallow(z);
4134 198152 : if (typ(z) != t_POL || varn(z) != v) z = scalarpol_shallow(z, v);
4135 198152 : gel(y,1) = z;
4136 198152 : gel(y,2) = simplify_shallow(gel(x,2)); return y;
4137 :
4138 68766228 : case t_POL:
4139 68766228 : lx = lg(x);
4140 68766228 : if (lx==2) return gen_0;
4141 61175346 : if (lx==3) return simplify_shallow(gel(x,2));
4142 192704997 : pari_APPLY_pol(simplify_shallow(gel(x,i)));
4143 :
4144 3276 : case t_SER:
4145 1063979 : pari_APPLY_ser(simplify_shallow(gel(x,i)));
4146 :
4147 649465 : case t_RFRAC: y = cgetg(3,t_RFRAC);
4148 649465 : gel(y,1) = simplify_shallow(gel(x,1));
4149 649465 : z = simplify_shallow(gel(x,2));
4150 649465 : if (typ(z) != t_POL) return gdiv(gel(y,1), z);
4151 649465 : gel(y,2) = z; return y;
4152 :
4153 12097797 : case t_VEC: case t_COL: case t_MAT:
4154 63644918 : pari_APPLY_same(simplify_shallow(gel(x,i)));
4155 : }
4156 0 : pari_err_BUG("simplify_shallow, type unknown");
4157 : return NULL; /* LCOV_EXCL_LINE */
4158 : }
4159 :
4160 : GEN
4161 10811669 : simplify(GEN x)
4162 : {
4163 10811669 : pari_sp av = avma;
4164 10811669 : GEN y = simplify_shallow(x);
4165 10811669 : return av == avma ? gcopy(y): gerepilecopy(av, y);
4166 : }
4167 :
4168 : /*******************************************************************/
4169 : /* */
4170 : /* EVALUATION OF SOME SIMPLE OBJECTS */
4171 : /* */
4172 : /*******************************************************************/
4173 : /* q is a real symetric matrix, x a RgV. Horner-type evaluation of q(x)
4174 : * using (n^2+3n-2)/2 mul */
4175 : GEN
4176 16415 : qfeval(GEN q, GEN x)
4177 : {
4178 16415 : pari_sp av = avma;
4179 16415 : long i, j, l = lg(q);
4180 : GEN z;
4181 16415 : if (lg(x) != l) pari_err_DIM("qfeval");
4182 16408 : if (l==1) return gen_0;
4183 16408 : if (lgcols(q) != l) pari_err_DIM("qfeval");
4184 : /* l = lg(x) = lg(q) > 1 */
4185 16401 : z = gmul(gcoeff(q,1,1), gsqr(gel(x,1)));
4186 71981 : for (i=2; i<l; i++)
4187 : {
4188 55580 : GEN c = gel(q,i), s;
4189 55580 : if (isintzero(gel(x,i))) continue;
4190 41580 : s = gmul(gel(c,1), gel(x,1));
4191 124887 : for (j=2; j<i; j++) s = gadd(s, gmul(gel(c,j),gel(x,j)));
4192 41580 : s = gadd(gshift(s,1), gmul(gel(c,i),gel(x,i)));
4193 41580 : z = gadd(z, gmul(gel(x,i), s));
4194 : }
4195 16401 : return gerepileupto(av,z);
4196 : }
4197 :
4198 : static GEN
4199 351400 : qfbeval(GEN q, GEN z)
4200 : {
4201 351400 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3), x = gel(z,1), y = gel(z,2);
4202 351400 : pari_sp av = avma;
4203 351400 : A = gadd(gmul(x, gadd(gmul(a,x), gmul(b,y))), gmul(c, gsqr(y)));
4204 351400 : return gerepileupto(av, A);
4205 : }
4206 : static GEN
4207 7 : qfbevalb(GEN q, GEN z, GEN z2)
4208 : {
4209 7 : GEN A, a = gel(q,1), b = gel(q,2), c = gel(q,3);
4210 7 : GEN x = gel(z,1), y = gel(z,2);
4211 7 : GEN X = gel(z2,1), Y = gel(z2,2);
4212 7 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4213 7 : pari_sp av = avma;
4214 : /* a2 x X + b (x Y + X y) + c2 y Y */
4215 7 : A = gadd(gmul(x, gadd(gmul(a2,X), gmul(b,Y))),
4216 : gmul(y, gadd(gmul(c2,Y), gmul(b,X))));
4217 7 : return gerepileupto(av, gmul2n(A, -1));
4218 : }
4219 : GEN
4220 42 : qfb_apply_ZM(GEN q, GEN M)
4221 : {
4222 42 : pari_sp av = avma;
4223 42 : GEN A, B, C, a = gel(q,1), b = gel(q,2), c = gel(q,3);
4224 42 : GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
4225 42 : GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
4226 42 : GEN by = mulii(b,y), bt = mulii(b,t), bz = mulii(b,z);
4227 42 : GEN a2 = shifti(a,1), c2 = shifti(c,1);
4228 :
4229 42 : A = addii(mulii(x, addii(mulii(a,x), by)), mulii(c, sqri(y)));
4230 42 : B = addii(mulii(x, addii(mulii(a2,z), bt)),
4231 : mulii(y, addii(mulii(c2,t), bz)));
4232 42 : C = addii(mulii(z, addii(mulii(a,z), bt)), mulii(c, sqri(t)));
4233 42 : q = leafcopy(q);
4234 42 : gel(q,1) = A; gel(q,2) = B; gel(q,3) = C;
4235 42 : return gerepilecopy(av, q);
4236 : }
4237 :
4238 : static GEN
4239 351463 : qfnorm0(GEN q, GEN x)
4240 : {
4241 351463 : if (!q) switch(typ(x))
4242 : {
4243 7 : case t_VEC: case t_COL: return RgV_dotsquare(x);
4244 7 : case t_MAT: return gram_matrix(x);
4245 7 : default: pari_err_TYPE("qfeval",x);
4246 : }
4247 351442 : switch(typ(q))
4248 : {
4249 28 : case t_MAT: break;
4250 351407 : case t_QFB:
4251 351407 : if (lg(x) == 3) switch(typ(x))
4252 : {
4253 351400 : case t_VEC:
4254 351400 : case t_COL: return qfbeval(q,x);
4255 7 : case t_MAT: if (RgM_is_ZM(x)) return qfb_apply_ZM(q,x);
4256 0 : default: pari_err_TYPE("qfeval",x);
4257 : }
4258 7 : default: pari_err_TYPE("qfeval",q);
4259 : }
4260 28 : switch(typ(x))
4261 : {
4262 21 : case t_VEC: case t_COL: break;
4263 7 : case t_MAT: return qf_apply_RgM(q, x);
4264 0 : default: pari_err_TYPE("qfeval",x);
4265 : }
4266 21 : return qfeval(q,x);
4267 : }
4268 : /* obsolete, use qfeval0 */
4269 : GEN
4270 0 : qfnorm(GEN x, GEN q) { return qfnorm0(q,x); }
4271 :
4272 : /* assume q is square, x~ * q * y using n^2+n mul */
4273 : GEN
4274 21 : qfevalb(GEN q, GEN x, GEN y)
4275 : {
4276 21 : pari_sp av = avma;
4277 21 : long l = lg(q);
4278 21 : if (lg(x) != l || lg(y) != l) pari_err_DIM("qfevalb");
4279 14 : return gerepileupto(av, RgV_dotproduct(RgV_RgM_mul(x,q), y));
4280 : }
4281 :
4282 : /* obsolete, use qfeval0 */
4283 : GEN
4284 0 : qfbil(GEN x, GEN y, GEN q)
4285 : {
4286 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfbil",x);
4287 0 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfbil",y);
4288 0 : if (!q) {
4289 0 : if (lg(x) != lg(y)) pari_err_DIM("qfbil");
4290 0 : return RgV_dotproduct(x,y);
4291 : }
4292 0 : if (typ(q) != t_MAT) pari_err_TYPE("qfbil",q);
4293 0 : return qfevalb(q,x,y);
4294 : }
4295 : GEN
4296 351519 : qfeval0(GEN q, GEN x, GEN y)
4297 : {
4298 351519 : if (!y) return qfnorm0(q, x);
4299 56 : if (!is_vec_t(typ(x))) pari_err_TYPE("qfeval",x);
4300 42 : if (!is_vec_t(typ(y))) pari_err_TYPE("qfeval",y);
4301 42 : if (!q) {
4302 14 : if (lg(x) != lg(y)) pari_err_DIM("qfeval");
4303 7 : return RgV_dotproduct(x,y);
4304 : }
4305 28 : switch(typ(q))
4306 : {
4307 21 : case t_MAT: break;
4308 7 : case t_QFB:
4309 7 : if (lg(x) == 3 && lg(y) == 3) return qfbevalb(q,x,y);
4310 0 : default: pari_err_TYPE("qfeval",q);
4311 : }
4312 21 : return qfevalb(q,x,y);
4313 : }
4314 :
4315 : /* q a hermitian complex matrix, x a RgV */
4316 : GEN
4317 0 : hqfeval(GEN q, GEN x)
4318 : {
4319 0 : pari_sp av = avma;
4320 0 : long i, j, l = lg(q);
4321 : GEN z, xc;
4322 :
4323 0 : if (lg(x) != l) pari_err_DIM("hqfeval");
4324 0 : if (l==1) return gen_0;
4325 0 : if (lgcols(q) != l) pari_err_DIM("hqfeval");
4326 0 : if (l == 2) return gerepileupto(av, gmul(gcoeff(q,1,1), gnorm(gel(x,1))));
4327 : /* l = lg(x) = lg(q) > 2 */
4328 0 : xc = conj_i(x);
4329 0 : z = mulreal(gcoeff(q,2,1), gmul(gel(x,2),gel(xc,1)));
4330 0 : for (i=3;i<l;i++)
4331 0 : for (j=1;j<i;j++)
4332 0 : z = gadd(z, mulreal(gcoeff(q,i,j), gmul(gel(x,i),gel(xc,j))));
4333 0 : z = gshift(z,1);
4334 0 : for (i=1;i<l;i++) z = gadd(z, gmul(gcoeff(q,i,i), gnorm(gel(x,i))));
4335 0 : return gerepileupto(av,z);
4336 : }
4337 :
4338 : static void
4339 498135 : init_qf_apply(GEN q, GEN M, long *l)
4340 : {
4341 498135 : long k = lg(M);
4342 498135 : *l = lg(q);
4343 498135 : if (*l == 1) { if (k == 1) return; }
4344 498128 : else { if (k != 1 && lgcols(M) == *l) return; }
4345 0 : pari_err_DIM("qf_apply_RgM");
4346 : }
4347 : /* Return X = M'.q.M, assuming q is a symetric matrix and M is a
4348 : * matrix of compatible dimensions. X_ij are X_ji identical, not copies */
4349 : GEN
4350 7461 : qf_apply_RgM(GEN q, GEN M)
4351 : {
4352 7461 : pari_sp av = avma;
4353 7461 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4354 7461 : return gerepileupto(av, RgM_transmultosym(M, RgM_mul(q, M)));
4355 : }
4356 : GEN
4357 490674 : qf_apply_ZM(GEN q, GEN M)
4358 : {
4359 490674 : pari_sp av = avma;
4360 490674 : long l; init_qf_apply(q, M, &l); if (l == 1) return cgetg(1, t_MAT);
4361 490667 : return gerepileupto(av, ZM_transmultosym(M, ZM_mul(q, M)));
4362 : }
4363 :
4364 : GEN
4365 2855160 : poleval(GEN x, GEN y)
4366 : {
4367 2855160 : long i, j, imin, tx = typ(x);
4368 2855160 : pari_sp av0 = avma, av;
4369 : GEN t, t2, r, s;
4370 :
4371 2855160 : if (is_scalar_t(tx)) return gcopy(x);
4372 2685161 : switch(tx)
4373 : {
4374 2636504 : case t_POL:
4375 2636504 : i = lg(x)-1; imin = 2; break;
4376 :
4377 1568 : case t_RFRAC:
4378 1568 : return gerepileupto(av0, gdiv(poleval(gel(x,1),y),
4379 1568 : poleval(gel(x,2),y)));
4380 :
4381 47089 : case t_VEC: case t_COL:
4382 47089 : i = lg(x)-1; imin = 1; break;
4383 0 : default: pari_err_TYPE("poleval",x);
4384 : return NULL; /* LCOV_EXCL_LINE */
4385 : }
4386 2683593 : if (i<=imin) return (i==imin)? gmul(gel(x,imin),Rg_get_1(y)): Rg_get_0(y);
4387 2543368 : if (isintzero(y)) return gcopy(gel(x,imin));
4388 :
4389 2536400 : t = gel(x,i); i--;
4390 2536400 : if (typ(y)!=t_COMPLEX)
4391 : {
4392 : #if 0 /* standard Horner's rule */
4393 : for ( ; i>=imin; i--)
4394 : t = gadd(gmul(t,y),gel(x,i));
4395 : #endif
4396 : /* specific attention to sparse polynomials */
4397 18328505 : for ( ; i>=imin; i=j-1)
4398 : {
4399 18574738 : for (j=i; isexactzero(gel(x,j)); j--)
4400 2696912 : if (j==imin)
4401 : {
4402 970644 : if (i!=j) y = gpowgs(y, i-j+1);
4403 970644 : return gerepileupto(av0, gmul(t,y));
4404 : }
4405 15877827 : r = (i==j)? y: gpowgs(y, i-j+1);
4406 15877827 : t = gadd(gmul(t,r), gel(x,j));
4407 15877752 : if (gc_needed(av0,2))
4408 : {
4409 55 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4410 55 : t = gerepileupto(av0, t);
4411 : }
4412 : }
4413 1480035 : return gerepileupto(av0, t);
4414 : }
4415 :
4416 85647 : t2 = gel(x,i); i--; r = gtrace(y); s = gneg_i(gnorm(y));
4417 85647 : av = avma;
4418 4098504 : for ( ; i>=imin; i--)
4419 : {
4420 4012857 : GEN p = gadd(t2, gmul(r, t));
4421 4012857 : t2 = gadd(gel(x,i), gmul(s, t)); t = p;
4422 4012857 : if (gc_needed(av0,2))
4423 : {
4424 0 : if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
4425 0 : gerepileall(av, 2, &t, &t2);
4426 : }
4427 : }
4428 85647 : return gerepileupto(av0, gadd(t2, gmul(y, t)));
4429 : }
4430 :
4431 : /* Evaluate a polynomial using Horner. Unstable!
4432 : * If ui != NULL, ui = 1/u, evaluate P(1/u)*u^(deg P): useful for |u|>1 */
4433 : GEN
4434 2223266 : RgX_cxeval(GEN T, GEN u, GEN ui)
4435 : {
4436 2223266 : pari_sp ltop = avma;
4437 : GEN S;
4438 2223266 : long n, lim = lg(T)-1;
4439 2223266 : if (lim == 1) return gen_0;
4440 2223266 : if (lim == 2) return gcopy(gel(T,2));
4441 2222195 : if (!ui)
4442 : {
4443 1495120 : n = lim; S = gel(T,n);
4444 14101947 : for (n--; n >= 2; n--) S = gadd(gmul(u,S), gel(T,n));
4445 : }
4446 : else
4447 : {
4448 727075 : n = 2; S = gel(T,2);
4449 4369667 : for (n++; n <= lim; n++) S = gadd(gmul(ui,S), gel(T,n));
4450 727073 : S = gmul(gpowgs(u, lim-2), S);
4451 : }
4452 2221940 : return gerepileupto(ltop, S);
4453 : }
4454 :
4455 : GEN
4456 63 : RgXY_cxevalx(GEN x, GEN u, GEN ui)
4457 : {
4458 427 : pari_APPLY_pol(typ(gel(x,i))==t_POL? RgX_cxeval(gel(x,i), u, ui): gel(x,i));
4459 : }
|
