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 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_pol
19 :
20 : /*******************************************************************/
21 : /* */
22 : /* GENERIC */
23 : /* */
24 : /*******************************************************************/
25 :
26 : /* Return optimal parameter l for the evaluation of n/m polynomials of degree d
27 : Fractional values can be used if the evaluations are done with different
28 : accuracies, and thus have different weights.
29 : */
30 : long
31 19692461 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19692461 : long pold=1, rold=n*(d-1);
35 92698601 : for(p=2; p<=d; p++)
36 : {
37 73006140 : r = m*(p-1) + n*((d-1)/p);
38 73006140 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19692461 : return pold;
41 : }
42 :
43 : static GEN
44 14485901 : gen_RgXQ_eval_powers(GEN P, GEN V, long a, long n, void *E, const struct bb_algebra *ff,
45 : GEN cmul(void *E, GEN P, long a, GEN x))
46 : {
47 14485901 : pari_sp av = avma;
48 : long i;
49 14485901 : GEN z = cmul(E,P,a,ff->one(E));
50 14458481 : if (!z) z = gen_0;
51 75529619 : for (i=1; i<=n; i++)
52 : {
53 61069716 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 61091240 : if (t) {
55 46165020 : z = ff->add(E, z, t);
56 46133890 : if (gc_needed(av,2)) z = gerepileupto(av, z);
57 : }
58 : }
59 14459903 : return ff->red(E,z);
60 : }
61 :
62 : /* Brent & Kung
63 : * (Fast algorithms for manipulating formal power series, JACM 25:581-595, 1978)
64 : *
65 : * V as output by FpXQ_powers(x,l,T,p). For optimal performance, l is as given
66 : * by brent_kung_optpow */
67 : GEN
68 12437825 : gen_bkeval_powers(GEN P, long d, GEN V, void *E, const struct bb_algebra *ff,
69 : GEN cmul(void *E, GEN P, long a, GEN x))
70 : {
71 12437825 : pari_sp av = avma;
72 12437825 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12437825 : if (d < 0) return ff->zero(E);
76 11830847 : if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1147717 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1147717 : if (DEBUGLEVEL>=8)
79 : {
80 0 : long cnt = 1 + (d - l) / (l-1);
81 0 : err_printf("RgX_RgXQV_eval(%ld/%ld): %ld RgXQ_mul\n", d, l-1, cnt);
82 : }
83 1147717 : d -= l;
84 1147717 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2654385 : while (d >= l-1)
86 : {
87 1505520 : d -= l-1;
88 1505520 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1505568 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1505484 : if (gc_needed(av,2))
91 105 : z = gerepileupto(av, z);
92 : }
93 1148865 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1148879 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1148902 : return gerepileupto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 866648 : gen_bkeval(GEN Q, long d, GEN x, int use_sqr, void *E, const struct bb_algebra *ff,
100 : GEN cmul(void *E, GEN P, long a, GEN x))
101 : {
102 866648 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 866648 : if (d < 0) return ff->zero(E);
106 865570 : rtd = (long) sqrt((double)d);
107 865570 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 865587 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 865579 : return gerepileupto(av, z);
110 : }
111 :
112 : static GEN
113 2027153 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19227235 : _gen_add(void *E, GEN x, GEN y) { (void)E; return gadd(x, y); }
116 : static GEN
117 0 : _gen_sub(void *E, GEN x, GEN y) { (void)E; return gsub(x, y); }
118 : static GEN
119 1842127 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 599490 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2068161 : _gen_one(void *E) { (void)E; return gen_1; }
124 : static GEN
125 24184 : _gen_zero(void *E) { (void)E; return gen_0; }
126 :
127 : static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
128 : _gen_mul, _gen_sqr,_gen_one,_gen_zero };
129 :
130 : static GEN
131 512818 : _gen_cmul(void *E, GEN P, long a, GEN x)
132 512818 : {(void)E; return gmul(gel(P,a+2), x);}
133 :
134 : GEN
135 166768 : RgX_RgV_eval(GEN Q, GEN x)
136 : {
137 166768 : return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
138 : }
139 :
140 : GEN
141 0 : RgX_Rg_eval_bk(GEN Q, GEN x)
142 : {
143 0 : return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
144 : }
145 :
146 : GEN
147 2947 : RgXV_RgV_eval(GEN Q, GEN x)
148 : {
149 2947 : long i, l = lg(Q), vQ = gvar(Q);
150 2947 : GEN v = cgetg(l, t_VEC);
151 248311 : for (i = 1; i < l; i++)
152 : {
153 245364 : GEN Qi = gel(Q, i);
154 245364 : gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
155 : }
156 2947 : return v;
157 : }
158 :
159 : GEN
160 666094 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 666094 : pari_sp av = avma, btop;
163 666094 : long i, d = degpol(P), o;
164 : GEN s;
165 666093 : if (signe(P)==0) return pol_0(varn(P));
166 666093 : s = gel(P, d+2);
167 666093 : if (d == 0) return gcopy(s);
168 662817 : o = RgX_deflate_order(P);
169 662840 : if (o > 1) A = gpowgs(A, o);
170 662840 : btop = avma;
171 2304844 : for (i = d-o; i >= 0; i-=o)
172 : {
173 1642005 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 1642004 : if (gc_needed(btop,1))
175 : {
176 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
177 13 : s = gerepileupto(btop, s);
178 : }
179 : }
180 662839 : return gerepileupto(av, s);
181 : }
182 :
183 : GEN
184 1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
185 : {
186 1652 : pari_sp av = avma;
187 1652 : long i, d = degpol(P), v = varn(A);
188 : GEN s;
189 1652 : if (signe(P)==0) return pol_0(v);
190 1652 : if (d == 0) return scalarpol(gel(P, d+2), v);
191 1232 : s = scalarpol_shallow(gel(P, d+2), v);
192 4963 : for (i = d-1; i >= 0; i--)
193 : {
194 3731 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
195 3731 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
196 3731 : if (gc_needed(av,1))
197 : {
198 0 : if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
199 0 : s = gerepileupto(av, s);
200 : }
201 : }
202 1232 : return gerepileupto(av, s);
203 : }
204 :
205 : const struct bb_algebra *
206 281035 : get_Rg_algebra(void)
207 : {
208 281035 : return &Rg_algebra;
209 : }
210 :
211 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
212 :
213 : static GEN
214 30352 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
215 : {
216 : (void) E;
217 30352 : return RgX_divrem(x, y, r);
218 : }
219 :
220 : GEN
221 6923 : RgX_digits(GEN x, GEN T)
222 : {
223 6923 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
224 6923 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
225 : }
226 :
227 : /*******************************************************************/
228 : /* */
229 : /* RgX */
230 : /* */
231 : /*******************************************************************/
232 :
233 : long
234 24741740 : RgX_equal(GEN x, GEN y)
235 : {
236 24741740 : long i = lg(x);
237 :
238 24741740 : if (i != lg(y)) return 0;
239 106408107 : for (i--; i > 1; i--)
240 81982916 : if (!gequal(gel(x,i),gel(y,i))) return 0;
241 24425191 : return 1;
242 : }
243 :
244 : /* Returns 1 in the base ring over which x is defined */
245 : /* HACK: this also works for t_SER */
246 : GEN
247 156674693 : Rg_get_1(GEN x)
248 : {
249 : GEN p, T;
250 156674693 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
251 156674693 : if (RgX_type_is_composite(tx))
252 11628293 : RgX_type_decode(tx, &i /*junk*/, &tx);
253 156674693 : switch(tx)
254 : {
255 795550 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
256 973 : case t_PADIC: return cvtop(gen_1, p, lx);
257 5516 : case t_FFELT: return FF_1(T);
258 155872654 : default: return gen_1;
259 : }
260 : }
261 : /* Returns 0 in the base ring over which x is defined */
262 : /* HACK: this also works for t_SER */
263 : GEN
264 6985678 : Rg_get_0(GEN x)
265 : {
266 : GEN p, T;
267 6985678 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
268 6985678 : if (RgX_type_is_composite(tx))
269 57274 : RgX_type_decode(tx, &i /*junk*/, &tx);
270 6985678 : switch(tx)
271 : {
272 532 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
273 42 : case t_PADIC: return zeropadic(p, lx);
274 231 : case t_FFELT: return FF_zero(T);
275 6984873 : default: return gen_0;
276 : }
277 : }
278 :
279 : GEN
280 6279 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
281 : {
282 6279 : long i, n = degpol(P);
283 : GEN z, dz, dP;
284 6279 : if (n < 0) return gen_0;
285 6279 : P = Q_remove_denom(P, &dP);
286 6279 : z = gel(P,2); if (n == 0) return icopy(z);
287 3619 : if (dV) z = mulii(dV, z); /* V[1] = dV */
288 3619 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
289 7210 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
290 3619 : dz = mul_denom(dP, dV);
291 3619 : return dz? RgX_Rg_div(z, dz): z;
292 : }
293 :
294 : /* Return P(h * x), not memory clean */
295 : GEN
296 19856 : RgX_unscale(GEN P, GEN h)
297 : {
298 19856 : long i, l = lg(P);
299 19856 : GEN hi = gen_1, Q = cgetg(l, t_POL);
300 19856 : Q[1] = P[1];
301 19856 : if (l == 2) return Q;
302 19828 : gel(Q,2) = gcopy(gel(P,2));
303 42494 : for (i=3; i<l; i++)
304 : {
305 22666 : hi = gmul(hi,h);
306 22667 : gel(Q,i) = gmul(gel(P,i), hi);
307 : }
308 19828 : return Q;
309 : }
310 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
311 : GEN
312 1463527 : ZX_z_unscale(GEN P, long h)
313 : {
314 1463527 : long i, l = lg(P);
315 1463527 : GEN Q = cgetg(l, t_POL);
316 1463526 : Q[1] = P[1];
317 1463526 : if (l == 2) return Q;
318 1374980 : gel(Q,2) = gel(P,2);
319 1374980 : if (l == 3) return Q;
320 1348674 : if (h == -1)
321 236354 : for (i = 3; i < l; i++)
322 : {
323 195943 : gel(Q,i) = negi(gel(P,i));
324 195943 : if (++i == l) break;
325 144991 : gel(Q,i) = gel(P,i);
326 : }
327 : else
328 : {
329 : GEN hi;
330 1257311 : gel(Q,3) = mulis(gel(P,3), h);
331 1257310 : hi = sqrs(h);
332 5397290 : for (i = 4; i < l; i++)
333 : {
334 4139981 : gel(Q,i) = mulii(gel(P,i), hi);
335 4139982 : if (i != l-1) hi = mulis(hi,h);
336 : }
337 : }
338 1348672 : return Q;
339 : }
340 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
341 : GEN
342 1009704 : ZX_unscale(GEN P, GEN h)
343 : {
344 : long i, l;
345 : GEN Q, hi;
346 1009704 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
347 881 : l = lg(P); Q = cgetg(l, t_POL);
348 881 : Q[1] = P[1];
349 881 : if (l == 2) return Q;
350 881 : gel(Q,2) = gel(P,2);
351 881 : if (l == 3) return Q;
352 881 : hi = h;
353 881 : gel(Q,3) = mulii(gel(P,3), hi);
354 2741 : for (i = 4; i < l; i++)
355 : {
356 1860 : hi = mulii(hi,h);
357 1860 : gel(Q,i) = mulii(gel(P,i), hi);
358 : }
359 881 : return Q;
360 : }
361 : /* P a ZX. Return P(x << n), not memory clean */
362 : GEN
363 926606 : ZX_unscale2n(GEN P, long n)
364 : {
365 926606 : long i, ni = n, l = lg(P);
366 926606 : GEN Q = cgetg(l, t_POL);
367 926608 : Q[1] = P[1];
368 926608 : if (l == 2) return Q;
369 926608 : gel(Q,2) = gel(P,2);
370 926608 : if (l == 3) return Q;
371 926608 : gel(Q,3) = shifti(gel(P,3), ni);
372 3168534 : for (i=4; i<l; i++)
373 : {
374 2241979 : ni += n;
375 2241979 : gel(Q,i) = shifti(gel(P,i), ni);
376 : }
377 926555 : return Q;
378 : }
379 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
380 : GEN
381 12997 : ZX_unscale_div(GEN P, GEN h)
382 : {
383 12997 : long i, l = lg(P);
384 12997 : GEN hi, Q = cgetg(l, t_POL);
385 12997 : Q[1] = P[1];
386 12997 : if (l == 2) return Q;
387 12997 : gel(Q,2) = diviiexact(gel(P,2), h);
388 12997 : if (l == 3) return Q;
389 12997 : gel(Q,3) = gel(P,3);
390 12997 : if (l == 4) return Q;
391 12997 : hi = h;
392 12997 : gel(Q,4) = mulii(gel(P,4), hi);
393 64121 : for (i=5; i<l; i++)
394 : {
395 51124 : hi = mulii(hi,h);
396 51124 : gel(Q,i) = mulii(gel(P,i), hi);
397 : }
398 12997 : return Q;
399 : }
400 : /* P(h*X) / h^k, assuming the result is a ZX */
401 : GEN
402 966 : ZX_unscale_divpow(GEN P, GEN h, long k)
403 : {
404 966 : long i, j, l = lg(P);
405 966 : GEN H, Q = cgetg(l, t_POL);
406 966 : Q[1] = P[1]; if (l == 2) return Q;
407 966 : H = gpowers(h, maxss(k, l - 3 - k));
408 3864 : for (i = 2, j = k+1; j > 1 && i < l; i++)
409 2898 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
410 966 : if (i == l) return Q;
411 966 : gel(Q, i) = gel(P, i); i++;
412 3542 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
413 966 : return Q;
414 : }
415 :
416 : GEN
417 5845 : RgXV_unscale(GEN x, GEN h)
418 : {
419 5845 : if (isint1(h)) return gcopy(x);
420 16739 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
421 : }
422 :
423 : /* Return h^degpol(P) P(x / h), not memory clean */
424 : GEN
425 4334226 : RgX_rescale(GEN P, GEN h)
426 : {
427 4334226 : long i, l = lg(P);
428 4334226 : GEN Q = cgetg(l,t_POL), hi = h;
429 4334224 : gel(Q,l-1) = gel(P,l-1);
430 11403664 : for (i=l-2; i>=2; i--)
431 : {
432 11401088 : gel(Q,i) = gmul(gel(P,i), hi);
433 11400895 : if (i == 2) break;
434 7069119 : hi = gmul(hi,h);
435 : }
436 4334352 : Q[1] = P[1]; return Q;
437 : }
438 :
439 : GEN
440 2401 : RgXV_rescale(GEN x, GEN h)
441 : {
442 2401 : if (isint1(h)) return RgX_copy(x);
443 16086 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
444 : }
445 :
446 : /* A(X^d) --> A(X) */
447 : GEN
448 984818 : RgX_deflate(GEN x0, long d)
449 : {
450 : GEN z, y, x;
451 984818 : long i,id, dy, dx = degpol(x0);
452 984815 : if (d == 1 || dx <= 0) return leafcopy(x0);
453 374984 : dy = dx/d;
454 374984 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
455 374980 : z = y + 2;
456 374980 : x = x0+ 2;
457 1509346 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
458 374980 : return y;
459 : }
460 :
461 : GEN
462 18732 : RgX_homogenize(GEN P, long v)
463 : {
464 : long i, l, d;
465 18732 : GEN Q = cgetg_copy(P, &l);
466 18732 : Q[1] = P[1]; d = l-3;
467 165998 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
468 18732 : return Q;
469 : }
470 :
471 : /* F a t_RFRAC */
472 : long
473 105 : rfrac_deflate_order(GEN F)
474 : {
475 105 : GEN N = gel(F,1), D = gel(F,2);
476 105 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
477 105 : if (m == 1) return 1;
478 42 : if (typ(N) == t_POL && varn(N) == varn(D))
479 28 : m = cgcd(m, RgX_deflate_order(N));
480 42 : return m;
481 : }
482 : /* F a t_RFRAC */
483 : GEN
484 105 : rfrac_deflate_max(GEN F, long *m)
485 : {
486 105 : *m = rfrac_deflate_order(F);
487 105 : return rfrac_deflate(F, *m);
488 : }
489 : /* F a t_RFRAC */
490 : GEN
491 105 : rfrac_deflate(GEN F, long m)
492 : {
493 105 : GEN N = gel(F,1), D = gel(F,2);
494 105 : if (m == 1) return F;
495 42 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
496 42 : D = RgX_deflate(D, m); return mkrfrac(N, D);
497 : }
498 :
499 : /* return x0(X^d) */
500 : GEN
501 798166 : RgX_inflate(GEN x0, long d)
502 : {
503 798166 : long i, id, dy, dx = degpol(x0);
504 798166 : GEN x = x0 + 2, z, y;
505 798166 : if (dx <= 0) return leafcopy(x0);
506 790564 : dy = dx*d;
507 790564 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
508 790564 : z = y + 2;
509 27844483 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
510 11404749 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
511 790564 : return y;
512 : }
513 :
514 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
515 : static GEN
516 5956781 : RgX_translate_basecase(GEN P, GEN c)
517 : {
518 5956781 : pari_sp av = avma;
519 : GEN Q, R;
520 : long i, k, n;
521 :
522 5956781 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
523 5954810 : Q = leafcopy(P);
524 5954807 : R = Q+2; n = degpol(P);
525 5954805 : if (isint1(c))
526 : {
527 26178 : for (i=1; i<=n; i++)
528 : {
529 135148 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
530 22336 : if (gc_needed(av,2))
531 : {
532 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
533 0 : Q = gerepilecopy(av, Q); R = Q+2;
534 : }
535 : }
536 : }
537 5950958 : else if (isintm1(c))
538 : {
539 901306 : for (i=1; i<=n; i++)
540 : {
541 2855358 : for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
542 683882 : if (gc_needed(av,2))
543 : {
544 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
545 0 : Q = gerepilecopy(av, Q); R = Q+2;
546 : }
547 : }
548 : }
549 : else
550 : {
551 19600501 : for (i=1; i<=n; i++)
552 : {
553 46008474 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
554 13866966 : if (gc_needed(av,2))
555 : {
556 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
557 0 : Q = gerepilecopy(av, Q); R = Q+2;
558 : }
559 : }
560 : }
561 5954476 : return gerepilecopy(av, Q);
562 : }
563 : GEN
564 5957212 : RgX_translate(GEN P, GEN c)
565 : {
566 5957212 : pari_sp av = avma;
567 5957212 : long n = degpol(P);
568 5957201 : if (n < 40)
569 5956780 : return RgX_translate_basecase(P, c);
570 : else
571 : {
572 421 : long d = n >> 1;
573 421 : GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
574 420 : GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
575 420 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
576 420 : return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
577 : }
578 : }
579 :
580 : /* P(ax + b) */
581 : GEN
582 0 : RgX_affine(GEN P, GEN a, GEN b)
583 : {
584 0 : if (signe(b)) P = RgX_translate(P, b);
585 0 : return RgX_unscale(P, a);
586 : }
587 :
588 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
589 : GEN
590 27529 : RgXQX_translate(GEN P, GEN c, GEN T)
591 : {
592 27529 : pari_sp av = avma;
593 : GEN Q, R;
594 : long i, k, n;
595 :
596 27529 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
597 27186 : Q = leafcopy(P);
598 27186 : R = Q+2; n = degpol(P);
599 89759 : for (i=1; i<=n; i++)
600 : {
601 270605 : for (k=n-i; k<n; k++)
602 : {
603 208032 : pari_sp av2 = avma;
604 208032 : gel(R,k) = gerepileupto(av2,
605 208032 : RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
606 : }
607 62573 : if (gc_needed(av,2))
608 : {
609 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
610 0 : Q = gerepilecopy(av, Q); R = Q+2;
611 : }
612 : }
613 27186 : return gerepilecopy(av, Q);
614 : }
615 :
616 : /********************************************************************/
617 : /** **/
618 : /** CONVERSIONS **/
619 : /** (not memory clean) **/
620 : /** **/
621 : /********************************************************************/
622 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
623 : * but everything else is */
624 : static GEN
625 167398 : QXQ_to_mod(GEN x, GEN T)
626 : {
627 : long d;
628 167398 : switch(typ(x))
629 : {
630 62577 : case t_INT: return icopy(x);
631 2079 : case t_FRAC: return gcopy(x);
632 102742 : case t_POL:
633 102742 : d = degpol(x);
634 102742 : if (d < 0) return gen_0;
635 100474 : if (d == 0) return gcopy(gel(x,2));
636 98335 : return mkpolmod(RgX_copy(x), T);
637 0 : default: pari_err_TYPE("QXQ_to_mod",x);
638 : return NULL;/* LCOV_EXCL_LINE */
639 : }
640 : }
641 : /* pure shallow version */
642 : GEN
643 780217 : QXQ_to_mod_shallow(GEN x, GEN T)
644 : {
645 : long d;
646 780217 : switch(typ(x))
647 : {
648 500747 : case t_INT:
649 500747 : case t_FRAC: return x;
650 279470 : case t_POL:
651 279470 : d = degpol(x);
652 279470 : if (d < 0) return gen_0;
653 232548 : if (d == 0) return gel(x,2);
654 216018 : return mkpolmod(x, T);
655 0 : default: pari_err_TYPE("QXQ_to_mod",x);
656 : return NULL;/* LCOV_EXCL_LINE */
657 : }
658 : }
659 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
660 : * Not memory clean because T not copied, but everything else is */
661 : static GEN
662 37366 : QXQX_to_mod(GEN z, GEN T)
663 : {
664 37366 : long i,l = lg(z);
665 37366 : GEN x = cgetg(l,t_POL);
666 186550 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
667 37366 : x[1] = z[1]; return normalizepol_lg(x,l);
668 : }
669 : /* pure shallow version */
670 : GEN
671 181115 : QXQX_to_mod_shallow(GEN z, GEN T)
672 : {
673 181115 : long i,l = lg(z);
674 181115 : GEN x = cgetg(l,t_POL);
675 878732 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
676 181115 : x[1] = z[1]; return normalizepol_lg(x,l);
677 : }
678 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
679 : GEN
680 12586 : QXQXV_to_mod(GEN V, GEN T)
681 : {
682 12586 : long i, l = lg(V);
683 12586 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
684 49952 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
685 12586 : return z;
686 : }
687 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
688 : GEN
689 18227 : QXQV_to_mod(GEN V, GEN T)
690 : {
691 18227 : long i, l = lg(V);
692 18227 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
693 36441 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
694 18227 : return z;
695 : }
696 :
697 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
698 : GEN
699 14336 : QXQC_to_mod_shallow(GEN V, GEN T)
700 : {
701 14336 : long i, l = lg(V);
702 14336 : GEN z = cgetg(l, t_COL);
703 96936 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
704 14336 : return z;
705 : }
706 :
707 : GEN
708 6573 : QXQM_to_mod_shallow(GEN V, GEN T)
709 : {
710 6573 : long i, l = lg(V);
711 6573 : GEN z = cgetg(l, t_MAT);
712 20909 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
713 6573 : return z;
714 : }
715 :
716 : GEN
717 6216999 : RgX_renormalize_lg(GEN x, long lx)
718 : {
719 : long i;
720 9240708 : for (i = lx-1; i>1; i--)
721 8788503 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
722 6216999 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
723 6216999 : setlg(x, i+1); setsigne(x, i != 1); return x;
724 : }
725 :
726 : GEN
727 1501916 : RgV_to_RgX(GEN x, long v)
728 : {
729 1501916 : long i, k = lg(x);
730 : GEN p;
731 :
732 4172978 : while (--k && gequal0(gel(x,k)));
733 1501916 : if (!k) return pol_0(v);
734 1494645 : i = k+2; p = cgetg(i,t_POL);
735 1494646 : p[1] = evalsigne(1) | evalvarn(v);
736 10270287 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
737 1494646 : return p;
738 : }
739 : GEN
740 190005 : RgV_to_RgX_reverse(GEN x, long v)
741 : {
742 190005 : long j, k, l = lg(x);
743 : GEN p;
744 :
745 191573 : for (k = 1; k < l; k++)
746 191573 : if (!gequal0(gel(x,k))) break;
747 190005 : if (k == l) return pol_0(v);
748 190005 : k -= 1;
749 190005 : l -= k;
750 190005 : x += k;
751 190005 : p = cgetg(l+1,t_POL);
752 190005 : p[1] = evalsigne(1) | evalvarn(v);
753 999045 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
754 190005 : return p;
755 : }
756 :
757 : /* return the (N-dimensional) vector of coeffs of p */
758 : GEN
759 13342621 : RgX_to_RgC(GEN x, long N)
760 : {
761 : long i, l;
762 : GEN z;
763 13342621 : l = lg(x)-1; x++;
764 13342621 : if (l > N+1) l = N+1; /* truncate higher degree terms */
765 13342621 : z = cgetg(N+1,t_COL);
766 80524020 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
767 24349278 : for ( ; i<=N; i++) gel(z,i) = gen_0;
768 13342672 : return z;
769 : }
770 : GEN
771 1349242 : Rg_to_RgC(GEN x, long N)
772 : {
773 1349242 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
774 : }
775 :
776 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
777 : GEN
778 302864 : RgM_to_RgXV(GEN x, long v)
779 1292825 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
780 : GEN
781 7202 : RgM_to_RgXV_reverse(GEN x, long v)
782 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
783 :
784 : /* matrix whose entries are given by the coeffs of the polynomials in
785 : * vector v (considered as degree n-1 polynomials) */
786 : GEN
787 334790 : RgV_to_RgM(GEN x, long n)
788 1681544 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
789 :
790 : GEN
791 77924 : RgXV_to_RgM(GEN x, long n)
792 394797 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
793 :
794 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
795 : GEN
796 23466 : RgM_to_RgXX(GEN x, long v,long w)
797 : {
798 23466 : long j, lx = lg(x);
799 23466 : GEN y = cgetg(lx+1, t_POL);
800 23466 : y[1] = evalsigne(1) | evalvarn(v);
801 23466 : y++;
802 130418 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
803 23466 : return normalizepol_lg(--y, lx+1);
804 : }
805 :
806 : /* matrix whose entries are given by the coeffs of the polynomial v in
807 : * two variables (considered as degree n-1 polynomials) */
808 : GEN
809 322 : RgXX_to_RgM(GEN v, long n)
810 : {
811 322 : long j, N = lg(v)-1;
812 322 : GEN y = cgetg(N, t_MAT);
813 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
814 322 : return y;
815 : }
816 :
817 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
818 : GEN
819 31304 : RgXY_swapspec(GEN x, long n, long w, long nx)
820 : {
821 31304 : long j, ly = n+3;
822 31304 : GEN y = cgetg(ly, t_POL);
823 31304 : y[1] = evalsigne(1);
824 384950 : for (j=2; j<ly; j++)
825 : {
826 : long k;
827 353646 : GEN a = cgetg(nx+2,t_POL);
828 353646 : a[1] = evalsigne(1) | evalvarn(w);
829 1809599 : for (k=0; k<nx; k++)
830 : {
831 1455953 : GEN xk = gel(x,k);
832 1455953 : if (typ(xk)==t_POL && varn(xk)==w)
833 1367437 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
834 : else
835 88516 : gel(a,k+2) = j==2 ? xk: gen_0;
836 : }
837 353646 : gel(y,j) = normalizepol_lg(a, nx+2);
838 : }
839 31304 : return normalizepol_lg(y,ly);
840 : }
841 :
842 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
843 : GEN
844 952 : RgXY_swap(GEN x, long n, long w)
845 : {
846 952 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
847 952 : setvarn(z, varn(x)); return z;
848 : }
849 :
850 : long
851 288 : RgXY_degreex(GEN b)
852 : {
853 288 : long deg = 0, i;
854 288 : if (!signe(b)) return -1;
855 1396 : for (i = 2; i < lg(b); ++i)
856 : {
857 1108 : GEN bi = gel(b, i);
858 1108 : if (typ(bi) == t_POL)
859 1030 : deg = maxss(deg, degpol(bi));
860 : }
861 288 : return deg;
862 : }
863 :
864 : GEN
865 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
866 :
867 : /* return (x % X^n). Shallow */
868 : GEN
869 7923949 : RgXn_red_shallow(GEN a, long n)
870 : {
871 7923949 : long i, L = n+2, l = lg(a);
872 : GEN b;
873 7923949 : if (L >= l) return a; /* deg(x) < n */
874 5735114 : b = cgetg(L, t_POL); b[1] = a[1];
875 36587483 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
876 5735116 : return normalizepol_lg(b,L);
877 : }
878 :
879 : GEN
880 483 : RgXnV_red_shallow(GEN x, long n)
881 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
882 :
883 : /* return (x * X^n). Shallow */
884 : GEN
885 175347143 : RgX_shift_shallow(GEN a, long n)
886 : {
887 175347143 : long i, l = lg(a);
888 : GEN b;
889 175347143 : if (l == 2 || !n) return a;
890 110347146 : l += n;
891 110347146 : if (n < 0)
892 : {
893 55492569 : if (l <= 2) return pol_0(varn(a));
894 54021034 : b = cgetg(l, t_POL); b[1] = a[1];
895 54021836 : a -= n;
896 165320302 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
897 : } else {
898 54854577 : b = cgetg(l, t_POL); b[1] = a[1];
899 54859253 : a -= n; n += 2;
900 119053889 : for (i=2; i<n; i++) gel(b,i) = gen_0;
901 211105585 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
902 : }
903 108881089 : return b;
904 : }
905 : /* return (x * X^n). */
906 : GEN
907 1567508 : RgX_shift(GEN a, long n)
908 : {
909 1567508 : long i, l = lg(a);
910 : GEN b;
911 1567508 : if (l == 2 || !n) return RgX_copy(a);
912 1567144 : l += n;
913 1567144 : if (n < 0)
914 : {
915 847 : if (l <= 2) return pol_0(varn(a));
916 777 : b = cgetg(l, t_POL); b[1] = a[1];
917 777 : a -= n;
918 2947 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
919 : } else {
920 1566297 : b = cgetg(l, t_POL); b[1] = a[1];
921 1566297 : a -= n; n += 2;
922 3915600 : for (i=2; i<n; i++) gel(b,i) = gen_0;
923 4319707 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
924 : }
925 1567074 : return b;
926 : }
927 :
928 : GEN
929 317037 : RgX_rotate_shallow(GEN P, long k, long p)
930 : {
931 317037 : long i, l = lgpol(P);
932 : GEN r;
933 317037 : if (signe(P)==0)
934 1365 : return pol_0(varn(P));
935 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
936 2100644 : for(i=0; i<p; i++)
937 : {
938 1784972 : long s = 2+(i+k)%p;
939 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
940 : }
941 315672 : return RgX_renormalize(r);
942 : }
943 :
944 : GEN
945 2984903 : RgX_mulXn(GEN x, long d)
946 : {
947 : pari_sp av;
948 : GEN z;
949 : long v;
950 2984903 : if (d >= 0) return RgX_shift(x, d);
951 1472286 : d = -d;
952 1472286 : v = RgX_val(x);
953 1472286 : if (v >= d) return RgX_shift(x, -d);
954 1472279 : av = avma;
955 1472279 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
956 1472279 : return gerepileupto(av, z);
957 : }
958 :
959 : long
960 588 : RgXV_maxdegree(GEN x)
961 : {
962 588 : long d = -1, i, l = lg(x);
963 4494 : for (i = 1; i < l; i++)
964 3906 : d = maxss(d, degpol(gel(x,i)));
965 588 : return d;
966 : }
967 :
968 : long
969 3647227 : RgX_val(GEN x)
970 : {
971 3647227 : long i, lx = lg(x);
972 3647227 : if (lx == 2) return LONG_MAX;
973 4577156 : for (i = 2; i < lx; i++)
974 4577100 : if (!isexactzero(gel(x,i))) break;
975 3647066 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
976 3647010 : return i - 2;
977 : }
978 : long
979 82399269 : RgX_valrem(GEN x, GEN *Z)
980 : {
981 82399269 : long v, i, lx = lg(x);
982 82399269 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
983 128190777 : for (i = 2; i < lx; i++)
984 128192187 : if (!isexactzero(gel(x,i))) break;
985 : /* possible with nonrational zeros */
986 82399384 : if (i == lx)
987 : {
988 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
989 14 : return LONG_MAX;
990 : }
991 82399370 : v = i - 2;
992 82399370 : *Z = RgX_shift_shallow(x, -v);
993 82411789 : return v;
994 : }
995 : long
996 851768 : RgX_valrem_inexact(GEN x, GEN *Z)
997 : {
998 : long v;
999 851768 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1000 868748 : for (v = 0;; v++)
1001 868748 : if (!gequal0(gel(x,2+v))) break;
1002 851761 : if (Z) *Z = RgX_shift_shallow(x, -v);
1003 851761 : return v;
1004 : }
1005 :
1006 : GEN
1007 67067 : RgXQC_red(GEN x, GEN T)
1008 409479 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1009 :
1010 : GEN
1011 1204 : RgXQV_red(GEN x, GEN T)
1012 27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1013 :
1014 : GEN
1015 12901 : RgXQM_red(GEN x, GEN T)
1016 79968 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1017 :
1018 : GEN
1019 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1020 : {
1021 322 : return RgXQM_red(RgM_mul(P, Q), T);
1022 : }
1023 :
1024 : GEN
1025 482680 : RgXQX_red(GEN P, GEN T)
1026 : {
1027 482680 : long i, l = lg(P);
1028 482680 : GEN Q = cgetg(l, t_POL);
1029 482680 : Q[1] = P[1];
1030 2546290 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1031 482680 : return normalizepol_lg(Q, l);
1032 : }
1033 :
1034 : GEN
1035 823277 : RgX_deriv(GEN x)
1036 : {
1037 823277 : long i,lx = lg(x)-1;
1038 : GEN y;
1039 :
1040 823277 : if (lx<3) return pol_0(varn(x));
1041 820190 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1042 3632210 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1043 820186 : y[1] = x[1]; return normalizepol_lg(y,i);
1044 : }
1045 :
1046 : GEN
1047 1379437 : RgX_recipspec_shallow(GEN x, long l, long n)
1048 : {
1049 : long i;
1050 1379437 : GEN z = cgetg(n+2,t_POL);
1051 1379445 : z[1] = 0; z += 2;
1052 36087822 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1053 1602696 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1054 1379445 : return normalizepol_lg(z-2,n+2);
1055 : }
1056 :
1057 : GEN
1058 639087 : RgXn_recip_shallow(GEN P, long n)
1059 : {
1060 639087 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1061 639097 : setvarn(Q, varn(P));
1062 639097 : return Q;
1063 : }
1064 :
1065 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1066 : GEN
1067 31248 : RgX_recip(GEN x)
1068 : {
1069 : long lx, i, j;
1070 31248 : GEN y = cgetg_copy(x, &lx);
1071 274029 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1072 31248 : return normalizepol_lg(y,lx);
1073 : }
1074 : /* shallow version */
1075 : GEN
1076 59367 : RgX_recip_shallow(GEN x)
1077 : {
1078 : long lx, i, j;
1079 59367 : GEN y = cgetg_copy(x, &lx);
1080 356041 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1081 59367 : return normalizepol_lg(y,lx);
1082 : }
1083 :
1084 : GEN
1085 5110314 : RgX_recip_i(GEN x)
1086 : {
1087 : long lx, i, j;
1088 5110314 : GEN y = cgetg_copy(x, &lx);
1089 25926007 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1090 5110310 : return y;
1091 : }
1092 : /*******************************************************************/
1093 : /* */
1094 : /* ADDITION / SUBTRACTION */
1095 : /* */
1096 : /*******************************************************************/
1097 : /* same variable */
1098 : GEN
1099 82376484 : RgX_add(GEN x, GEN y)
1100 : {
1101 82376484 : long i, lx = lg(x), ly = lg(y);
1102 : GEN z;
1103 82376484 : if (ly <= lx) {
1104 70051303 : z = cgetg(lx,t_POL); z[1] = x[1];
1105 279850125 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1106 123394407 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1107 70027687 : z = normalizepol_lg(z, lx);
1108 : } else {
1109 12325181 : z = cgetg(ly,t_POL); z[1] = y[1];
1110 49014732 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1111 36289164 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1112 12326508 : z = normalizepol_lg(z, ly);
1113 : }
1114 82367485 : return z;
1115 : }
1116 : GEN
1117 63717218 : RgX_sub(GEN x, GEN y)
1118 : {
1119 63717218 : long i, lx = lg(x), ly = lg(y);
1120 : GEN z;
1121 63717218 : if (ly <= lx) {
1122 29991122 : z = cgetg(lx,t_POL); z[1] = x[1];
1123 157163954 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1124 56351066 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1125 29967256 : z = normalizepol_lg(z, lx);
1126 : } else {
1127 33726096 : z = cgetg(ly,t_POL); z[1] = y[1];
1128 125453623 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1129 74656487 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1130 33695438 : z = normalizepol_lg(z, ly);
1131 : }
1132 63714387 : return z;
1133 : }
1134 : GEN
1135 7771632 : RgX_neg(GEN x)
1136 : {
1137 7771632 : long i, lx = lg(x);
1138 7771632 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
1139 49895261 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
1140 7771548 : return y;
1141 : }
1142 :
1143 : GEN
1144 22045249 : RgX_Rg_add(GEN y, GEN x)
1145 : {
1146 : GEN z;
1147 22045249 : long lz = lg(y), i;
1148 22045249 : if (lz == 2) return scalarpol(x,varn(y));
1149 18998370 : z = cgetg(lz,t_POL); z[1] = y[1];
1150 18998340 : gel(z,2) = gadd(gel(y,2),x);
1151 69154619 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1152 : /* probably useless unless lz = 3, but cannot be skipped if y is
1153 : * an inexact 0 */
1154 18998341 : return normalizepol_lg(z,lz);
1155 : }
1156 : GEN
1157 64885 : RgX_Rg_add_shallow(GEN y, GEN x)
1158 : {
1159 : GEN z;
1160 64885 : long lz = lg(y), i;
1161 64885 : if (lz == 2) return scalarpol(x,varn(y));
1162 64885 : z = cgetg(lz,t_POL); z[1] = y[1];
1163 64885 : gel(z,2) = gadd(gel(y,2),x);
1164 129930 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1165 64885 : return normalizepol_lg(z,lz);
1166 : }
1167 : GEN
1168 179650 : RgX_Rg_sub(GEN y, GEN x)
1169 : {
1170 : GEN z;
1171 179650 : long lz = lg(y), i;
1172 179650 : if (lz == 2)
1173 : { /* scalarpol(gneg(x),varn(y)) optimized */
1174 2548 : long v = varn(y);
1175 2548 : if (isrationalzero(x)) return pol_0(v);
1176 14 : z = cgetg(3,t_POL);
1177 14 : z[1] = gequal0(x)? evalvarn(v)
1178 14 : : evalvarn(v) | evalsigne(1);
1179 14 : gel(z,2) = gneg(x); return z;
1180 : }
1181 177102 : z = cgetg(lz,t_POL); z[1] = y[1];
1182 177102 : gel(z,2) = gsub(gel(y,2),x);
1183 466106 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1184 177099 : return normalizepol_lg(z,lz);
1185 : }
1186 : GEN
1187 4855913 : Rg_RgX_sub(GEN x, GEN y)
1188 : {
1189 : GEN z;
1190 4855913 : long lz = lg(y), i;
1191 4855913 : if (lz == 2) return scalarpol(x,varn(y));
1192 4854835 : z = cgetg(lz,t_POL); z[1] = y[1];
1193 4854824 : gel(z,2) = gsub(x, gel(y,2));
1194 7280525 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1195 4854828 : return normalizepol_lg(z,lz);
1196 : }
1197 : /*******************************************************************/
1198 : /* */
1199 : /* KARATSUBA MULTIPLICATION */
1200 : /* */
1201 : /*******************************************************************/
1202 : #if 0
1203 : /* to debug Karatsuba-like routines */
1204 : GEN
1205 : zx_debug_spec(GEN x, long nx)
1206 : {
1207 : GEN z = cgetg(nx+2,t_POL);
1208 : long i;
1209 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1210 : z[1] = evalsigne(1); return z;
1211 : }
1212 :
1213 : GEN
1214 : RgX_debug_spec(GEN x, long nx)
1215 : {
1216 : GEN z = cgetg(nx+2,t_POL);
1217 : long i;
1218 : for (i=0; i<nx; i++) z[i+2] = x[i];
1219 : z[1] = evalsigne(1); return z;
1220 : }
1221 : #endif
1222 :
1223 : /* generic multiplication */
1224 : GEN
1225 8885447 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1226 : {
1227 : GEN z, t;
1228 : long i;
1229 8885447 : if (nx == ny) {
1230 1529175 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1231 4933624 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1232 1529165 : return normalizepol_lg(z, nx+2);
1233 : }
1234 7356272 : if (ny < nx) {
1235 7169439 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1236 26388001 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1237 17141690 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1238 7168965 : return normalizepol_lg(z, nx+2);
1239 : } else {
1240 186833 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1241 3675482 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1242 451993 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1243 186845 : return normalizepol_lg(z, ny+2);
1244 : }
1245 : }
1246 : GEN
1247 223517 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1248 : {
1249 : GEN z, t;
1250 : long i;
1251 223517 : if (nx == ny) {
1252 12810 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1253 2183524 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1254 12810 : return normalizepol_lg(z, nx+2);
1255 : }
1256 210707 : if (ny < nx) {
1257 208929 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1258 3726848 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1259 2383505 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1260 208929 : return normalizepol_lg(z, nx+2);
1261 : } else {
1262 1778 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1263 330904 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1264 12222 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1265 1778 : return normalizepol_lg(z, ny+2);
1266 : }
1267 : }
1268 :
1269 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1270 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1271 : * t_VECSMALL, to hold this, which can be left on stack: gerepile
1272 : * will not crash on it. The returned vector itself is not a proper GEN,
1273 : * we access the coefficients as x[i], i = 0..deg(x) */
1274 : static GEN
1275 56504379 : RgXspec_kill0(GEN x, long lx)
1276 : {
1277 56504379 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1278 : long i;
1279 191968326 : for (i=0; i <lx; i++)
1280 : {
1281 135463926 : GEN c = gel(x,i);
1282 135463926 : z[i] = (long)(isrationalzero(c)? NULL: c);
1283 : }
1284 56504400 : return z;
1285 : }
1286 :
1287 : INLINE GEN
1288 94268767 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1289 : {
1290 94268767 : pari_sp av = avma;
1291 94268767 : GEN s = NULL;
1292 : long i;
1293 :
1294 356806775 : for (i=a; i<b; i++)
1295 262546962 : if (gel(y,i) && gel(x,-i))
1296 : {
1297 184524095 : GEN t = gmul(gel(y,i), gel(x,-i));
1298 184518836 : s = s? gadd(s, t): t;
1299 : }
1300 94259813 : return s? gerepileupto(av, s): gen_0;
1301 : }
1302 :
1303 : /* assume nx >= ny > 0, return x * y * t^v */
1304 : static GEN
1305 21620032 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1306 : {
1307 : long i, lz, nz;
1308 : GEN z;
1309 :
1310 21620032 : x = RgXspec_kill0(x,nx);
1311 21619999 : y = RgXspec_kill0(y,ny);
1312 21620002 : lz = nx + ny + 1; nz = lz-2;
1313 21620002 : lz += v;
1314 21620002 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1315 33359785 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1316 57283031 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1317 41348455 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1318 35663052 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1319 21619370 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1320 : }
1321 :
1322 : /* return (x * X^d) + y. Assume d > 0 */
1323 : GEN
1324 10037477 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1325 : {
1326 : GEN x, y, xd, yd, zd;
1327 : long a, lz, nx, ny;
1328 :
1329 10037477 : if (!signe(x0)) return y0;
1330 9403647 : ny = lgpol(y0);
1331 9403654 : nx = lgpol(x0);
1332 9403752 : zd = (GEN)avma;
1333 9403752 : x = x0 + 2; y = y0 + 2; a = ny-d;
1334 9403752 : if (a <= 0)
1335 : {
1336 1476739 : lz = nx+d+2;
1337 1476739 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1338 3392087 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1339 1476738 : x = zd + a;
1340 1494439 : while (zd > x) gel(--zd,0) = gen_0;
1341 : }
1342 : else
1343 : {
1344 7927013 : xd = new_chunk(d); yd = y+d;
1345 7927012 : x = RgX_addspec_shallow(x,yd, nx,a);
1346 7926976 : lz = (a>nx)? ny+2: lg(x)+d;
1347 39518164 : x += 2; while (xd > x) *--zd = *--xd;
1348 : }
1349 22910751 : while (yd > y) *--zd = *--yd;
1350 9403714 : *--zd = x0[1];
1351 9403714 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1352 : }
1353 : GEN
1354 516076 : RgX_addmulXn(GEN x0, GEN y0, long d)
1355 : {
1356 : GEN x, y, xd, yd, zd;
1357 : long a, lz, nx, ny;
1358 :
1359 516076 : if (!signe(x0)) return RgX_copy(y0);
1360 515306 : nx = lgpol(x0);
1361 515306 : ny = lgpol(y0);
1362 515306 : zd = (GEN)avma;
1363 515306 : x = x0 + 2; y = y0 + 2; a = ny-d;
1364 515306 : if (a <= 0)
1365 : {
1366 291789 : lz = nx+d+2;
1367 291789 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1368 4294362 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1369 291789 : x = zd + a;
1370 757955 : while (zd > x) gel(--zd,0) = gen_0;
1371 : }
1372 : else
1373 : {
1374 223517 : xd = new_chunk(d); yd = y+d;
1375 223517 : x = RgX_addspec(x,yd, nx,a);
1376 223517 : lz = (a>nx)? ny+2: lg(x)+d;
1377 8426296 : x += 2; while (xd > x) *--zd = *--xd;
1378 : }
1379 2623651 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1380 515306 : *--zd = x0[1];
1381 515306 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1382 : }
1383 :
1384 : /* return x * y mod t^n */
1385 : static GEN
1386 6627805 : RgXn_mul_basecase(GEN x, GEN y, long n)
1387 : {
1388 6627805 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1389 : GEN z;
1390 6627805 : if (lx < 0) return pol_0(varn(x));
1391 6627805 : if (ly < 0) return pol_0(varn(x));
1392 6627805 : z = cgetg(lz, t_POL) + 2;
1393 6627805 : x+=2; if (lx > n) lx = n;
1394 6627805 : y+=2; if (ly > n) ly = n;
1395 6627805 : z[-1] = x[-1];
1396 6627805 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1397 6627805 : x = RgXspec_kill0(x, lx);
1398 6627805 : y = RgXspec_kill0(y, ly);
1399 : /* x:y:z [i] = term of degree i */
1400 26217431 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1401 11825932 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1402 6673807 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1403 6627805 : return normalizepol_lg(z - 2, lz);
1404 : }
1405 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1406 : static GEN
1407 9413660 : RgXn_mul2(GEN f, GEN g, long n)
1408 : {
1409 9413660 : pari_sp av = avma;
1410 : GEN fe,fo, ge,go, l,h,m;
1411 : long n0, n1;
1412 9413660 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1413 6663225 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1414 35420 : n0 = n>>1; n1 = n-n0;
1415 35420 : RgX_even_odd(f, &fe, &fo);
1416 35420 : RgX_even_odd(g, &ge, &go);
1417 35420 : l = RgXn_mul2(fe,ge,n1);
1418 35420 : h = RgXn_mul2(fo,go,n0);
1419 35420 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1420 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1421 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1422 35420 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1423 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1424 35420 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1425 35420 : m = RgX_inflate(m,2);
1426 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1427 35420 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1428 35420 : h = RgX_inflate(h,2);
1429 35420 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1430 35420 : return gerepileupto(av, h);
1431 : }
1432 : /* (f*g) \/ x^n */
1433 : static GEN
1434 1589258 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1435 : {
1436 1589258 : long d = degpol(f)+degpol(g) + 1 - n;
1437 : GEN h;
1438 1589258 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1439 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1440 : RgX_recip_i(g), d));
1441 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1442 : }
1443 :
1444 : /* (f*g) \/ x^n */
1445 : static GEN
1446 0 : RgX_sqrhigh_i2(GEN f, long n)
1447 : {
1448 0 : long d = 2*degpol(f)+ 1 - n;
1449 : GEN h;
1450 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1451 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1452 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1453 : }
1454 :
1455 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1456 : * b+2 were sent instead. na, nb = number of terms of a, b.
1457 : * Only c, c0, c1, c2 are genuine GEN.
1458 : */
1459 : GEN
1460 22349262 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1461 : {
1462 : GEN a0, c, c0;
1463 22349262 : long n0, n0a, i, v = 0;
1464 : pari_sp av;
1465 :
1466 28732803 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1467 28239649 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1468 22349220 : if (na < nb) swapspec(a,b, na,nb);
1469 22349220 : if (!nb) return pol_0(0);
1470 :
1471 22100628 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1472 480606 : RgX_shift_inplace_init(v);
1473 480617 : i = (na>>1); n0 = na-i; na = i;
1474 480617 : av = avma; a0 = a+n0; n0a = n0;
1475 1336133 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1476 :
1477 480617 : if (nb > n0)
1478 : {
1479 : GEN b0,c1,c2;
1480 : long n0b;
1481 :
1482 479218 : nb -= n0; b0 = b+n0; n0b = n0;
1483 1447197 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1484 479218 : c = RgX_mulspec(a,b,n0a,n0b);
1485 479218 : c0 = RgX_mulspec(a0,b0, na,nb);
1486 :
1487 479218 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1488 479218 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1489 :
1490 479218 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1491 479218 : c2 = RgX_sub(c1, RgX_add(c0,c));
1492 479218 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1493 : }
1494 : else
1495 : {
1496 1399 : c = RgX_mulspec(a,b,n0a,nb);
1497 1399 : c0 = RgX_mulspec(a0,b,na,nb);
1498 : }
1499 480617 : c0 = RgX_addmulXn(c0,c,n0);
1500 480617 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1501 : }
1502 :
1503 : INLINE GEN
1504 50522 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1505 : {
1506 50522 : pari_sp av = avma;
1507 50522 : GEN s = NULL;
1508 50522 : long j, l = (i+1)>>1;
1509 139734 : for (j=a; j<l; j++)
1510 : {
1511 89212 : GEN xj = gel(x,j), xx = gel(x,i-j);
1512 89212 : if (xj && xx)
1513 : {
1514 81568 : GEN t = gmul(xj, xx);
1515 81568 : s = s? gadd(s, t): t;
1516 : }
1517 : }
1518 50522 : if (s) s = gshift(s,1);
1519 50522 : if ((i&1) == 0)
1520 : {
1521 29715 : GEN t = gel(x, i>>1);
1522 29715 : if (t) {
1523 27321 : t = gsqr(t);
1524 27321 : s = s? gadd(s, t): t;
1525 : }
1526 : }
1527 50522 : return s? gerepileupto(av,s): gen_0;
1528 : }
1529 : static GEN
1530 8908 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1531 : {
1532 : long i, lz, nz;
1533 : GEN z;
1534 :
1535 8908 : if (!nx) return pol_0(0);
1536 8908 : x = RgXspec_kill0(x,nx);
1537 8908 : lz = (nx << 1) + 1, nz = lz-2;
1538 8908 : lz += v;
1539 8908 : z = cgetg(lz,t_POL) + 2;
1540 13080 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1541 38623 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1542 29715 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1543 8908 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1544 : }
1545 : /* return x^2 mod t^n */
1546 : static GEN
1547 0 : RgXn_sqr_basecase(GEN x, long n)
1548 : {
1549 0 : long i, lz = n+2, lx = lgpol(x);
1550 : GEN z;
1551 0 : if (lx < 0) return pol_0(varn(x));
1552 0 : z = cgetg(lz, t_POL);
1553 0 : z[1] = x[1];
1554 0 : x+=2; if (lx > n) lx = n;
1555 0 : x = RgXspec_kill0(x,lx);
1556 0 : z+=2;/* x:z [i] = term of degree i */
1557 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1558 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1559 0 : z -= 2; return normalizepol_lg(z, lz);
1560 : }
1561 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1562 : static GEN
1563 315 : RgXn_sqr2(GEN f, long n)
1564 : {
1565 315 : pari_sp av = avma;
1566 : GEN fe,fo, l,h,m;
1567 : long n0, n1;
1568 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1569 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1570 0 : n0 = n>>1; n1 = n-n0;
1571 0 : RgX_even_odd(f, &fe, &fo);
1572 0 : l = RgXn_sqr(fe,n1);
1573 0 : h = RgXn_sqr(fo,n0);
1574 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1575 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1576 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1577 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1578 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1579 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1580 0 : m = RgX_inflate(m,2);
1581 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1582 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1583 0 : h = RgX_inflate(h,2);
1584 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1585 0 : return gerepileupto(av, h);
1586 : }
1587 : GEN
1588 8947 : RgX_sqrspec(GEN a, long na)
1589 : {
1590 : GEN a0, c, c0, c1;
1591 8947 : long n0, n0a, i, v = 0;
1592 : pari_sp av;
1593 :
1594 11033 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1595 8947 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1596 39 : RgX_shift_inplace_init(v);
1597 39 : i = (na>>1); n0 = na-i; na = i;
1598 39 : av = avma; a0 = a+n0; n0a = n0;
1599 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1600 :
1601 39 : c = RgX_sqrspec(a,n0a);
1602 39 : c0 = RgX_sqrspec(a0,na);
1603 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1604 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1605 39 : c0 = RgX_addmulXn(c0,c,n0);
1606 39 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1607 : }
1608 :
1609 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1610 : GEN
1611 1711152 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1612 : {
1613 1711152 : GEN z = RgX_mul(A, B);
1614 1711136 : if (a < b)
1615 10445 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1616 1700691 : else if (a > b)
1617 1098945 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1618 : else
1619 601746 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1620 1711146 : return z;
1621 : }
1622 :
1623 : GEN
1624 20908822 : RgX_mul_i(GEN x, GEN y)
1625 : {
1626 20908822 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1627 20908360 : setvarn(z, varn(x)); return z;
1628 : }
1629 :
1630 : GEN
1631 8869 : RgX_sqr_i(GEN x)
1632 : {
1633 8869 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1634 8869 : setvarn(z,varn(x)); return z;
1635 : }
1636 :
1637 : /*******************************************************************/
1638 : /* */
1639 : /* DIVISION */
1640 : /* */
1641 : /*******************************************************************/
1642 : GEN
1643 7161801 : RgX_Rg_divexact(GEN x, GEN y) {
1644 7161801 : long i, lx = lg(x);
1645 : GEN z;
1646 7161801 : if (lx == 2) return gcopy(x);
1647 7111040 : switch(typ(y))
1648 : {
1649 6951319 : case t_INT:
1650 6951319 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1651 6576777 : break;
1652 5215 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1653 : }
1654 6731283 : z = cgetg(lx, t_POL); z[1] = x[1];
1655 34373445 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1656 6730008 : return z;
1657 : }
1658 : GEN
1659 30169370 : RgX_Rg_div(GEN x, GEN y) {
1660 30169370 : long i, lx = lg(x);
1661 : GEN z;
1662 30169370 : if (lx == 2) return gcopy(x);
1663 29909312 : switch(typ(y))
1664 : {
1665 20715694 : case t_INT:
1666 20715694 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1667 4116186 : break;
1668 13125 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1669 : }
1670 13296679 : z = cgetg(lx, t_POL); z[1] = x[1];
1671 49573978 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1672 13296200 : return normalizepol_lg(z, lx);
1673 : }
1674 : GEN
1675 39018 : RgX_normalize(GEN x)
1676 : {
1677 39018 : GEN z, d = NULL;
1678 39018 : long i, n = lg(x)-1;
1679 39018 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1680 39018 : if (i == 1) return pol_0(varn(x));
1681 39018 : if (i == n && isint1(d)) return x;
1682 17584 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1683 31570 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1684 17584 : gel(z,n) = Rg_get_1(d); return z;
1685 : }
1686 : GEN
1687 10423 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1688 : GEN
1689 250634 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1690 : {
1691 250634 : long l = lg(a), i;
1692 : GEN a0, z0, z;
1693 :
1694 250634 : if (l <= 3)
1695 : {
1696 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1697 0 : return pol_0(varn(a));
1698 : }
1699 250634 : z = cgetg(l-1, t_POL);
1700 250635 : z[1] = a[1];
1701 250635 : a0 = a + l-1;
1702 250635 : z0 = z + l-2; *z0 = *a0--;
1703 2954279 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1704 : {
1705 2703642 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1706 2703644 : gel(z0,0) = t;
1707 : }
1708 250637 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1709 250637 : return z;
1710 : }
1711 : /* Polynomial division x / y:
1712 : * if pr = ONLY_REM return remainder, otherwise return quotient
1713 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1714 : * if pr != NULL set *pr to remainder, as the last object on stack */
1715 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1716 : static GEN
1717 23610744 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1718 : {
1719 : pari_sp avy, av, av1;
1720 : long dx,dy,dz,i,j,sx,lr;
1721 : GEN z,p1,p2,rem,y_lead,mod,p;
1722 : GEN (*f)(GEN,GEN);
1723 :
1724 23610744 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1725 :
1726 23610744 : dy = degpol(y);
1727 23610695 : y_lead = gel(y,dy+2);
1728 23610695 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1729 : {
1730 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1731 0 : for (dy--; dy>=0; dy--)
1732 : {
1733 0 : y_lead = gel(y,dy+2);
1734 0 : if (!gequal0(y_lead)) break;
1735 : }
1736 : }
1737 23610696 : if (!dy) /* y is constant */
1738 : {
1739 6843 : if (pr == ONLY_REM) return pol_0(varn(x));
1740 6836 : z = RgX_Rg_div(x, y_lead);
1741 6836 : if (pr == ONLY_DIVIDES) return z;
1742 2622 : if (pr) *pr = pol_0(varn(x));
1743 2622 : return z;
1744 : }
1745 23603853 : dx = degpol(x);
1746 23603788 : if (dx < dy)
1747 : {
1748 3858635 : if (pr == ONLY_REM) return RgX_copy(x);
1749 374785 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1750 374764 : z = pol_0(varn(x));
1751 374764 : if (pr) *pr = RgX_copy(x);
1752 374764 : return z;
1753 : }
1754 :
1755 : /* x,y in R[X], y non constant */
1756 19745153 : av = avma;
1757 19745153 : p = NULL;
1758 19745153 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1759 : {
1760 238784 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1761 238784 : if (!z) return gc_NULL(av);
1762 238784 : z = FpX_to_mod(z, p);
1763 238784 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1764 122528 : return gerepileupto(av, z);
1765 116256 : *pr = FpX_to_mod(*pr, p);
1766 116256 : return gc_all(av, 2, &z, pr);
1767 : }
1768 19506527 : switch(typ(y_lead))
1769 : {
1770 0 : case t_REAL:
1771 0 : y_lead = ginv(y_lead);
1772 0 : f = gmul; mod = NULL;
1773 0 : break;
1774 2074 : case t_INTMOD:
1775 2074 : case t_POLMOD: y_lead = ginv(y_lead);
1776 2074 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1777 2074 : break;
1778 19504453 : default: if (gequal1(y_lead)) y_lead = NULL;
1779 19504438 : f = gdiv; mod = NULL;
1780 : }
1781 :
1782 19506512 : if (y_lead == NULL)
1783 17492016 : p2 = gel(x,dx+2);
1784 : else {
1785 : for(;;) {
1786 2014496 : p2 = f(gel(x,dx+2),y_lead);
1787 2014497 : p2 = simplify_shallow(p2);
1788 2014497 : if (!isexactzero(p2) || (--dx < 0)) break;
1789 : }
1790 2014497 : if (dx < dy) /* leading coeff of x was in fact zero */
1791 : {
1792 0 : if (pr == ONLY_DIVIDES) {
1793 0 : set_avma(av);
1794 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1795 : }
1796 0 : if (pr == ONLY_REM)
1797 : {
1798 0 : if (dx < 0)
1799 0 : return gerepilecopy(av, scalarpol(p2, varn(x)));
1800 : else
1801 : {
1802 : GEN t;
1803 0 : set_avma(av);
1804 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1805 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1806 0 : return t;
1807 : }
1808 : }
1809 0 : if (pr) /* cf ONLY_REM above */
1810 : {
1811 0 : if (dx < 0)
1812 : {
1813 0 : p2 = gclone(p2);
1814 0 : set_avma(av);
1815 0 : z = pol_0(varn(x));
1816 0 : x = scalarpol(p2, varn(x));
1817 0 : gunclone(p2);
1818 : }
1819 : else
1820 : {
1821 : GEN t;
1822 0 : set_avma(av);
1823 0 : z = pol_0(varn(x));
1824 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1825 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1826 0 : x = t;
1827 : }
1828 0 : *pr = x;
1829 : }
1830 : else
1831 : {
1832 0 : set_avma(av);
1833 0 : z = pol_0(varn(x));
1834 : }
1835 0 : return z;
1836 : }
1837 : }
1838 : /* dx >= dy */
1839 19506513 : avy = avma;
1840 19506513 : dz = dx-dy;
1841 19506513 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1842 19506480 : x += 2;
1843 19506480 : z += 2;
1844 19506480 : y += 2;
1845 19506480 : gel(z,dz) = gcopy(p2);
1846 :
1847 54027443 : for (i=dx-1; i>=dy; i--)
1848 : {
1849 34521150 : av1=avma; p1=gel(x,i);
1850 1140354277 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1851 34516305 : if (y_lead) p1 = simplify(f(p1,y_lead));
1852 :
1853 34516305 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1854 : else
1855 24526641 : p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1856 34520523 : gel(z,i-dy) = p1;
1857 : }
1858 19506293 : if (!pr) return gerepileupto(av,z-2);
1859 :
1860 12531526 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1861 12894805 : for (sx=0; ; i--)
1862 : {
1863 12894805 : p1 = gel(x,i);
1864 : /* we always enter this loop at least once */
1865 32376431 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1866 12893728 : if (mod && avma==av1) p1 = gmul(p1,mod);
1867 12893728 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1868 3585399 : if (!isexactzero(p1)) break;
1869 3436185 : if (!i) break;
1870 363259 : set_avma(av1);
1871 : }
1872 12530801 : if (pr == ONLY_DIVIDES)
1873 : {
1874 7917 : if (sx) return gc_NULL(av);
1875 7896 : set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1876 : }
1877 12522884 : lr=i+3; rem -= lr;
1878 12522884 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1879 12479265 : else p1 = gerepileupto((pari_sp)rem,p1);
1880 12523582 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1881 12523582 : rem[1] = z[-1];
1882 12523582 : rem += 2;
1883 12523582 : gel(rem,i) = p1;
1884 17210189 : for (i--; i>=0; i--)
1885 : {
1886 4686627 : av1=avma; p1 = gel(x,i);
1887 13306147 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1888 4686181 : if (mod && avma==av1) p1 = gmul(p1,mod);
1889 4686375 : gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1890 : }
1891 12523562 : rem -= 2;
1892 12523562 : if (!sx) (void)normalizepol_lg(rem, lr);
1893 12523793 : if (pr == ONLY_REM) return gerepileupto(av,rem);
1894 6758628 : z -= 2;
1895 : {
1896 6758628 : GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1897 6758628 : gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1898 : }
1899 : }
1900 :
1901 : GEN
1902 14361746 : RgX_divrem(GEN x, GEN y, GEN *pr)
1903 : {
1904 14361746 : if (pr == ONLY_REM) return RgX_rem(x, y);
1905 14361746 : return RgX_divrem_i(x, y, pr);
1906 : }
1907 :
1908 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1909 : GEN
1910 143338 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1911 : {
1912 : long vx, dx, dy, dz, i, j, sx, lr;
1913 : pari_sp av0, av, tetpil;
1914 : GEN z,p1,rem,lead;
1915 :
1916 143338 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1917 143338 : vx = varn(x);
1918 143338 : dx = degpol(x);
1919 143338 : dy = degpol(y);
1920 143338 : if (dx < dy)
1921 : {
1922 31466 : if (pr)
1923 : {
1924 31438 : av0 = avma; x = RgXQX_red(x, T);
1925 31438 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1926 31431 : if (pr == ONLY_REM) return x;
1927 0 : *pr = x;
1928 : }
1929 28 : return pol_0(vx);
1930 : }
1931 111872 : lead = leading_coeff(y);
1932 111872 : if (!dy) /* y is constant */
1933 : {
1934 602 : if (pr && pr != ONLY_DIVIDES)
1935 : {
1936 0 : if (pr == ONLY_REM) return pol_0(vx);
1937 0 : *pr = pol_0(vx);
1938 : }
1939 602 : if (gequal1(lead)) return RgX_copy(x);
1940 0 : av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
1941 0 : return gerepile(av0,tetpil,RgXQX_red(x,T));
1942 : }
1943 111270 : av0 = avma; dz = dx-dy;
1944 111270 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
1945 111270 : set_avma(av0);
1946 111270 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1947 111270 : x += 2; y += 2; z += 2;
1948 :
1949 111270 : p1 = gel(x,dx); av = avma;
1950 111270 : gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
1951 571003 : for (i=dx-1; i>=dy; i--)
1952 : {
1953 459733 : av=avma; p1=gel(x,i);
1954 2395195 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1955 459693 : if (lead) p1 = gmul(grem(p1, T), lead);
1956 459698 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
1957 : }
1958 111270 : if (!pr) { guncloneNULL(lead); return z-2; }
1959 :
1960 109996 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
1961 185634 : for (sx=0; ; i--)
1962 : {
1963 185634 : p1 = gel(x,i);
1964 642520 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1965 185631 : tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
1966 102833 : if (!i) break;
1967 75638 : set_avma(av);
1968 : }
1969 109996 : if (pr == ONLY_DIVIDES)
1970 : {
1971 27190 : guncloneNULL(lead);
1972 27190 : if (sx) return gc_NULL(av0);
1973 24619 : return gc_const((pari_sp)rem, z-2);
1974 : }
1975 82806 : lr=i+3; rem -= lr;
1976 82806 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1977 82806 : rem[1] = z[-1];
1978 82806 : p1 = gerepile((pari_sp)rem,tetpil,p1);
1979 82806 : rem += 2; gel(rem,i) = p1;
1980 179107 : for (i--; i>=0; i--)
1981 : {
1982 96301 : av=avma; p1 = gel(x,i);
1983 322820 : for (j=0; j<=i && j<=dz; j++)
1984 226519 : p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1985 96301 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
1986 : }
1987 82806 : rem -= 2;
1988 82806 : guncloneNULL(lead);
1989 82806 : if (!sx) (void)normalizepol_lg(rem, lr);
1990 82806 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
1991 161 : *pr = rem; return z-2;
1992 : }
1993 :
1994 : /*******************************************************************/
1995 : /* */
1996 : /* PSEUDO-DIVISION */
1997 : /* */
1998 : /*******************************************************************/
1999 : INLINE GEN
2000 4039539 : rem(GEN c, GEN T)
2001 : {
2002 4039539 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2003 4039539 : return c;
2004 : }
2005 :
2006 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2007 : int
2008 17282 : ZXQX_dvd(GEN x, GEN y, GEN T)
2009 : {
2010 : long dx, dy, i, T_ismonic;
2011 17282 : pari_sp av = avma, av2;
2012 : GEN y_lead;
2013 :
2014 17282 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2015 17282 : dy = degpol(y); y_lead = gel(y,dy+2);
2016 17282 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2017 : /* if monic, no point in using pseudo-division */
2018 17282 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2019 14741 : T_ismonic = gequal1(leading_coeff(T));
2020 14741 : dx = degpol(x);
2021 14741 : if (dx < dy) return !signe(x);
2022 14741 : (void)new_chunk(2);
2023 14741 : x = RgX_recip_i(x)+2;
2024 14741 : y = RgX_recip_i(y)+2;
2025 : /* pay attention to sparse divisors */
2026 29983 : for (i = 1; i <= dy; i++)
2027 15242 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2028 14741 : av2 = avma;
2029 : for (;;)
2030 72872 : {
2031 87613 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2032 87613 : x0 = gneg(x0);
2033 87613 : m = gcdii(cx, y0);
2034 87613 : if (!equali1(m))
2035 : {
2036 85118 : x0 = gdiv(x0, m);
2037 85118 : y0 = diviiexact(y0, m);
2038 85118 : if (equali1(y0)) y0 = NULL;
2039 : }
2040 178999 : for (i=1; i<=dy; i++)
2041 : {
2042 91386 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2043 91386 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2044 91386 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2045 91386 : gel(x,i) = c;
2046 : }
2047 687617 : for ( ; i<=dx; i++)
2048 : {
2049 600004 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2050 600004 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2051 600004 : gel(x,i) = c;
2052 : }
2053 102110 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2054 87613 : if (dx < dy) break;
2055 72872 : if (gc_needed(av2,1))
2056 : {
2057 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2058 28 : gerepilecoeffs(av2,x,dx+1);
2059 : }
2060 : }
2061 14741 : return gc_bool(av, dx < 0);
2062 : }
2063 :
2064 : /* T either NULL or a t_POL. */
2065 : GEN
2066 629959 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2067 : {
2068 629959 : long vx = varn(x), dx, dy, dz, i, lx, p;
2069 629959 : pari_sp av = avma, av2;
2070 : GEN y_lead;
2071 :
2072 629959 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2073 629959 : dy = degpol(y); y_lead = gel(y,dy+2);
2074 : /* if monic, no point in using pseudo-division */
2075 629960 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2076 580584 : dx = degpol(x);
2077 580584 : if (dx < dy) return RgX_copy(x);
2078 580584 : (void)new_chunk(2);
2079 580583 : x = RgX_recip_i(x)+2;
2080 580584 : y = RgX_recip_i(y)+2;
2081 : /* pay attention to sparse divisors */
2082 1956001 : for (i = 1; i <= dy; i++)
2083 1375417 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2084 580584 : dz = dx-dy; p = dz+1;
2085 580584 : av2 = avma;
2086 : for (;;)
2087 : {
2088 1247240 : gel(x,0) = gneg(gel(x,0)); p--;
2089 4193328 : for (i=1; i<=dy; i++)
2090 : {
2091 2946096 : GEN c = gmul(y_lead, gel(x,i));
2092 2946007 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2093 2946071 : gel(x,i) = rem(c, T);
2094 : }
2095 2179077 : for ( ; i<=dx; i++)
2096 : {
2097 931839 : GEN c = gmul(y_lead, gel(x,i));
2098 931818 : gel(x,i) = rem(c, T);
2099 : }
2100 1293441 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2101 1247237 : if (dx < dy) break;
2102 666656 : if (gc_needed(av2,1))
2103 : {
2104 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2105 0 : gerepilecoeffs(av2,x,dx+1);
2106 : }
2107 : }
2108 580581 : if (dx < 0) return pol_0(vx);
2109 579398 : lx = dx+3; x -= 2;
2110 579398 : x[0] = evaltyp(t_POL) | _evallg(lx);
2111 579398 : x[1] = evalsigne(1) | evalvarn(vx);
2112 579398 : x = RgX_recip_i(x);
2113 579398 : if (p)
2114 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2115 28612 : GEN t = y_lead;
2116 28612 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2117 0 : t = RgXQ_powu(t, p, T);
2118 : else
2119 28612 : t = gpowgs(t, p);
2120 88535 : for (i=2; i<lx; i++)
2121 : {
2122 59923 : GEN c = gmul(gel(x,i), t);
2123 59923 : gel(x,i) = rem(c,T);
2124 : }
2125 28612 : if (!T) return gerepileupto(av, x);
2126 : }
2127 550786 : return gerepilecopy(av, x);
2128 : }
2129 :
2130 : GEN
2131 629959 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2132 :
2133 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2134 : GEN
2135 12084 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2136 : {
2137 12084 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2138 12084 : pari_sp av = avma, av2;
2139 : GEN z, r, ypow, y_lead;
2140 :
2141 12084 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2142 12084 : dy = degpol(y); y_lead = gel(y,dy+2);
2143 12084 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2144 7589 : dx = degpol(x);
2145 7589 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2146 7589 : if (dx == dy)
2147 : {
2148 98 : GEN x_lead = gel(x,lg(x)-1);
2149 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2150 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2151 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2152 98 : *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
2153 : }
2154 7491 : (void)new_chunk(2);
2155 7491 : x = RgX_recip_i(x)+2;
2156 7491 : y = RgX_recip_i(y)+2;
2157 : /* pay attention to sparse divisors */
2158 39049 : for (i = 1; i <= dy; i++)
2159 31558 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2160 7491 : dz = dx-dy; p = dz+1;
2161 7491 : lz = dz+3;
2162 7491 : z = cgetg(lz, t_POL);
2163 7491 : z[1] = evalsigne(1) | evalvarn(vx);
2164 28745 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2165 7491 : ypow = new_chunk(dz+1);
2166 7491 : gel(ypow,0) = gen_1;
2167 7491 : gel(ypow,1) = y_lead;
2168 13763 : for (i=2; i<=dz; i++)
2169 : {
2170 6272 : GEN c = gmul(gel(ypow,i-1), y_lead);
2171 6272 : gel(ypow,i) = rem(c,T);
2172 : }
2173 7491 : av2 = avma;
2174 7491 : for (iz=2;;)
2175 : {
2176 15682 : p--;
2177 15682 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2178 70034 : for (i=1; i<=dy; i++)
2179 : {
2180 54352 : GEN c = gmul(y_lead, gel(x,i));
2181 54352 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2182 54352 : gel(x,i) = rem(c, T);
2183 : }
2184 41184 : for ( ; i<=dx; i++)
2185 : {
2186 25502 : GEN c = gmul(y_lead, gel(x,i));
2187 25502 : gel(x,i) = rem(c,T);
2188 : }
2189 15682 : x++; dx--;
2190 21254 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2191 15682 : if (dx < dy) break;
2192 8191 : if (gc_needed(av2,1))
2193 : {
2194 0 : GEN X = x-2;
2195 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2196 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2197 0 : gerepileall(av2,2, &X, &z); x = X+2;
2198 : }
2199 : }
2200 14036 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2201 7491 : if (dx < 0)
2202 182 : x = pol_0(vx);
2203 : else
2204 : {
2205 7309 : lx = dx+3; x -= 2;
2206 7309 : x[0] = evaltyp(t_POL) | _evallg(lx);
2207 7309 : x[1] = evalsigne(1) | evalvarn(vx);
2208 7309 : x = RgX_recip_i(x);
2209 : }
2210 7491 : z = RgX_recip_i(z);
2211 7491 : r = x;
2212 7491 : if (p)
2213 : {
2214 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2215 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2216 : }
2217 7491 : *ptr = r; return gc_all(av, 2, &z, ptr);
2218 : }
2219 : GEN
2220 11846 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2221 11846 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2222 :
2223 : GEN
2224 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2225 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2226 : GEN
2227 523769181 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2228 : GEN
2229 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2230 : GEN
2231 26376 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2232 : GEN
2233 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2234 : GEN
2235 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2236 :
2237 : GEN
2238 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2239 :
2240 : GEN
2241 0 : RgXQX_powers(GEN P, long n, GEN T)
2242 : {
2243 0 : GEN v = cgetg(n+2, t_VEC);
2244 : long i;
2245 0 : gel(v, 1) = pol_1(varn(T));
2246 0 : if (n==0) return v;
2247 0 : gel(v, 2) = gcopy(P);
2248 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2249 0 : return v;
2250 : }
2251 :
2252 : static GEN
2253 613863 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2254 : static GEN
2255 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2256 : static GEN
2257 67352 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2258 : static GEN
2259 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2260 : static GEN
2261 274852 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2262 : static GEN
2263 1023155 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2264 : static GEN
2265 1001252 : _one(void *data) { return pol_1(varn((GEN)data)); }
2266 : static GEN
2267 973 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2268 : static GEN
2269 713483 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2270 :
2271 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2272 : _mul, _sqr, _one, _zero };
2273 :
2274 : GEN
2275 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2276 : {
2277 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2278 : }
2279 :
2280 : GEN
2281 409296 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2282 : {
2283 409296 : int use_sqr = 2*degpol(x) >= degpol(T);
2284 409296 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2285 : }
2286 :
2287 : /* mod X^n */
2288 : struct modXn {
2289 : long v; /* varn(X) */
2290 : long n;
2291 : } ;
2292 : static GEN
2293 11893 : _sqrXn(void *data, GEN x) {
2294 11893 : struct modXn *S = (struct modXn*)data;
2295 11893 : return RgXn_sqr(x, S->n);
2296 : }
2297 : static GEN
2298 4528 : _mulXn(void *data, GEN x, GEN y) {
2299 4528 : struct modXn *S = (struct modXn*)data;
2300 4528 : return RgXn_mul(x,y, S->n);
2301 : }
2302 : static GEN
2303 1939 : _oneXn(void *data) {
2304 1939 : struct modXn *S = (struct modXn*)data;
2305 1939 : return pol_1(S->v);
2306 : }
2307 : static GEN
2308 0 : _zeroXn(void *data) {
2309 0 : struct modXn *S = (struct modXn*)data;
2310 0 : return pol_0(S->v);
2311 : }
2312 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2313 : _oneXn, _zeroXn };
2314 :
2315 : GEN
2316 483 : RgXn_powers(GEN x, long m, long n)
2317 : {
2318 483 : long d = degpol(x);
2319 483 : int use_sqr = (d<<1) >= n;
2320 : struct modXn S;
2321 483 : S.v = varn(x); S.n = n;
2322 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2323 : }
2324 :
2325 : GEN
2326 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2327 : {
2328 : struct modXn S;
2329 : long v;
2330 2286 : if (n == 0) return x;
2331 2286 : v = RgX_valrem(x, &x);
2332 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2333 2265 : S.v = varn(x); S.n = n;
2334 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2335 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2336 2265 : return x;
2337 : }
2338 : GEN
2339 0 : RgXn_powu(GEN x, ulong m, long n)
2340 : {
2341 : pari_sp av;
2342 0 : if (n == 0) return gcopy(x);
2343 0 : av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
2344 : }
2345 :
2346 : GEN
2347 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2348 : {
2349 : struct modXn S;
2350 966 : S.v = varn(gel(x,2)); S.n = n;
2351 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2352 : }
2353 :
2354 : GEN
2355 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2356 : {
2357 0 : int use_sqr = 2*degpol(x) >= n;
2358 : struct modXn S;
2359 0 : S.v = varn(x); S.n = n;
2360 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2361 : }
2362 :
2363 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2364 : GEN
2365 5145 : RgXn_eval(GEN Q, GEN x, long n)
2366 : {
2367 5145 : long d = degpol(x);
2368 : int use_sqr;
2369 : struct modXn S;
2370 5145 : if (d == 1 && isrationalzero(gel(x,2)))
2371 : {
2372 5138 : GEN y = RgX_unscale(Q, gel(x,3));
2373 5138 : setvarn(y, varn(x)); return y;
2374 : }
2375 7 : S.v = varn(x);
2376 7 : S.n = n;
2377 7 : use_sqr = (d<<1) >= n;
2378 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2379 : }
2380 :
2381 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2382 : static GEN
2383 2614484 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2384 : {
2385 2614484 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2386 2614484 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2387 : }
2388 :
2389 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2390 : static GEN
2391 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2392 : {
2393 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2394 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2395 : }
2396 :
2397 : static GEN
2398 2188899 : RgXn_div_gen(GEN g, GEN f, long e)
2399 : {
2400 : pari_sp av;
2401 : ulong mask;
2402 : GEN W, a;
2403 2188899 : long v = varn(f), n = 1;
2404 :
2405 2188899 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2406 2188892 : a = ginv(gel(f,2));
2407 2188897 : if (e == 1 && !g) return scalarpol(a, v);
2408 2184025 : else if (e == 2 && !g)
2409 : {
2410 : GEN b;
2411 836929 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2412 246249 : b = gneg(b);
2413 246247 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2414 246247 : return deg1pol(b, a, v);
2415 : }
2416 1347096 : av = avma;
2417 1347096 : W = scalarpol_shallow(a,v);
2418 1347096 : mask = quadratic_prec_mask(e);
2419 3955623 : while (mask > 1)
2420 : {
2421 : GEN u, fr;
2422 2608527 : long n2 = n;
2423 2608527 : n<<=1; if (mask & 1) n--;
2424 2608527 : mask >>= 1;
2425 2608527 : fr = RgXn_red_shallow(f, n);
2426 2608527 : if (mask>1 || !g)
2427 : {
2428 1324547 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2429 1324547 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2430 : }
2431 : else
2432 : {
2433 1283980 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2434 1283980 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2435 1283980 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2436 : }
2437 2608527 : if (gc_needed(av,2))
2438 : {
2439 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2440 0 : W = gerepileupto(av, W);
2441 : }
2442 : }
2443 1347096 : return W;
2444 : }
2445 :
2446 : static GEN
2447 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2448 : {
2449 : GEN r;
2450 119 : if (lgefint(p) == 3)
2451 : {
2452 119 : ulong pp = uel(p, 2);
2453 119 : if (pp == 2)
2454 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2455 : else
2456 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2457 : }
2458 : else
2459 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2460 119 : return FpX_to_mod(r, p);
2461 : }
2462 :
2463 : static GEN
2464 7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2465 : {
2466 7 : GEN r, T = RgX_to_FpX(pol, p);
2467 7 : if (signe(T) == 0) pari_err_OP("/", x, y);
2468 7 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2469 7 : return FpXQX_to_mod(r, T, p);
2470 : }
2471 :
2472 : static GEN
2473 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2474 : {
2475 : GEN r;
2476 91 : if (lgefint(p) == 3)
2477 : {
2478 91 : ulong pp = uel(p, 2);
2479 91 : if (pp == 2)
2480 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2481 : else
2482 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2483 : }
2484 : else
2485 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2486 91 : return FpX_to_mod(r, p);
2487 : }
2488 :
2489 : static GEN
2490 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2491 : {
2492 0 : GEN r, T = RgX_to_FpX(pol, p);
2493 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2494 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2495 0 : return FpXQX_to_mod(r, T, p);
2496 : }
2497 :
2498 : #define code(t1,t2) ((t1 << 6) | t2)
2499 :
2500 : static GEN
2501 905018 : RgXn_inv_fast(GEN x, long e)
2502 : {
2503 : GEN p, pol;
2504 : long pa;
2505 905018 : long t = RgX_type(x,&p,&pol,&pa);
2506 905019 : switch(t)
2507 : {
2508 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2509 0 : case code(t_POLMOD, t_INTMOD):
2510 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2511 904928 : default: return NULL;
2512 : }
2513 : }
2514 :
2515 : static GEN
2516 1284106 : RgXn_div_fast(GEN x, GEN y, long e)
2517 : {
2518 : GEN p, pol;
2519 : long pa;
2520 1284106 : long t = RgX_type2(x,y,&p,&pol,&pa);
2521 1284106 : switch(t)
2522 : {
2523 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2524 7 : case code(t_POLMOD, t_INTMOD):
2525 7 : return RgXn_div_FpXQX(x, y, e, pol, p);
2526 1283980 : default: return NULL;
2527 : }
2528 : }
2529 : #undef code
2530 :
2531 : GEN
2532 1284106 : RgXn_div_i(GEN g, GEN f, long e)
2533 : {
2534 1284106 : GEN h = RgXn_div_fast(g, f, e);
2535 1284106 : if (h) return h;
2536 1283980 : return RgXn_div_gen(g, f, e);
2537 : }
2538 :
2539 : GEN
2540 584443 : RgXn_div(GEN g, GEN f, long e)
2541 : {
2542 584443 : pari_sp av = avma;
2543 584443 : return gerepileupto(av, RgXn_div_i(g, f, e));
2544 : }
2545 :
2546 : GEN
2547 905019 : RgXn_inv_i(GEN f, long e)
2548 : {
2549 905019 : GEN h = RgXn_inv_fast(f, e);
2550 905019 : if (h) return h;
2551 904928 : return RgXn_div_gen(NULL, f, e);
2552 : }
2553 :
2554 : GEN
2555 807451 : RgXn_inv(GEN f, long e)
2556 : {
2557 807451 : pari_sp av = avma;
2558 807451 : return gerepileupto(av, RgXn_inv_i(f, e));
2559 : }
2560 :
2561 : /* Compute intformal(x^n*S)/x^(n+1) */
2562 : static GEN
2563 56373 : RgX_integXn(GEN x, long n)
2564 : {
2565 56373 : long i, lx = lg(x);
2566 : GEN y;
2567 56373 : if (lx == 2) return RgX_copy(x);
2568 52782 : y = cgetg(lx, t_POL); y[1] = x[1];
2569 110759 : for (i=2; i<lx; i++)
2570 57977 : gel(y,i) = gdivgs(gel(x,i), n+i-1);
2571 52782 : return RgX_renormalize_lg(y, lx);;
2572 : }
2573 :
2574 : GEN
2575 52918 : RgXn_expint(GEN h, long e)
2576 : {
2577 52918 : pari_sp av = avma, av2;
2578 52918 : long v = varn(h), n;
2579 52918 : GEN f = pol_1(v), g;
2580 : ulong mask;
2581 :
2582 52917 : if (!signe(h)) return f;
2583 50425 : g = pol_1(v);
2584 50424 : n = 1; mask = quadratic_prec_mask(e);
2585 50424 : av2 = avma;
2586 56374 : for (;mask>1;)
2587 : {
2588 : GEN u, w;
2589 56374 : long n2 = n;
2590 56374 : n<<=1; if (mask & 1) n--;
2591 56374 : mask >>= 1;
2592 56374 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2593 56375 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2594 56373 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2595 56373 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2596 56376 : if (mask<=1) break;
2597 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2598 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2599 5950 : if (gc_needed(av2,2))
2600 : {
2601 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2602 0 : gerepileall(av2, 2, &f, &g);
2603 : }
2604 : }
2605 50426 : return gerepileupto(av, f);
2606 : }
2607 :
2608 : GEN
2609 0 : RgXn_exp(GEN h, long e)
2610 : {
2611 0 : long d = degpol(h);
2612 0 : if (d < 0) return pol_1(varn(h));
2613 0 : if (!d || !gequal0(gel(h,2)))
2614 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2615 0 : return RgXn_expint(RgX_deriv(h), e);
2616 : }
2617 :
2618 : GEN
2619 154 : RgXn_reverse(GEN f, long e)
2620 : {
2621 154 : pari_sp av = avma, av2;
2622 : ulong mask;
2623 : GEN fi, a, df, W, an;
2624 154 : long v = varn(f), n=1;
2625 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2626 0 : pari_err_INV("serreverse",f);
2627 154 : fi = ginv(gel(f,3));
2628 154 : a = deg1pol_shallow(fi,gen_0,v);
2629 154 : if (e <= 2) return gerepilecopy(av, a);
2630 133 : W = scalarpol(fi,v);
2631 133 : df = RgX_deriv(f);
2632 133 : mask = quadratic_prec_mask(e);
2633 133 : av2 = avma;
2634 616 : for (;mask>1;)
2635 : {
2636 : GEN u, fa, fr;
2637 483 : long n2 = n, rt;
2638 483 : n<<=1; if (mask & 1) n--;
2639 483 : mask >>= 1;
2640 483 : fr = RgXn_red_shallow(f, n);
2641 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2642 483 : an = RgXn_powers(a, rt, n);
2643 483 : if (n>1)
2644 : {
2645 483 : long n4 = (n2+1)>>1;
2646 483 : GEN dfr = RgXn_red_shallow(df, n2);
2647 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2648 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2649 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2650 : }
2651 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2652 483 : fa = RgX_shift(fa, -n2);
2653 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2654 483 : if (gc_needed(av2,2))
2655 : {
2656 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2657 0 : gerepileall(av2, 2, &a, &W);
2658 : }
2659 : }
2660 133 : return gerepileupto(av, a);
2661 : }
2662 :
2663 : GEN
2664 7 : RgXn_sqrt(GEN h, long e)
2665 : {
2666 7 : pari_sp av = avma, av2;
2667 7 : long v = varn(h), n = 1;
2668 7 : GEN f = scalarpol(gen_1, v), df = f;
2669 7 : ulong mask = quadratic_prec_mask(e);
2670 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2671 0 : pari_err_SQRTN("RgXn_sqrt",h);
2672 7 : av2 = avma;
2673 : while(1)
2674 7 : {
2675 14 : long n2 = n, m;
2676 : GEN g;
2677 14 : n<<=1; if (mask & 1) n--;
2678 14 : mask >>= 1;
2679 14 : m = n-n2;
2680 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2681 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2682 14 : if (mask==1) return gerepileupto(av, f);
2683 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2684 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2685 7 : if (gc_needed(av2,2))
2686 : {
2687 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2688 0 : gerepileall(av2, 2, &f, &df);
2689 : }
2690 : }
2691 : }
2692 :
2693 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2694 : GEN
2695 115807 : RgXQ_powu(GEN x, ulong n, GEN T)
2696 : {
2697 115807 : pari_sp av = avma;
2698 :
2699 115807 : if (!n) return pol_1(varn(x));
2700 67952 : if (n == 1) return RgX_copy(x);
2701 19740 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2702 19740 : return gerepilecopy(av, x);
2703 : }
2704 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2705 : GEN
2706 102412 : RgXQ_pow(GEN x, GEN n, GEN T)
2707 : {
2708 : pari_sp av;
2709 102412 : long s = signe(n);
2710 :
2711 102412 : if (!s) return pol_1(varn(x));
2712 102412 : if (is_pm1(n) == 1)
2713 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2714 14105 : av = avma;
2715 14105 : if (s < 0) x = RgXQ_inv(x, T);
2716 14105 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2717 14105 : return gerepilecopy(av, x);
2718 : }
2719 : static GEN
2720 200568 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2721 : static GEN
2722 111173 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2723 :
2724 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2725 : GEN
2726 11410 : ZXQ_powers(GEN x, long l, GEN T)
2727 : {
2728 11410 : int use_sqr = 2*degpol(x) >= degpol(T);
2729 11410 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2730 : }
2731 :
2732 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2733 : GEN
2734 168510 : ZXQ_powu(GEN x, ulong n, GEN T)
2735 : {
2736 168510 : pari_sp av = avma;
2737 :
2738 168510 : if (!n) return pol_1(varn(x));
2739 168510 : if (n == 1) return ZX_copy(x);
2740 113874 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2741 113873 : return gerepilecopy(av, x);
2742 : }
2743 :
2744 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2745 : GEN
2746 8547 : RgXQ_powers(GEN x, long l, GEN T)
2747 : {
2748 8547 : int use_sqr = 2*degpol(x) >= degpol(T);
2749 8547 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2750 : }
2751 :
2752 : GEN
2753 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2754 : {
2755 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2756 : }
2757 :
2758 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2759 : GEN
2760 8848 : QXQ_powers(GEN a, long n, GEN T)
2761 : {
2762 : GEN den, v;
2763 8848 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2764 8827 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2765 : /* den*a integral; v[i+1] = (den*a)^i in K */
2766 8827 : if (den)
2767 : { /* restore denominators */
2768 5435 : GEN d = den;
2769 : long i;
2770 5435 : gel(v,2) = a;
2771 27006 : for (i=3; i<=n+1; i++) {
2772 21571 : d = mulii(d,den);
2773 21571 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2774 : }
2775 : }
2776 8827 : return v;
2777 : }
2778 :
2779 : static GEN
2780 2898 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2781 : {
2782 2898 : long l, i, m = 0;
2783 : GEN dz, z;
2784 2898 : GEN V = cgetg_copy(v, &l);
2785 9884 : for (i = imin; i < l; i++)
2786 : {
2787 6986 : GEN c = gel(v, i);
2788 6986 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2789 : }
2790 2898 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2791 3129 : for (i = 1; i < imin; i++) V[i] = v[i];
2792 9884 : for (i = imin; i < l; i++)
2793 : {
2794 6986 : GEN c = gel(v,i);
2795 6986 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2796 6986 : gel(V,i) = c;
2797 : }
2798 2898 : return V;
2799 : }
2800 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2801 : GEN
2802 2667 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2803 2667 : { return do_QXQ_eval(v, 1, a, T); }
2804 :
2805 : GEN
2806 231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2807 231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2808 :
2809 : GEN
2810 2940 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2811 : {
2812 2940 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2813 : }
2814 :
2815 : GEN
2816 5144 : RgXQ_norm(GEN x, GEN T)
2817 : {
2818 : pari_sp av;
2819 5144 : long dx = degpol(x);
2820 : GEN L, y;
2821 5144 : if (degpol(T)==0) return gpowgs(x,0);
2822 5137 : av = avma; y = resultant(T, x);
2823 5137 : L = leading_coeff(T);
2824 5137 : if (gequal1(L) || !signe(x)) return y;
2825 0 : return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2826 : }
2827 :
2828 : GEN
2829 476 : RgXQ_trace(GEN x, GEN T)
2830 : {
2831 476 : pari_sp av = avma;
2832 : GEN dT, z;
2833 : long n;
2834 476 : if (degpol(T)==0) return gmulgs(x,0);
2835 469 : dT = RgX_deriv(T); n = degpol(dT);
2836 469 : z = RgXQ_mul(x, dT, T);
2837 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2838 420 : return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2839 : }
2840 :
2841 : GEN
2842 2787609 : RgX_blocks(GEN P, long n, long m)
2843 : {
2844 2787609 : GEN z = cgetg(m+1,t_VEC);
2845 2787609 : long i,j, k=2, l = lg(P);
2846 8556405 : for(i=1; i<=m; i++)
2847 : {
2848 5768796 : GEN zi = cgetg(n+2,t_POL);
2849 5768796 : zi[1] = P[1];
2850 5768796 : gel(z,i) = zi;
2851 16990226 : for(j=2; j<n+2; j++)
2852 11221430 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2853 5768796 : zi = RgX_renormalize_lg(zi, n+2);
2854 : }
2855 2787609 : return z;
2856 : }
2857 :
2858 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2859 : void
2860 49653117 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2861 : {
2862 49653117 : long n = degpol(p), v = varn(p), n0, n1, i;
2863 : GEN p0, p1;
2864 :
2865 49653791 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2866 :
2867 49560514 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2868 49560514 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2869 49562406 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2870 146477716 : for (i=0; i<n1; i++)
2871 : {
2872 96913969 : p0[2+i] = p[2+(i<<1)];
2873 96913969 : p1[2+i] = p[3+(i<<1)];
2874 : }
2875 49563747 : if (n1 != n0)
2876 17960898 : p0[2+i] = p[2+(i<<1)];
2877 49563747 : *pe = normalizepol(p0);
2878 49566672 : *po = normalizepol(p1);
2879 : }
2880 :
2881 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2882 : GEN
2883 43631 : RgX_splitting(GEN p, long k)
2884 : {
2885 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
2886 : GEN r;
2887 :
2888 43631 : m = n/k;
2889 43631 : r = cgetg(k+1,t_VEC);
2890 234059 : for(i=1; i<=k; i++)
2891 : {
2892 190428 : gel(r,i) = cgetg(m+3, t_POL);
2893 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2894 : }
2895 571991 : for (j=1, i=0, l=2; i<=n; i++)
2896 : {
2897 528360 : gmael(r,j,l) = gel(p,2+i);
2898 528360 : if (j==k) { j=1; l++; } else j++;
2899 : }
2900 234059 : for(i=1; i<=k; i++)
2901 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2902 43631 : return r;
2903 : }
2904 :
2905 : /*******************************************************************/
2906 : /* */
2907 : /* Kronecker form */
2908 : /* */
2909 : /*******************************************************************/
2910 :
2911 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2912 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2913 : * normalized coefficients */
2914 : GEN
2915 186343 : Kronecker_to_mod(GEN z, GEN T)
2916 : {
2917 186343 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2918 186343 : GEN x, t = cgetg(N,t_POL);
2919 186343 : t[1] = T[1];
2920 186343 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2921 186343 : x[1] = z[1];
2922 186343 : T = RgX_copy(T);
2923 1405143 : for (i=2; i<lx+2; i++, z+= N-2)
2924 : {
2925 5792558 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2926 1218800 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2927 : }
2928 186343 : N = (l-2) % (N-2) + 2;
2929 473273 : for (j=2; j<N; j++) t[j] = z[j];
2930 186343 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2931 186343 : return normalizepol_lg(x, i+1);
2932 : }
2933 :
2934 : /*******************************************************************/
2935 : /* */
2936 : /* Domain detection */
2937 : /* */
2938 : /*******************************************************************/
2939 :
2940 : static GEN
2941 152790 : zero_FpX_mod(GEN p, long v)
2942 : {
2943 152790 : GEN r = cgetg(3,t_POL);
2944 152790 : r[1] = evalvarn(v);
2945 152790 : gel(r,2) = mkintmod(gen_0, icopy(p));
2946 152790 : return r;
2947 : }
2948 :
2949 : static GEN
2950 818173 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2951 : {
2952 818173 : pari_sp av = avma;
2953 : GEN r;
2954 818173 : if (lgefint(p) == 3)
2955 : {
2956 766817 : ulong pp = uel(p, 2);
2957 766817 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2958 : RgX_to_Flx(y, pp), pp));
2959 : }
2960 : else
2961 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2962 818173 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2963 703120 : return gerepileupto(av, FpX_to_mod(r, p));
2964 : }
2965 :
2966 : static GEN
2967 532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
2968 : {
2969 532 : GEN r = cgetg(3,t_POL);
2970 532 : r[1] = evalvarn(v);
2971 532 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
2972 532 : return r;
2973 : }
2974 :
2975 : static GEN
2976 2044 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
2977 : {
2978 2044 : pari_sp av = avma;
2979 : long dT;
2980 : GEN kx, ky, r;
2981 2044 : GEN T = RgX_to_FpX(pol, p);
2982 2044 : if (signe(T)==0) pari_err_OP("*", x, y);
2983 2037 : dT = degpol(T);
2984 2037 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
2985 2037 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
2986 2037 : r = FpX_mul(kx, ky, p);
2987 2037 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
2988 1512 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
2989 : }
2990 :
2991 : static GEN
2992 440280 : RgX_liftred(GEN x, GEN T)
2993 440280 : { return RgXQX_red(liftpol_shallow(x), T); }
2994 :
2995 : static GEN
2996 171210 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
2997 : {
2998 171210 : pari_sp av = avma;
2999 171210 : long dT = degpol(T);
3000 171210 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3001 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3002 171210 : return gerepileupto(av, Kronecker_to_mod(r, T));
3003 : }
3004 :
3005 : static GEN
3006 1421 : RgX_sqr_FpX(GEN x, GEN p)
3007 : {
3008 1421 : pari_sp av = avma;
3009 : GEN r;
3010 1421 : if (lgefint(p) == 3)
3011 : {
3012 1217 : ulong pp = uel(p, 2);
3013 1217 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3014 : }
3015 : else
3016 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3017 1421 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3018 1281 : return gerepileupto(av, FpX_to_mod(r, p));
3019 : }
3020 :
3021 : static GEN
3022 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3023 : {
3024 196 : pari_sp av = avma;
3025 : long dT;
3026 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3027 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3028 189 : dT = degpol(T);
3029 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3030 189 : r = FpX_sqr(kx, p);
3031 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3032 189 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3033 : }
3034 :
3035 : static GEN
3036 13432 : RgX_sqr_QXQX(GEN x, GEN T)
3037 : {
3038 13432 : pari_sp av = avma;
3039 13432 : long dT = degpol(T);
3040 13432 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3041 13432 : return gerepileupto(av, Kronecker_to_mod(r, T));
3042 : }
3043 :
3044 : static GEN
3045 68670 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3046 : {
3047 68670 : pari_sp av = avma;
3048 : GEN r;
3049 68670 : if (lgefint(p) == 3)
3050 : {
3051 51880 : ulong pp = uel(p, 2);
3052 51880 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3053 : RgX_to_Flx(y, pp), pp));
3054 : }
3055 : else
3056 16790 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3057 68670 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3058 31073 : return gerepileupto(av, FpX_to_mod(r, p));
3059 : }
3060 :
3061 : static GEN
3062 42214 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3063 : {
3064 42214 : pari_sp av = avma;
3065 : GEN r;
3066 42214 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3067 42214 : return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
3068 : }
3069 : static GEN
3070 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3071 : {
3072 70 : pari_sp av = avma;
3073 : GEN r;
3074 70 : GEN T = RgX_to_FpX(pol, p);
3075 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3076 63 : if (lgefint(p) == 3)
3077 : {
3078 55 : ulong pp = uel(p, 2);
3079 55 : GEN Tp = ZX_to_Flx(T, pp);
3080 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3081 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3082 : }
3083 : else
3084 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3085 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3086 56 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3087 : }
3088 :
3089 : #define code(t1,t2) ((t1 << 6) | t2)
3090 : static GEN
3091 77694048 : RgX_mul_fast(GEN x, GEN y)
3092 : {
3093 : GEN p, pol;
3094 : long pa;
3095 77694048 : long t = RgX_type2(x,y, &p,&pol,&pa);
3096 77694354 : switch(t)
3097 : {
3098 41495210 : case t_INT: return ZX_mul(x,y);
3099 3327190 : case t_FRAC: return QX_mul(x,y);
3100 104735 : case t_FFELT: return FFX_mul(x, y, pol);
3101 818173 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3102 171210 : case code(t_POLMOD, t_INT):
3103 : case code(t_POLMOD, t_FRAC):
3104 171210 : return RgX_mul_QXQX(x, y, pol);
3105 2007 : case code(t_POLMOD, t_INTMOD):
3106 2007 : return RgX_mul_FpXQX(x, y, pol, p);
3107 31775829 : default: return NULL;
3108 : }
3109 : }
3110 : static GEN
3111 1359308 : RgX_sqr_fast(GEN x)
3112 : {
3113 : GEN p, pol;
3114 : long pa;
3115 1359308 : long t = RgX_type(x,&p,&pol,&pa);
3116 1359311 : switch(t)
3117 : {
3118 1253467 : case t_INT: return ZX_sqr(x);
3119 79395 : case t_FRAC: return QX_sqr(x);
3120 2534 : case t_FFELT: return FFX_sqr(x, pol);
3121 1421 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3122 13432 : case code(t_POLMOD, t_INT):
3123 : case code(t_POLMOD, t_FRAC):
3124 13432 : return RgX_sqr_QXQX(x, pol);
3125 193 : case code(t_POLMOD, t_INTMOD):
3126 193 : return RgX_sqr_FpXQX(x, pol, p);
3127 8869 : default: return NULL;
3128 : }
3129 : }
3130 :
3131 : static GEN
3132 14520246 : RgX_rem_fast(GEN x, GEN y)
3133 : {
3134 : GEN p, pol;
3135 : long pa;
3136 14520246 : long t = RgX_type2(x,y, &p,&pol,&pa);
3137 14520387 : switch(t)
3138 : {
3139 3381638 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3140 1793998 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3141 84 : case t_FFELT: return FFX_rem(x, y, pol);
3142 68670 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3143 42214 : case code(t_POLMOD, t_INT):
3144 : case code(t_POLMOD, t_FRAC):
3145 42214 : return RgX_rem_QXQX(x, y, pol);
3146 95 : case code(t_POLMOD, t_INTMOD):
3147 95 : return RgX_rem_FpXQX(x, y, pol, p);
3148 9233688 : default: return NULL;
3149 : }
3150 : }
3151 :
3152 : #undef code
3153 :
3154 : GEN
3155 63962019 : RgX_mul(GEN x, GEN y)
3156 : {
3157 63962019 : GEN z = RgX_mul_fast(x,y);
3158 63962211 : if (!z) z = RgX_mul_i(x,y);
3159 63961743 : return z;
3160 : }
3161 :
3162 : GEN
3163 1347156 : RgX_sqr(GEN x)
3164 : {
3165 1347156 : GEN z = RgX_sqr_fast(x);
3166 1347153 : if (!z) z = RgX_sqr_i(x);
3167 1347153 : return z;
3168 : }
3169 :
3170 : GEN
3171 14520253 : RgX_rem(GEN x, GEN y)
3172 : {
3173 14520253 : GEN z = RgX_rem_fast(x, y);
3174 14520363 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3175 14520280 : return z;
3176 : }
3177 :
3178 : static GEN
3179 0 : _RgX_mul(void* E, GEN x, GEN y)
3180 0 : { (void) E; return RgX_mul(x, y); }
3181 :
3182 : GEN
3183 0 : RgXV_prod(GEN V)
3184 0 : { return gen_product(V, NULL, &_RgX_mul); }
3185 :
3186 : GEN
3187 11061157 : RgXn_mul(GEN f, GEN g, long n)
3188 : {
3189 11061157 : pari_sp av = avma;
3190 11061157 : GEN h = RgX_mul_fast(f,g);
3191 11061155 : if (!h) return RgXn_mul2(f,g,n);
3192 1783409 : if (degpol(h) < n) return h;
3193 382941 : return gerepilecopy(av, RgXn_red_shallow(h, n));
3194 : }
3195 :
3196 : GEN
3197 12152 : RgXn_sqr(GEN f, long n)
3198 : {
3199 12152 : pari_sp av = avma;
3200 12152 : GEN g = RgX_sqr_fast(f);
3201 12152 : if (!g) return RgXn_sqr2(f,n);
3202 11837 : if (degpol(g) < n) return g;
3203 10640 : return gerepilecopy(av, RgXn_red_shallow(g, n));
3204 : }
3205 :
3206 : /* (f*g) \/ x^n */
3207 : GEN
3208 2670872 : RgX_mulhigh_i(GEN f, GEN g, long n)
3209 : {
3210 2670872 : GEN h = RgX_mul_fast(f,g);
3211 2670873 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3212 : }
3213 :
3214 : /* (f*g) \/ x^n */
3215 : GEN
3216 0 : RgX_sqrhigh_i(GEN f, long n)
3217 : {
3218 0 : GEN h = RgX_sqr_fast(f);
3219 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3220 : }
|
