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 19590324 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19590324 : long pold=1, rold=n*(d-1);
35 91766944 : for(p=2; p<=d; p++)
36 : {
37 72176620 : r = m*(p-1) + n*((d-1)/p);
38 72176620 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19590324 : return pold;
41 : }
42 :
43 : static GEN
44 14622312 : 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 14622312 : pari_sp av = avma;
48 : long i;
49 14622312 : GEN z = cmul(E,P,a,ff->one(E));
50 14594894 : if (!z) z = gen_0;
51 76538763 : for (i=1; i<=n; i++)
52 : {
53 61942821 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 61970765 : if (t) {
55 47038511 : z = ff->add(E, z, t);
56 47002717 : if (gc_needed(av,2)) z = gc_upto(av, z);
57 : }
58 : }
59 14595942 : 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 12526699 : 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 12526699 : pari_sp av = avma;
72 12526699 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12526699 : if (d < 0) return ff->zero(E);
76 11916547 : if (d < l) return gc_upto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1158676 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1158676 : 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 1158676 : d -= l;
84 1158676 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2704921 : while (d >= l-1)
86 : {
87 1545006 : d -= l-1;
88 1545006 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1545032 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1544967 : if (gc_needed(av,2))
91 91 : z = gc_upto(av, z);
92 : }
93 1159915 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1159929 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1159917 : return gc_upto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 875617 : 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 875617 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 875617 : if (d < 0) return ff->zero(E);
106 874322 : rtd = (long) sqrt((double)d);
107 874322 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 874338 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 874325 : return gc_upto(av, z);
110 : }
111 :
112 : static GEN
113 2033649 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19257653 : _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 1844896 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 593313 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2074862 : _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 887905 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 887905 : pari_sp av = avma, btop;
163 887905 : long i, d = degpol(P), o;
164 : GEN s;
165 887902 : if (signe(P)==0) return pol_0(varn(P));
166 887900 : s = gel(P, d+2);
167 887900 : if (d == 0) return gcopy(s);
168 884624 : o = RgX_deflate_order(P);
169 884639 : if (o > 1) A = gpowgs(A, o);
170 884652 : btop = avma;
171 2971452 : for (i = d-o; i >= 0; i-=o)
172 : {
173 2086807 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 2086800 : if (gc_needed(btop,1))
175 : {
176 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
177 13 : s = gc_upto(btop, s);
178 : }
179 : }
180 884645 : return gc_upto(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 = gc_upto(av, s);
200 : }
201 : }
202 1232 : return gc_upto(av, s);
203 : }
204 :
205 : const struct bb_algebra *
206 282220 : get_Rg_algebra(void)
207 : {
208 282220 : return &Rg_algebra;
209 : }
210 :
211 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
212 :
213 : static GEN
214 11333 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
215 : {
216 : (void) E;
217 11333 : return RgX_divrem(x, y, r);
218 : }
219 :
220 : GEN
221 3157 : RgX_digits(GEN x, GEN T)
222 : {
223 3157 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
224 3157 : if (signe(x)==0) return(cgetg(1, t_VEC));
225 3157 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
226 : }
227 :
228 : /*******************************************************************/
229 : /* */
230 : /* RgX */
231 : /* */
232 : /*******************************************************************/
233 :
234 : long
235 24892503 : RgX_equal(GEN x, GEN y)
236 : {
237 24892503 : long i = lg(x);
238 :
239 24892503 : if (i != lg(y)) return 0;
240 107136958 : for (i--; i > 1; i--)
241 82572629 : if (!gequal(gel(x,i),gel(y,i))) return 0;
242 24564329 : return 1;
243 : }
244 :
245 : /* Returns 1 in the base ring over which x is defined */
246 : /* HACK: this also works for t_SER */
247 : GEN
248 156950330 : Rg_get_1(GEN x)
249 : {
250 : GEN p, T;
251 156950330 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
252 156950331 : if (RgX_type_is_composite(tx))
253 11657616 : RgX_type_decode(tx, &i /*junk*/, &tx);
254 156950331 : switch(tx)
255 : {
256 427406 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
257 945 : case t_PADIC: return cvtop(gen_1, p, lx);
258 5950 : case t_FFELT: return FF_1(T);
259 156516030 : default: return gen_1;
260 : }
261 : }
262 : /* Returns 0 in the base ring over which x is defined */
263 : /* HACK: this also works for t_SER */
264 : GEN
265 6942779 : Rg_get_0(GEN x)
266 : {
267 : GEN p, T;
268 6942779 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
269 6942779 : if (RgX_type_is_composite(tx))
270 61341 : RgX_type_decode(tx, &i /*junk*/, &tx);
271 6942779 : switch(tx)
272 : {
273 497 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
274 42 : case t_PADIC: return zeropadic(p, lx);
275 210 : case t_FFELT: return FF_zero(T);
276 6942030 : default: return gen_0;
277 : }
278 : }
279 :
280 : GEN
281 7595 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
282 : {
283 7595 : long i, n = degpol(P);
284 : GEN z, dz, dP;
285 7595 : if (n < 0) return gen_0;
286 7595 : P = Q_remove_denom(P, &dP);
287 7595 : z = gel(P,2); if (n == 0) return icopy(z);
288 4361 : if (dV) z = mulii(dV, z); /* V[1] = dV */
289 4361 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
290 8169 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
291 4361 : dz = mul_denom(dP, dV);
292 4361 : return dz? RgX_Rg_div(z, dz): z;
293 : }
294 :
295 : /* Return P(h * x), not memory clean */
296 : GEN
297 51885 : RgX_unscale(GEN P, GEN h)
298 : {
299 51885 : long i, l = lg(P);
300 51885 : GEN hi = gen_1, Q = cgetg(l, t_POL);
301 51885 : Q[1] = P[1];
302 51885 : if (l == 2) return Q;
303 40643 : gel(Q,2) = gcopy(gel(P,2));
304 95893 : for (i=3; i<l; i++)
305 : {
306 55251 : hi = gmul(hi,h);
307 55250 : gel(Q,i) = gmul(gel(P,i), hi);
308 : }
309 40642 : return Q;
310 : }
311 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
312 : GEN
313 1485103 : ZX_z_unscale(GEN P, long h)
314 : {
315 1485103 : long i, l = lg(P);
316 1485103 : GEN Q = cgetg(l, t_POL);
317 1485104 : Q[1] = P[1];
318 1485104 : if (l == 2) return Q;
319 1396558 : gel(Q,2) = gel(P,2);
320 1396558 : if (l == 3) return Q;
321 1370252 : if (h == -1)
322 241760 : for (i = 3; i < l; i++)
323 : {
324 199017 : gel(Q,i) = negi(gel(P,i));
325 199016 : if (++i == l) break;
326 148019 : gel(Q,i) = gel(P,i);
327 : }
328 : else
329 : {
330 : GEN hi;
331 1276511 : gel(Q,3) = mulis(gel(P,3), h);
332 1276510 : hi = sqrs(h);
333 5447143 : for (i = 4; i < l; i++)
334 : {
335 4170632 : gel(Q,i) = mulii(gel(P,i), hi);
336 4170633 : if (i != l-1) hi = mulis(hi,h);
337 : }
338 : }
339 1370251 : return Q;
340 : }
341 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
342 : GEN
343 1013197 : ZX_unscale(GEN P, GEN h)
344 : {
345 : long i, l;
346 : GEN Q, hi;
347 1013197 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
348 881 : l = lg(P); Q = cgetg(l, t_POL);
349 881 : Q[1] = P[1];
350 881 : if (l == 2) return Q;
351 881 : gel(Q,2) = gel(P,2);
352 881 : if (l == 3) return Q;
353 881 : hi = h;
354 881 : gel(Q,3) = mulii(gel(P,3), hi);
355 2741 : for (i = 4; i < l; i++)
356 : {
357 1860 : hi = mulii(hi,h);
358 1860 : gel(Q,i) = mulii(gel(P,i), hi);
359 : }
360 881 : return Q;
361 : }
362 : /* P a ZX. Return P(x << n), not memory clean */
363 : GEN
364 1113428 : ZX_unscale2n(GEN P, long n)
365 : {
366 1113428 : long i, ni = n, l = lg(P);
367 1113428 : GEN Q = cgetg(l, t_POL);
368 1113427 : Q[1] = P[1];
369 1113427 : if (l == 2) return Q;
370 1113427 : gel(Q,2) = gel(P,2);
371 1113427 : if (l == 3) return Q;
372 1113427 : gel(Q,3) = shifti(gel(P,3), ni);
373 3603528 : for (i=4; i<l; i++)
374 : {
375 2490154 : ni += n;
376 2490154 : gel(Q,i) = shifti(gel(P,i), ni);
377 : }
378 1113374 : return Q;
379 : }
380 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
381 : GEN
382 12997 : ZX_unscale_div(GEN P, GEN h)
383 : {
384 12997 : long i, l = lg(P);
385 12997 : GEN hi, Q = cgetg(l, t_POL);
386 12997 : Q[1] = P[1];
387 12997 : if (l == 2) return Q;
388 12997 : gel(Q,2) = diviiexact(gel(P,2), h);
389 12997 : if (l == 3) return Q;
390 12997 : gel(Q,3) = gel(P,3);
391 12997 : if (l == 4) return Q;
392 12997 : hi = h;
393 12997 : gel(Q,4) = mulii(gel(P,4), hi);
394 64121 : for (i=5; i<l; i++)
395 : {
396 51124 : hi = mulii(hi,h);
397 51124 : gel(Q,i) = mulii(gel(P,i), hi);
398 : }
399 12997 : return Q;
400 : }
401 : /* P(h*X) / h^k, assuming the result is a ZX */
402 : GEN
403 1393 : ZX_unscale_divpow(GEN P, GEN h, long k)
404 : {
405 1393 : long i, j, l = lg(P);
406 1393 : GEN H, Q = cgetg(l, t_POL);
407 1393 : Q[1] = P[1]; if (l == 2) return Q;
408 1393 : H = gpowers(h, maxss(k, l - 3 - k));
409 5572 : for (i = 2, j = k+1; j > 1 && i < l; i++)
410 4179 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
411 1393 : if (i == l) return Q;
412 1393 : gel(Q, i) = gel(P, i); i++;
413 5082 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
414 1393 : return Q;
415 : }
416 :
417 : GEN
418 6559 : RgXV_unscale(GEN x, GEN h)
419 : {
420 6559 : if (isint1(h)) return gcopy(x);
421 18735 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
422 : }
423 :
424 : /* Return h^degpol(P) P(x / h), not memory clean */
425 : GEN
426 4343462 : RgX_rescale(GEN P, GEN h)
427 : {
428 4343462 : long i, l = lg(P);
429 4343462 : GEN Q = cgetg(l,t_POL), hi = h;
430 4343454 : gel(Q,l-1) = gel(P,l-1);
431 11430490 : for (i=l-2; i>=2; i--)
432 : {
433 11427886 : gel(Q,i) = gmul(gel(P,i), hi);
434 11427779 : if (i == 2) break;
435 7086750 : hi = gmul(hi,h);
436 : }
437 4343633 : Q[1] = P[1]; return Q;
438 : }
439 :
440 : GEN
441 2401 : RgXV_rescale(GEN x, GEN h)
442 : {
443 2401 : if (isint1(h)) return RgX_copy(x);
444 16086 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
445 : }
446 :
447 : /* A(X^d) --> A(X) */
448 : GEN
449 1181415 : RgX_deflate(GEN x0, long d)
450 : {
451 : GEN z, y, x;
452 1181415 : long i,id, dy, dx = degpol(x0);
453 1181414 : if (d == 1 || dx <= 0) return x0;
454 461795 : dy = dx/d;
455 461795 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
456 461794 : z = y + 2;
457 461794 : x = x0+ 2;
458 1771917 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
459 461794 : return y;
460 : }
461 :
462 : GEN
463 1260 : RgX_homogenize_deg(GEN P, long d, long v)
464 : {
465 : long i, l;
466 1260 : GEN Q = cgetg_copy(P, &l);
467 1260 : Q[1] = P[1];
468 4018 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
469 1260 : return Q;
470 : }
471 :
472 : GEN
473 18732 : RgX_homogenize(GEN P, long v)
474 : {
475 : long i, l, d;
476 18732 : GEN Q = cgetg_copy(P, &l);
477 18732 : Q[1] = P[1]; d = l-3;
478 165998 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
479 18732 : return Q;
480 : }
481 :
482 : /* F a t_RFRAC */
483 : long
484 140 : rfrac_deflate_order(GEN F)
485 : {
486 140 : GEN N = gel(F,1), D = gel(F,2);
487 140 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
488 140 : if (m == 1) return 1;
489 49 : if (typ(N) == t_POL && varn(N) == varn(D))
490 28 : m = cgcd(m, RgX_deflate_order(N));
491 49 : return m;
492 : }
493 : /* F a t_RFRAC */
494 : GEN
495 140 : rfrac_deflate_max(GEN F, long *m)
496 : {
497 140 : *m = rfrac_deflate_order(F);
498 140 : return rfrac_deflate(F, *m);
499 : }
500 : /* F a t_RFRAC */
501 : GEN
502 140 : rfrac_deflate(GEN F, long m)
503 : {
504 140 : GEN N = gel(F,1), D = gel(F,2);
505 140 : if (m == 1) return F;
506 49 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
507 49 : D = RgX_deflate(D, m); return mkrfrac(N, D);
508 : }
509 :
510 : /* return x0(X^d) */
511 : GEN
512 883187 : RgX_inflate(GEN x0, long d)
513 : {
514 883187 : long i, id, dy, dx = degpol(x0);
515 883187 : GEN x = x0 + 2, z, y;
516 883187 : if (dx <= 0) return leafcopy(x0);
517 805890 : dy = dx*d;
518 805890 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
519 805890 : z = y + 2;
520 28095155 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
521 11496982 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
522 805890 : return y;
523 : }
524 :
525 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
526 : static GEN
527 5702420 : RgX_Rg_translate_basecase(GEN P, GEN c)
528 : {
529 5702420 : pari_sp av = avma;
530 : GEN Q;
531 : long i, k, n;
532 :
533 5702420 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
534 5700488 : Q = leafcopy(P); n = degpol(P);
535 5700488 : if (isint1(c))
536 : {
537 7595 : for (i=1; i<=n; i++)
538 : {
539 20265 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gel(Q,2+k+1));
540 5390 : if (gc_needed(av,2))
541 : {
542 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(1), i = %ld/%ld", i,n);
543 0 : Q = gc_GEN(av, Q);
544 : }
545 : }
546 : }
547 5698284 : else if (isintm1(c))
548 : {
549 15974 : for (i=1; i<=n; i++)
550 : {
551 49630 : for (k=n-i; k<n; k++) gel(Q,2+k) = gsub(gel(Q,2+k), gel(Q,2+k+1));
552 12299 : if (gc_needed(av,2))
553 : {
554 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(-1), i = %ld/%ld", i,n);
555 0 : Q = gc_GEN(av, Q);
556 : }
557 : }
558 : }
559 : else
560 : {
561 19427156 : for (i=1; i<=n; i++)
562 : {
563 45455159 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1)));
564 13732548 : if (gc_needed(av,2))
565 : {
566 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate, i = %ld/%ld", i,n);
567 0 : Q = gc_GEN(av, Q);
568 : }
569 : }
570 : }
571 5700203 : return gc_GEN(av, Q);
572 : }
573 :
574 : static GEN
575 102627 : zero_FpX_mod(GEN p, long v)
576 : {
577 102627 : GEN r = cgetg(3,t_POL);
578 102627 : r[1] = evalvarn(v);
579 102627 : gel(r,2) = mkintmod(gen_0, icopy(p));
580 102627 : return r;
581 : }
582 :
583 : static GEN
584 616 : zero_FpXQX_mod(GEN pol, GEN p, long v)
585 : {
586 616 : GEN r = cgetg(3,t_POL);
587 616 : r[1] = evalvarn(v);
588 616 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
589 616 : return r;
590 : }
591 :
592 : static GEN
593 0 : RgX_Rg_translate_FpX(GEN P, GEN c, GEN p)
594 : {
595 0 : pari_sp av = avma;
596 : GEN r;
597 : #if 0
598 : /* 'divide by 0' error if p is not prime and c not invertible */
599 : if (lgefint(p) == 3)
600 : {
601 : ulong pp = uel(p, 2);
602 : r = Flx_to_ZX_inplace(Flx_Fl_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
603 : }
604 : else
605 : #endif
606 0 : r = FpX_Fp_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
607 0 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
608 0 : return gc_upto(av, FpX_to_mod(r, p));
609 : }
610 :
611 : static GEN
612 7 : RgX_Rg_translate_FpXQX(GEN x, GEN c, GEN pol, GEN p)
613 : {
614 7 : pari_sp av = avma;
615 7 : GEN r, T = RgX_to_FpX(pol, p);
616 7 : if (signe(T)==0) pari_err_OP("subst", x, c);
617 7 : r = FpXQX_FpXQ_translate(RgX_to_FpXQX(x, T, p), Rg_to_FpXQ(c, T, p), T, p);
618 7 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
619 7 : return gc_upto(av, FpXQX_to_mod(r, T, p));
620 : }
621 :
622 : static GEN
623 6064157 : RgX_Rg_translate_fast(GEN P, GEN c)
624 : {
625 : GEN p, pol;
626 : long pa;
627 6064157 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
628 6064154 : switch(t)
629 : {
630 361918 : case t_INT: return ZX_Z_translate(P, c);
631 0 : case t_INTMOD: return RgX_Rg_translate_FpX(P, c, p);
632 6 : case RgX_type_code(t_POLMOD, t_INTMOD):
633 6 : return RgX_Rg_translate_FpXQX(P, c, pol, p);
634 5702232 : default: return NULL;
635 : }
636 : }
637 :
638 : static GEN
639 5702608 : RgX_Rg_translate_i(GEN P, GEN c)
640 : {
641 5702608 : pari_sp av = avma;
642 : long n;
643 5702608 : n = degpol(P);
644 5702609 : if (n < 40)
645 5702420 : return RgX_Rg_translate_basecase(P, c);
646 : else
647 : {
648 189 : long d = n >> 1;
649 189 : GEN Q = RgX_Rg_translate_i(RgX_shift_shallow(P, -d), c);
650 189 : GEN R = RgX_Rg_translate_i(RgXn_red_shallow(P, d), c);
651 189 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
652 189 : return gc_upto(av, RgX_add(RgX_mul(Q, S), R));
653 : }
654 : }
655 :
656 : GEN
657 6064158 : RgX_Rg_translate(GEN P, GEN c)
658 : {
659 6064158 : GEN R = RgX_Rg_translate_fast(P, c);
660 6064147 : return R ? R: RgX_Rg_translate_i(P,c);
661 : }
662 : /* P(ax + b) */
663 : GEN
664 30212 : RgX_affine(GEN P, GEN a, GEN b)
665 : {
666 30212 : if (!gequal0(b)) P = RgX_Rg_translate(P, b);
667 30212 : return RgX_unscale(P, a);
668 : }
669 :
670 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
671 : GEN
672 33584 : RgXQX_RgXQ_translate(GEN P, GEN c, GEN T)
673 : {
674 33584 : pari_sp av = avma;
675 : GEN Q;
676 : long i, k, n;
677 :
678 33584 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
679 33241 : Q = leafcopy(P); n = degpol(P);
680 105778 : for (i=1; i<=n; i++)
681 : {
682 303686 : for (k=n-i; k<n; k++)
683 : {
684 231149 : pari_sp av2 = avma;
685 231149 : gel(Q,2+k) = gc_upto(av2,
686 231149 : RgX_rem(gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1))), T));
687 : }
688 72537 : if (gc_needed(av,2))
689 : {
690 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_RgXQ_translate, i = %ld/%ld", i,n);
691 0 : Q = gc_GEN(av, Q);
692 : }
693 : }
694 33241 : return gc_GEN(av, Q);
695 : }
696 :
697 : /********************************************************************/
698 : /** **/
699 : /** CONVERSIONS **/
700 : /** (not memory clean) **/
701 : /** **/
702 : /********************************************************************/
703 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
704 : * but everything else is */
705 : static GEN
706 169946 : QXQ_to_mod(GEN x, GEN T)
707 : {
708 : long d;
709 169946 : switch(typ(x))
710 : {
711 64796 : case t_INT: return icopy(x);
712 2079 : case t_FRAC: return gcopy(x);
713 103071 : case t_POL:
714 103071 : d = degpol(x);
715 103071 : if (d < 0) return gen_0;
716 100810 : if (d == 0) return gcopy(gel(x,2));
717 98671 : return mkpolmod(RgX_copy(x), T);
718 0 : default: pari_err_TYPE("QXQ_to_mod",x);
719 : return NULL;/* LCOV_EXCL_LINE */
720 : }
721 : }
722 : /* pure shallow version */
723 : GEN
724 840905 : QXQ_to_mod_shallow(GEN x, GEN T)
725 : {
726 : long d;
727 840905 : switch(typ(x))
728 : {
729 544756 : case t_INT:
730 544756 : case t_FRAC: return x;
731 296149 : case t_POL:
732 296149 : d = degpol(x);
733 296149 : if (d < 0) return gen_0;
734 249933 : if (d == 0) return gel(x,2);
735 233361 : return mkpolmod(x, T);
736 0 : default: pari_err_TYPE("QXQ_to_mod",x);
737 : return NULL;/* LCOV_EXCL_LINE */
738 : }
739 : }
740 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
741 : * Not memory clean because T not copied, but everything else is */
742 : static GEN
743 37471 : QXQX_to_mod(GEN z, GEN T)
744 : {
745 37471 : long i,l = lg(z);
746 37471 : GEN x = cgetg(l,t_POL);
747 187922 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
748 37471 : x[1] = z[1]; return normalizepol_lg(x,l);
749 : }
750 : /* pure shallow version */
751 : GEN
752 204271 : QXQX_to_mod_shallow(GEN z, GEN T)
753 : {
754 204271 : long i,l = lg(z);
755 204271 : GEN x = cgetg(l,t_POL);
756 961722 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
757 204271 : x[1] = z[1]; return normalizepol_lg(x,l);
758 : }
759 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
760 : GEN
761 12733 : QXQXV_to_mod(GEN V, GEN T)
762 : {
763 12733 : long i, l = lg(V);
764 12733 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
765 50204 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
766 12733 : return z;
767 : }
768 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
769 : GEN
770 22339 : QXQV_to_mod(GEN V, GEN T)
771 : {
772 22339 : long i, l = lg(V);
773 22339 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
774 41834 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
775 22339 : return z;
776 : }
777 :
778 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
779 : GEN
780 14854 : QXQC_to_mod_shallow(GEN V, GEN T)
781 : {
782 14854 : long i, l = lg(V);
783 14854 : GEN z = cgetg(l, t_COL);
784 98308 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
785 14854 : return z;
786 : }
787 :
788 : GEN
789 6720 : QXQM_to_mod_shallow(GEN V, GEN T)
790 : {
791 6720 : long i, l = lg(V);
792 6720 : GEN z = cgetg(l, t_MAT);
793 21574 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
794 6720 : return z;
795 : }
796 :
797 : GEN
798 7886625 : RgX_renormalize_lg(GEN x, long lx)
799 : {
800 : long i;
801 10938473 : for (i = lx-1; i>1; i--)
802 10485212 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
803 7886625 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
804 7886625 : setlg(x, i+1); setsigne(x, i != 1); return x;
805 : }
806 :
807 : GEN
808 1501742 : RgV_to_RgX(GEN x, long v)
809 : {
810 1501742 : long i, k = lg(x);
811 : GEN p;
812 :
813 4191163 : while (--k && gequal0(gel(x,k)));
814 1501741 : if (!k) return pol_0(v);
815 1494127 : i = k+2; p = cgetg(i,t_POL);
816 1494123 : p[1] = evalsigne(1) | evalvarn(v);
817 10289689 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
818 1494123 : return p;
819 : }
820 : GEN
821 188283 : RgV_to_RgX_reverse(GEN x, long v)
822 : {
823 188283 : long j, k, l = lg(x);
824 : GEN p;
825 :
826 189851 : for (k = 1; k < l; k++)
827 189851 : if (!gequal0(gel(x,k))) break;
828 188283 : if (k == l) return pol_0(v);
829 188283 : k -= 1;
830 188283 : l -= k;
831 188283 : x += k;
832 188283 : p = cgetg(l+1,t_POL);
833 188283 : p[1] = evalsigne(1) | evalvarn(v);
834 985038 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
835 188283 : return p;
836 : }
837 :
838 : /* return the (N-dimensional) vector of coeffs of p */
839 : GEN
840 13426169 : RgX_to_RgC(GEN x, long N)
841 : {
842 : long i, l;
843 : GEN z;
844 13426169 : l = lg(x)-1; x++;
845 13426169 : if (l > N+1) l = N+1; /* truncate higher degree terms */
846 13426169 : z = cgetg(N+1,t_COL);
847 80980378 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
848 24514606 : for ( ; i<=N; i++) gel(z,i) = gen_0;
849 13426258 : return z;
850 : }
851 : GEN
852 1360770 : Rg_to_RgC(GEN x, long N)
853 : {
854 1360770 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
855 : }
856 :
857 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
858 : GEN
859 303235 : RgM_to_RgXV(GEN x, long v)
860 1290955 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
861 : GEN
862 7202 : RgM_to_RgXV_reverse(GEN x, long v)
863 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
864 :
865 : /* matrix whose entries are given by the coeffs of the polynomials in
866 : * vector v (considered as degree n-1 polynomials) */
867 : GEN
868 336703 : RgV_to_RgM(GEN x, long n)
869 1692276 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
870 :
871 : GEN
872 79032 : RgXV_to_RgM(GEN x, long n)
873 401865 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
874 :
875 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
876 : GEN
877 23714 : RgM_to_RgXX(GEN x, long v,long w)
878 : {
879 23714 : long j, lx = lg(x);
880 23714 : GEN y = cgetg(lx+1, t_POL);
881 23714 : y[1] = evalsigne(1) | evalvarn(v);
882 23714 : y++;
883 131910 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
884 23714 : return normalizepol_lg(--y, lx+1);
885 : }
886 :
887 : /* matrix whose entries are given by the coeffs of the polynomial v in
888 : * two variables (considered as degree n-1 polynomials) */
889 : GEN
890 322 : RgXX_to_RgM(GEN v, long n)
891 : {
892 322 : long j, N = lg(v)-1;
893 322 : GEN y = cgetg(N, t_MAT);
894 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
895 322 : return y;
896 : }
897 :
898 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
899 : GEN
900 32767 : RgXY_swapspec(GEN x, long n, long w, long nx)
901 : {
902 32767 : long j, ly = n+3;
903 32767 : GEN y = cgetg(ly, t_POL);
904 32767 : y[1] = evalsigne(1);
905 401015 : for (j=2; j<ly; j++)
906 : {
907 : long k;
908 368248 : GEN a = cgetg(nx+2,t_POL);
909 368248 : a[1] = evalsigne(1) | evalvarn(w);
910 2028272 : for (k=0; k<nx; k++)
911 : {
912 1660024 : GEN xk = gel(x,k);
913 1660024 : if (typ(xk)==t_POL && varn(xk)==w)
914 1565649 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
915 : else
916 94375 : gel(a,k+2) = j==2 ? xk: gen_0;
917 : }
918 368248 : gel(y,j) = normalizepol_lg(a, nx+2);
919 : }
920 32767 : return normalizepol_lg(y,ly);
921 : }
922 :
923 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
924 : GEN
925 2051 : RgXY_swap(GEN x, long n, long w)
926 : {
927 2051 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
928 2051 : setvarn(z, varn(x)); return z;
929 : }
930 :
931 : long
932 1387 : RgXY_degreex(GEN b)
933 : {
934 1387 : long deg = 0, i;
935 1387 : if (!signe(b)) return -1;
936 15592 : for (i = 2; i < lg(b); ++i)
937 : {
938 14205 : GEN bi = gel(b, i);
939 14205 : if (typ(bi) == t_POL)
940 13742 : deg = maxss(deg, degpol(bi));
941 : }
942 1387 : return deg;
943 : }
944 :
945 : GEN
946 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
947 :
948 : /* return (x % X^n). Shallow */
949 : GEN
950 7944007 : RgXn_red_shallow(GEN a, long n)
951 : {
952 7944007 : long i, L = n+2, l = lg(a);
953 : GEN b;
954 7944007 : if (L >= l) return a; /* deg(x) < n */
955 5751644 : b = cgetg(L, t_POL); b[1] = a[1];
956 37327382 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
957 5751644 : return normalizepol_lg(b,L);
958 : }
959 :
960 : GEN
961 483 : RgXnV_red_shallow(GEN x, long n)
962 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
963 :
964 : /* return (x * X^n). Shallow */
965 : GEN
966 177405247 : RgX_shift_shallow(GEN a, long n)
967 : {
968 177405247 : long i, l = lg(a);
969 : GEN b;
970 177405247 : if (l == 2 || !n) return a;
971 109195485 : l += n;
972 109195485 : if (n < 0)
973 : {
974 54216095 : if (l <= 2) return pol_0(varn(a));
975 52741984 : b = cgetg(l, t_POL); b[1] = a[1];
976 52742822 : a -= n;
977 164032541 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
978 : } else {
979 54979390 : b = cgetg(l, t_POL); b[1] = a[1];
980 54983938 : a -= n; n += 2;
981 119507833 : for (i=2; i<n; i++) gel(b,i) = gen_0;
982 212056510 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
983 : }
984 107726760 : return b;
985 : }
986 : /* return (x * X^n). */
987 : GEN
988 1578036 : RgX_shift(GEN a, long n)
989 : {
990 1578036 : long i, l = lg(a);
991 : GEN b;
992 1578036 : if (l == 2 || !n) return RgX_copy(a);
993 1577084 : l += n;
994 1577084 : if (n < 0)
995 : {
996 1442 : if (l <= 2) return pol_0(varn(a));
997 1372 : b = cgetg(l, t_POL); b[1] = a[1];
998 1372 : a -= n;
999 9534 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1000 : } else {
1001 1575642 : b = cgetg(l, t_POL); b[1] = a[1];
1002 1575642 : a -= n; n += 2;
1003 3979090 : for (i=2; i<n; i++) gel(b,i) = gen_0;
1004 4346111 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1005 : }
1006 1577014 : return b;
1007 : }
1008 :
1009 : GEN
1010 317037 : RgX_rotate_shallow(GEN P, long k, long p)
1011 : {
1012 317037 : long i, l = lgpol(P);
1013 : GEN r;
1014 317037 : if (signe(P)==0)
1015 1365 : return pol_0(varn(P));
1016 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
1017 2100644 : for(i=0; i<p; i++)
1018 : {
1019 1784972 : long s = 2+(i+k)%p;
1020 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
1021 : }
1022 315672 : return RgX_renormalize(r);
1023 : }
1024 :
1025 : GEN
1026 2997279 : RgX_mulXn(GEN x, long d)
1027 : {
1028 : pari_sp av;
1029 : GEN z;
1030 : long v;
1031 2997279 : if (d >= 0) return RgX_shift(x, d);
1032 1472517 : d = -d;
1033 1472517 : v = RgX_val(x);
1034 1472517 : if (v >= d) return RgX_shift(x, -d);
1035 1472503 : av = avma;
1036 1472503 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1037 1472503 : return gc_upto(av, z);
1038 : }
1039 :
1040 : long
1041 588 : RgXV_maxdegree(GEN x)
1042 : {
1043 588 : long d = -1, i, l = lg(x);
1044 4494 : for (i = 1; i < l; i++)
1045 3906 : d = maxss(d, degpol(gel(x,i)));
1046 588 : return d;
1047 : }
1048 :
1049 : long
1050 3667905 : RgX_val(GEN x)
1051 : {
1052 3667905 : long i, lx = lg(x);
1053 3667905 : if (lx == 2) return LONG_MAX;
1054 4619737 : for (i = 2; i < lx; i++)
1055 4619681 : if (!isexactzero(gel(x,i))) break;
1056 3656978 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1057 3656922 : return i - 2;
1058 : }
1059 : long
1060 80796199 : RgX_valrem(GEN x, GEN *Z)
1061 : {
1062 80796199 : long v, i, lx = lg(x);
1063 80796199 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1064 125213453 : for (i = 2; i < lx; i++)
1065 125214719 : if (!isexactzero(gel(x,i))) break;
1066 : /* possible with nonrational zeros */
1067 80796361 : if (i == lx)
1068 : {
1069 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1070 14 : return LONG_MAX;
1071 : }
1072 80796347 : v = i - 2;
1073 80796347 : *Z = RgX_shift_shallow(x, -v);
1074 80808436 : return v;
1075 : }
1076 : long
1077 852252 : RgX_valrem_inexact(GEN x, GEN *Z)
1078 : {
1079 : long v;
1080 852252 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1081 869239 : for (v = 0;; v++)
1082 869239 : if (!gequal0(gel(x,2+v))) break;
1083 852245 : if (Z) *Z = RgX_shift_shallow(x, -v);
1084 852245 : return v;
1085 : }
1086 :
1087 : GEN
1088 68187 : RgXQC_red(GEN x, GEN T)
1089 414505 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1090 :
1091 : GEN
1092 1211 : RgXQV_red(GEN x, GEN T)
1093 27832 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1094 :
1095 : GEN
1096 13195 : RgXQM_red(GEN x, GEN T)
1097 81382 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1098 :
1099 : GEN
1100 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1101 : {
1102 322 : return RgXQM_red(RgM_mul(P, Q), T);
1103 : }
1104 :
1105 : GEN
1106 508979 : RgXQX_red(GEN P, GEN T)
1107 : {
1108 508979 : long i, l = lg(P);
1109 508979 : GEN Q = cgetg(l, t_POL);
1110 508979 : Q[1] = P[1];
1111 2628994 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1112 508979 : return normalizepol_lg(Q, l);
1113 : }
1114 :
1115 : GEN
1116 830470 : RgX_deriv(GEN x)
1117 : {
1118 830470 : long i,lx = lg(x)-1;
1119 : GEN y;
1120 :
1121 830470 : if (lx<3) return pol_0(varn(x));
1122 827369 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1123 3652459 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1124 827364 : y[1] = x[1]; return normalizepol_lg(y,i);
1125 : }
1126 :
1127 : GEN
1128 2584095 : RgX_recipspec_shallow(GEN x, long l, long n)
1129 : {
1130 : long i;
1131 2584095 : GEN z = cgetg(n+2,t_POL);
1132 2584103 : z[1] = 0; z += 2;
1133 144332725 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1134 2833838 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1135 2584103 : return normalizepol_lg(z-2,n+2);
1136 : }
1137 :
1138 : GEN
1139 643038 : RgXn_recip_shallow(GEN P, long n)
1140 : {
1141 643038 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1142 643044 : setvarn(Q, varn(P));
1143 643044 : return Q;
1144 : }
1145 :
1146 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1147 : GEN
1148 31388 : RgX_recip(GEN x)
1149 : {
1150 : long lx, i, j;
1151 31388 : GEN y = cgetg_copy(x, &lx);
1152 274519 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1153 31388 : return normalizepol_lg(y,lx);
1154 : }
1155 : /* shallow version */
1156 : GEN
1157 59388 : RgX_recip_shallow(GEN x)
1158 : {
1159 : long lx, i, j;
1160 59388 : GEN y = cgetg_copy(x, &lx);
1161 356118 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1162 59388 : return normalizepol_lg(y,lx);
1163 : }
1164 :
1165 : GEN
1166 5206430 : RgX_recip_i(GEN x)
1167 : {
1168 : long lx, i, j;
1169 5206430 : GEN y = cgetg_copy(x, &lx);
1170 26412289 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1171 5206429 : return y;
1172 : }
1173 : /*******************************************************************/
1174 : /* */
1175 : /* ADDITION / SUBTRACTION */
1176 : /* */
1177 : /*******************************************************************/
1178 : /* same variable */
1179 : GEN
1180 145709193 : RgX_add(GEN x, GEN y)
1181 : {
1182 145709193 : long i, lx = lg(x), ly = lg(y);
1183 : GEN z;
1184 145709193 : if (ly <= lx) {
1185 130999185 : z = cgetg(lx,t_POL); z[1] = x[1];
1186 509829858 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1187 190852403 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1188 130973603 : z = normalizepol_lg(z, lx);
1189 : } else {
1190 14710008 : z = cgetg(ly,t_POL); z[1] = y[1];
1191 55444308 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1192 42728436 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1193 14711581 : z = normalizepol_lg(z, ly);
1194 : }
1195 145699735 : return z;
1196 : }
1197 : GEN
1198 72827995 : RgX_sub(GEN x, GEN y)
1199 : {
1200 72827995 : long i, lx = lg(x), ly = lg(y);
1201 : GEN z;
1202 72827995 : if (ly <= lx) {
1203 38986733 : z = cgetg(lx,t_POL); z[1] = x[1];
1204 188371718 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1205 66482697 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1206 38963366 : z = normalizepol_lg(z, lx);
1207 : } else {
1208 33841262 : z = cgetg(ly,t_POL); z[1] = y[1];
1209 125743811 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1210 74956873 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1211 33808529 : z = normalizepol_lg(z, ly);
1212 : }
1213 72825565 : return z;
1214 : }
1215 : GEN
1216 9247367 : RgX_neg(GEN x)
1217 54563977 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1218 :
1219 : GEN
1220 45363744 : RgX_Rg_add(GEN y, GEN x)
1221 : {
1222 : GEN z;
1223 45363744 : long lz = lg(y), i;
1224 45363744 : if (lz == 2) return scalarpol(x,varn(y));
1225 23611204 : z = cgetg(lz,t_POL); z[1] = y[1];
1226 23611183 : gel(z,2) = gadd(gel(y,2),x);
1227 80862943 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1228 : /* probably useless unless lz = 3, but cannot be skipped if y is
1229 : * an inexact 0 */
1230 23611210 : return normalizepol_lg(z,lz);
1231 : }
1232 : GEN
1233 65102 : RgX_Rg_add_shallow(GEN y, GEN x)
1234 : {
1235 : GEN z;
1236 65102 : long lz = lg(y), i;
1237 65102 : if (lz == 2) return scalarpol(x,varn(y));
1238 65102 : z = cgetg(lz,t_POL); z[1] = y[1];
1239 65102 : gel(z,2) = gadd(gel(y,2),x);
1240 130340 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1241 65102 : return normalizepol_lg(z,lz);
1242 : }
1243 : GEN
1244 1226042 : RgX_Rg_sub(GEN y, GEN x)
1245 : {
1246 : GEN z;
1247 1226042 : long lz = lg(y), i;
1248 1226042 : if (lz == 2)
1249 : { /* scalarpol(gneg(x),varn(y)) optimized */
1250 922867 : long v = varn(y);
1251 922867 : if (isrationalzero(x)) return pol_0(v);
1252 2235 : z = cgetg(3,t_POL);
1253 2235 : z[1] = gequal0(x)? evalvarn(v)
1254 2235 : : evalvarn(v) | evalsigne(1);
1255 2235 : gel(z,2) = gneg(x); return z;
1256 : }
1257 303175 : z = cgetg(lz,t_POL); z[1] = y[1];
1258 303176 : gel(z,2) = gsub(gel(y,2),x);
1259 697262 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1260 303178 : return normalizepol_lg(z,lz);
1261 : }
1262 : GEN
1263 4953322 : Rg_RgX_sub(GEN x, GEN y)
1264 : {
1265 : GEN z;
1266 4953322 : long lz = lg(y), i;
1267 4953322 : if (lz == 2) return scalarpol(x,varn(y));
1268 4882928 : z = cgetg(lz,t_POL); z[1] = y[1];
1269 4882914 : gel(z,2) = gsub(x, gel(y,2));
1270 7336442 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1271 4882890 : return normalizepol_lg(z,lz);
1272 : }
1273 : /*******************************************************************/
1274 : /* */
1275 : /* KARATSUBA MULTIPLICATION */
1276 : /* */
1277 : /*******************************************************************/
1278 : #if 0
1279 : /* to debug Karatsuba-like routines */
1280 : GEN
1281 : zx_debug_spec(GEN x, long nx)
1282 : {
1283 : GEN z = cgetg(nx+2,t_POL);
1284 : long i;
1285 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1286 : z[1] = evalsigne(1); return z;
1287 : }
1288 :
1289 : GEN
1290 : RgX_debug_spec(GEN x, long nx)
1291 : {
1292 : GEN z = cgetg(nx+2,t_POL);
1293 : long i;
1294 : for (i=0; i<nx; i++) z[i+2] = x[i];
1295 : z[1] = evalsigne(1); return z;
1296 : }
1297 : #endif
1298 :
1299 : /* generic multiplication */
1300 : GEN
1301 8939618 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1302 : {
1303 : GEN z, t;
1304 : long i;
1305 8939618 : if (nx == ny) {
1306 1567387 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1307 5015633 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1308 1567376 : return normalizepol_lg(z, nx+2);
1309 : }
1310 7372231 : if (ny < nx) {
1311 7187056 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1312 26454559 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1313 17173643 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1314 7186572 : return normalizepol_lg(z, nx+2);
1315 : } else {
1316 185175 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1317 3671822 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1318 447113 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1319 185192 : return normalizepol_lg(z, ny+2);
1320 : }
1321 : }
1322 : GEN
1323 222294 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1324 : {
1325 : GEN z, t;
1326 : long i;
1327 222294 : if (nx == ny) {
1328 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1329 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1330 12824 : return normalizepol_lg(z, nx+2);
1331 : }
1332 209470 : if (ny < nx) {
1333 207685 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1334 3724613 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1335 2371419 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1336 207685 : return normalizepol_lg(z, nx+2);
1337 : } else {
1338 1785 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1339 331478 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1340 12236 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1341 1785 : return normalizepol_lg(z, ny+2);
1342 : }
1343 : }
1344 :
1345 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1346 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1347 : * t_VECSMALL, to hold this, which can be left on stack: GC functions
1348 : * will not crash on it. The returned vector itself is not a proper GEN,
1349 : * we access the coefficients as x[i], i = 0..deg(x) */
1350 : static GEN
1351 81702183 : RgXspec_kill0(GEN x, long lx)
1352 : {
1353 81702183 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit GC-wise */
1354 : long i;
1355 246598954 : for (i=0; i <lx; i++)
1356 : {
1357 164896717 : GEN c = gel(x,i);
1358 164896717 : z[i] = (long)(isrationalzero(c)? NULL: c);
1359 : }
1360 81702237 : return z;
1361 : }
1362 :
1363 : INLINE GEN
1364 111071815 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1365 : {
1366 111071815 : pari_sp av = avma;
1367 111071815 : GEN s = NULL;
1368 : long i;
1369 :
1370 391695955 : for (i=a; i<b; i++)
1371 280633492 : if (gel(y,i) && gel(x,-i))
1372 : {
1373 202308803 : GEN t = gmul(gel(y,i), gel(x,-i));
1374 202304555 : s = s? gadd(s, t): t;
1375 : }
1376 111062463 : return s? gc_upto(av, s): gen_0;
1377 : }
1378 :
1379 : /* assume nx >= ny > 0, return x * y * t^v */
1380 : static GEN
1381 34204196 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1382 : {
1383 : long i, lz, nz;
1384 : GEN z;
1385 :
1386 34204196 : x = RgXspec_kill0(x,nx);
1387 34204204 : y = RgXspec_kill0(y,ny);
1388 34204213 : lz = nx + ny + 1; nz = lz-2;
1389 34204213 : lz += v;
1390 34204213 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1391 72781359 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1392 83204465 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1393 56638680 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1394 49000357 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1395 34203614 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1396 : }
1397 :
1398 : /* return (x * X^d) + y. Assume d > 0 */
1399 : GEN
1400 10131643 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1401 : {
1402 : GEN x, y, xd, yd, zd;
1403 : long a, lz, nx, ny;
1404 :
1405 10131643 : if (!signe(x0)) return y0;
1406 9463414 : ny = lgpol(y0);
1407 9463408 : nx = lgpol(x0);
1408 9463482 : zd = (GEN)avma;
1409 9463482 : x = x0 + 2; y = y0 + 2; a = ny-d;
1410 9463482 : if (a <= 0)
1411 : {
1412 1479607 : lz = nx+d+2;
1413 1479607 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1414 3397754 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1415 1479609 : x = zd + a;
1416 1497306 : while (zd > x) gel(--zd,0) = gen_0;
1417 : }
1418 : else
1419 : {
1420 7983875 : xd = new_chunk(d); yd = y+d;
1421 7983874 : x = RgX_addspec_shallow(x,yd, nx,a);
1422 7983819 : lz = (a>nx)? ny+2: lg(x)+d;
1423 39690281 : x += 2; while (xd > x) *--zd = *--xd;
1424 : }
1425 23038282 : while (yd > y) *--zd = *--yd;
1426 9463428 : *--zd = x0[1];
1427 9463428 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1428 : }
1429 : GEN
1430 514764 : RgX_addmulXn(GEN x0, GEN y0, long d)
1431 : {
1432 : GEN x, y, xd, yd, zd;
1433 : long a, lz, nx, ny;
1434 :
1435 514764 : if (!signe(x0)) return RgX_copy(y0);
1436 513980 : nx = lgpol(x0);
1437 513980 : ny = lgpol(y0);
1438 513980 : zd = (GEN)avma;
1439 513980 : x = x0 + 2; y = y0 + 2; a = ny-d;
1440 513980 : if (a <= 0)
1441 : {
1442 291686 : lz = nx+d+2;
1443 291686 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1444 4293251 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1445 291686 : x = zd + a;
1446 757796 : while (zd > x) gel(--zd,0) = gen_0;
1447 : }
1448 : else
1449 : {
1450 222294 : xd = new_chunk(d); yd = y+d;
1451 222294 : x = RgX_addspec(x,yd, nx,a);
1452 222294 : lz = (a>nx)? ny+2: lg(x)+d;
1453 8416054 : x += 2; while (xd > x) *--zd = *--xd;
1454 : }
1455 2614912 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1456 513980 : *--zd = x0[1];
1457 513980 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1458 : }
1459 :
1460 : /* return x * y mod t^n */
1461 : static GEN
1462 6627987 : RgXn_mul_basecase(GEN x, GEN y, long n)
1463 : {
1464 6627987 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1465 : GEN z;
1466 6627987 : if (lx < 0) return pol_0(varn(x));
1467 6627987 : if (ly < 0) return pol_0(varn(x));
1468 6627987 : z = cgetg(lz, t_POL) + 2;
1469 6627987 : x+=2; if (lx > n) lx = n;
1470 6627987 : y+=2; if (ly > n) ly = n;
1471 6627987 : z[-1] = x[-1];
1472 6627987 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1473 6627987 : x = RgXspec_kill0(x, lx);
1474 6627987 : y = RgXspec_kill0(y, ly);
1475 : /* x:y:z [i] = term of degree i */
1476 26223941 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1477 11826583 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1478 6674059 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1479 6627987 : return normalizepol_lg(z - 2, lz);
1480 : }
1481 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1482 : static GEN
1483 9414549 : RgXn_mul2(GEN f, GEN g, long n)
1484 : {
1485 9414549 : pari_sp av = avma;
1486 : GEN fe,fo, ge,go, l,h,m;
1487 : long n0, n1;
1488 9414549 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1489 6663484 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1490 35497 : n0 = n>>1; n1 = n-n0;
1491 35497 : RgX_even_odd(f, &fe, &fo);
1492 35497 : RgX_even_odd(g, &ge, &go);
1493 35497 : l = RgXn_mul2(fe,ge,n1);
1494 35497 : h = RgXn_mul2(fo,go,n0);
1495 35497 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1496 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1497 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1498 35497 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1499 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1500 35497 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1501 35497 : m = RgX_inflate(m,2);
1502 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1503 35497 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1504 35497 : h = RgX_inflate(h,2);
1505 35497 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1506 35497 : return gc_upto(av, h);
1507 : }
1508 : /* (f*g) \/ x^n */
1509 : static GEN
1510 1589510 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1511 : {
1512 1589510 : long d = degpol(f)+degpol(g) + 1 - n;
1513 : GEN h;
1514 1589510 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1515 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1516 : RgX_recip_i(g), d));
1517 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1518 : }
1519 :
1520 : /* (f*g) \/ x^n */
1521 : static GEN
1522 0 : RgX_sqrhigh_i2(GEN f, long n)
1523 : {
1524 0 : long d = 2*degpol(f)+ 1 - n;
1525 : GEN h;
1526 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1527 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1528 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1529 : }
1530 :
1531 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1532 : * b+2 were sent instead. na, nb = number of terms of a, b.
1533 : * Only c, c0, c1, c2 are genuine GEN.
1534 : */
1535 : GEN
1536 40270428 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1537 : {
1538 : GEN a0, c, c0;
1539 40270428 : long n0, n0a, i, v = 0;
1540 : pari_sp av;
1541 :
1542 65268736 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1543 61961141 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1544 40270395 : if (na < nb) swapspec(a,b, na,nb);
1545 40270395 : if (!nb) return pol_0(0);
1546 :
1547 34683391 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1548 479212 : RgX_shift_inplace_init(v);
1549 479228 : i = (na>>1); n0 = na-i; na = i;
1550 479228 : av = avma; a0 = a+n0; n0a = n0;
1551 1334087 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1552 :
1553 479228 : if (nb > n0)
1554 : {
1555 : GEN b0,c1,c2;
1556 : long n0b;
1557 :
1558 477869 : nb -= n0; b0 = b+n0; n0b = n0;
1559 1444999 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1560 477869 : c = RgX_mulspec(a,b,n0a,n0b);
1561 477869 : c0 = RgX_mulspec(a0,b0, na,nb);
1562 :
1563 477869 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1564 477869 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1565 :
1566 477869 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1567 477869 : c2 = RgX_sub(c1, RgX_add(c0,c));
1568 477869 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1569 : }
1570 : else
1571 : {
1572 1359 : c = RgX_mulspec(a,b,n0a,nb);
1573 1359 : c0 = RgX_mulspec(a0,b,na,nb);
1574 : }
1575 479228 : c0 = RgX_addmulXn(c0,c,n0);
1576 479228 : return RgX_shift_inplace(gc_upto(av,c0), v);
1577 : }
1578 :
1579 : INLINE GEN
1580 100082 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1581 : {
1582 100082 : pari_sp av = avma;
1583 100082 : GEN s = NULL;
1584 100082 : long j, l = (i+1)>>1;
1585 203924 : for (j=a; j<l; j++)
1586 : {
1587 103842 : GEN xj = gel(x,j), xx = gel(x,i-j);
1588 103842 : if (xj && xx)
1589 : {
1590 94707 : GEN t = gmul(xj, xx);
1591 94707 : s = s? gadd(s, t): t;
1592 : }
1593 : }
1594 100082 : if (s) s = gshift(s,1);
1595 100082 : if ((i&1) == 0)
1596 : {
1597 69041 : GEN t = gel(x, i>>1);
1598 69041 : if (t) {
1599 65982 : t = gsqr(t);
1600 65982 : s = s? gadd(s, t): t;
1601 : }
1602 : }
1603 100082 : return s? gc_upto(av,s): gen_0;
1604 : }
1605 : static GEN
1606 38000 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1607 : {
1608 : long i, lz, nz;
1609 : GEN z;
1610 :
1611 38000 : if (!nx) return pol_0(0);
1612 38000 : x = RgXspec_kill0(x,nx);
1613 38000 : lz = (nx << 1) + 1, nz = lz-2;
1614 38000 : lz += v;
1615 38000 : z = cgetg(lz,t_POL) + 2;
1616 94602 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1617 107041 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1618 69041 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1619 38000 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1620 : }
1621 : /* return x^2 mod t^n */
1622 : static GEN
1623 0 : RgXn_sqr_basecase(GEN x, long n)
1624 : {
1625 0 : long i, lz = n+2, lx = lgpol(x);
1626 : GEN z;
1627 0 : if (lx < 0) return pol_0(varn(x));
1628 0 : z = cgetg(lz, t_POL);
1629 0 : z[1] = x[1];
1630 0 : x+=2; if (lx > n) lx = n;
1631 0 : x = RgXspec_kill0(x,lx);
1632 0 : z+=2;/* x:z [i] = term of degree i */
1633 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1634 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1635 0 : z -= 2; return normalizepol_lg(z, lz);
1636 : }
1637 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1638 : static GEN
1639 315 : RgXn_sqr2(GEN f, long n)
1640 : {
1641 315 : pari_sp av = avma;
1642 : GEN fe,fo, l,h,m;
1643 : long n0, n1;
1644 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1645 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1646 0 : n0 = n>>1; n1 = n-n0;
1647 0 : RgX_even_odd(f, &fe, &fo);
1648 0 : l = RgXn_sqr(fe,n1);
1649 0 : h = RgXn_sqr(fo,n0);
1650 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1651 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1652 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1653 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1654 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1655 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1656 0 : m = RgX_inflate(m,2);
1657 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1658 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1659 0 : h = RgX_inflate(h,2);
1660 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1661 0 : return gc_upto(av, h);
1662 : }
1663 : GEN
1664 38039 : RgX_sqrspec(GEN a, long na)
1665 : {
1666 : GEN a0, c, c0, c1;
1667 38039 : long n0, n0a, i, v = 0;
1668 : pari_sp av;
1669 :
1670 66340 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1671 38039 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1672 39 : RgX_shift_inplace_init(v);
1673 39 : i = (na>>1); n0 = na-i; na = i;
1674 39 : av = avma; a0 = a+n0; n0a = n0;
1675 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1676 :
1677 39 : c = RgX_sqrspec(a,n0a);
1678 39 : c0 = RgX_sqrspec(a0,na);
1679 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1680 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1681 39 : c0 = RgX_addmulXn(c0,c,n0);
1682 39 : return RgX_shift_inplace(gc_upto(av,c0), v);
1683 : }
1684 :
1685 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1686 : GEN
1687 1719207 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1688 : {
1689 1719207 : GEN z = RgX_mul(A, B);
1690 1719199 : if (a < b)
1691 10428 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1692 1708771 : else if (a > b)
1693 1104945 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1694 : else
1695 603826 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1696 1719201 : return z;
1697 : }
1698 :
1699 : GEN
1700 38834121 : RgX_mul_i(GEN x, GEN y)
1701 : {
1702 38834121 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1703 38833789 : setvarn(z, varn(x)); return z;
1704 : }
1705 :
1706 : GEN
1707 37961 : RgX_sqr_i(GEN x)
1708 : {
1709 37961 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1710 37961 : setvarn(z,varn(x)); return z;
1711 : }
1712 :
1713 : /*******************************************************************/
1714 : /* */
1715 : /* DIVISION */
1716 : /* */
1717 : /*******************************************************************/
1718 : GEN
1719 7483412 : RgX_Rg_divexact(GEN x, GEN y) {
1720 7483412 : long i, lx = lg(x);
1721 : GEN z;
1722 7483412 : if (lx == 2) return gcopy(x);
1723 7432447 : switch(typ(y))
1724 : {
1725 7270294 : case t_INT:
1726 7270294 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1727 6887103 : break;
1728 5243 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1729 : }
1730 7044013 : z = cgetg(lx, t_POL); z[1] = x[1];
1731 35452303 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1732 7042837 : return z;
1733 : }
1734 : GEN
1735 30719529 : RgX_Rg_div(GEN x, GEN y) {
1736 30719529 : long i, lx = lg(x);
1737 : GEN z;
1738 30719529 : if (lx == 2) return gcopy(x);
1739 30445728 : switch(typ(y))
1740 : {
1741 21212923 : case t_INT:
1742 21212923 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1743 4270041 : break;
1744 6069 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1745 : }
1746 13496777 : z = cgetg(lx, t_POL); z[1] = x[1];
1747 50360525 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1748 13496372 : return normalizepol_lg(z, lx);
1749 : }
1750 : GEN
1751 39130 : RgX_normalize(GEN x)
1752 : {
1753 39130 : GEN z, d = NULL;
1754 39130 : long i, n = lg(x)-1;
1755 39130 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1756 39130 : if (i == 1) return pol_0(varn(x));
1757 39130 : if (i == n && isint1(d)) return x;
1758 17612 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1759 31612 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1760 17612 : gel(z,n) = Rg_get_1(d); return z;
1761 : }
1762 : GEN
1763 10458 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1764 : GEN
1765 265423 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1766 : {
1767 265423 : long l = lg(a), i;
1768 : GEN a0, z0, z;
1769 :
1770 265423 : if (l <= 3)
1771 : {
1772 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1773 0 : return pol_0(varn(a));
1774 : }
1775 265423 : z = cgetg(l-1, t_POL);
1776 265423 : z[1] = a[1];
1777 265423 : a0 = a + l-1;
1778 265423 : z0 = z + l-2; *z0 = *a0--;
1779 3024610 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1780 : {
1781 2759188 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1782 2759187 : gel(z0,0) = t;
1783 : }
1784 265422 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1785 265422 : return z;
1786 : }
1787 : /* Polynomial division x / y:
1788 : * if pr = ONLY_REM return remainder, otherwise return quotient
1789 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1790 : * if pr != NULL set *pr to remainder, as the last object on stack */
1791 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1792 : static GEN
1793 23519675 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1794 : {
1795 : pari_sp avy, av, av1;
1796 : long dx,dy,dz,i,j,sx,lr;
1797 : GEN z,p1,p2,rem,y_lead,mod,p;
1798 : GEN (*f)(GEN,GEN);
1799 :
1800 23519675 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1801 :
1802 23519675 : dy = degpol(y);
1803 23519650 : y_lead = gel(y,dy+2);
1804 23519650 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1805 : {
1806 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1807 0 : for (dy--; dy>=0; dy--)
1808 : {
1809 0 : y_lead = gel(y,dy+2);
1810 0 : if (!gequal0(y_lead)) break;
1811 : }
1812 : }
1813 23519654 : if (!dy) /* y is constant */
1814 : {
1815 6844 : if (pr == ONLY_REM) return pol_0(varn(x));
1816 6844 : z = RgX_Rg_div(x, y_lead);
1817 6844 : if (pr == ONLY_DIVIDES) return z;
1818 2630 : if (pr) *pr = pol_0(varn(x));
1819 2630 : return z;
1820 : }
1821 23512810 : dx = degpol(x);
1822 23512776 : if (dx < dy)
1823 : {
1824 3847700 : if (pr == ONLY_REM) return RgX_copy(x);
1825 348591 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1826 348570 : z = pol_0(varn(x));
1827 348570 : if (pr) *pr = RgX_copy(x);
1828 348570 : return z;
1829 : }
1830 :
1831 : /* x,y in R[X], y non constant */
1832 19665076 : av = avma;
1833 19665076 : p = NULL;
1834 19665076 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1835 : {
1836 130277 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1837 130277 : if (!z) return gc_NULL(av);
1838 130277 : z = FpX_to_mod(z, p);
1839 130277 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1840 71134 : return gc_upto(av, z);
1841 59143 : *pr = FpX_to_mod(*pr, p);
1842 59143 : return gc_all(av, 2, &z, pr);
1843 : }
1844 19534949 : switch(typ(y_lead))
1845 : {
1846 0 : case t_REAL:
1847 0 : y_lead = ginv(y_lead);
1848 0 : f = gmul; mod = NULL;
1849 0 : break;
1850 1878 : case t_INTMOD:
1851 1878 : case t_POLMOD: y_lead = ginv(y_lead);
1852 1878 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1853 1878 : break;
1854 19533071 : default: if (gequal1(y_lead)) y_lead = NULL;
1855 19533056 : f = gdiv; mod = NULL;
1856 : }
1857 :
1858 19534934 : if (y_lead == NULL)
1859 17521994 : p2 = gel(x,dx+2);
1860 : else {
1861 : for(;;) {
1862 2012940 : p2 = f(gel(x,dx+2),y_lead);
1863 2012942 : p2 = simplify_shallow(p2);
1864 2012942 : if (!isexactzero(p2) || (--dx < 0)) break;
1865 : }
1866 2012942 : if (dx < dy) /* leading coeff of x was in fact zero */
1867 : {
1868 0 : if (pr == ONLY_DIVIDES) {
1869 0 : set_avma(av);
1870 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1871 : }
1872 0 : if (pr == ONLY_REM)
1873 : {
1874 0 : if (dx < 0)
1875 0 : return gc_GEN(av, scalarpol(p2, varn(x)));
1876 : else
1877 : {
1878 : GEN t;
1879 0 : set_avma(av);
1880 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1881 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1882 0 : return t;
1883 : }
1884 : }
1885 0 : if (pr) /* cf ONLY_REM above */
1886 : {
1887 0 : if (dx < 0)
1888 : {
1889 0 : p2 = gclone(p2);
1890 0 : set_avma(av);
1891 0 : z = pol_0(varn(x));
1892 0 : x = scalarpol(p2, varn(x));
1893 0 : gunclone(p2);
1894 : }
1895 : else
1896 : {
1897 : GEN t;
1898 0 : set_avma(av);
1899 0 : z = pol_0(varn(x));
1900 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1901 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1902 0 : x = t;
1903 : }
1904 0 : *pr = x;
1905 : }
1906 : else
1907 : {
1908 0 : set_avma(av);
1909 0 : z = pol_0(varn(x));
1910 : }
1911 0 : return z;
1912 : }
1913 : }
1914 : /* dx >= dy */
1915 19534936 : avy = avma;
1916 19534936 : dz = dx-dy;
1917 19534936 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1918 19534923 : x += 2;
1919 19534923 : z += 2;
1920 19534923 : y += 2;
1921 19534923 : gel(z,dz) = gcopy(p2);
1922 :
1923 54100326 : for (i=dx-1; i>=dy; i--)
1924 : {
1925 34565577 : av1=avma; p1=gel(x,i);
1926 1140275538 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1927 34560422 : if (y_lead) p1 = simplify(f(p1,y_lead));
1928 :
1929 34560424 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1930 : else
1931 24571410 : p1 = avma==av1? gcopy(p1): gc_upto(av1,p1);
1932 34564928 : gel(z,i-dy) = p1;
1933 : }
1934 19534749 : if (!pr) return gc_upto(av,z-2);
1935 :
1936 12556418 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1937 12913540 : for (sx=0; ; i--)
1938 : {
1939 12913540 : p1 = gel(x,i);
1940 : /* we always enter this loop at least once */
1941 32235613 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1942 12912287 : if (mod && avma==av1) p1 = gmul(p1,mod);
1943 12912287 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1944 3576696 : if (!isexactzero(p1)) break;
1945 3426583 : if (!i) break;
1946 357066 : set_avma(av1);
1947 : }
1948 12555615 : if (pr == ONLY_DIVIDES)
1949 : {
1950 7931 : if (sx) return gc_NULL(av);
1951 7910 : set_avma((pari_sp)rem); return gc_upto(av,z-2);
1952 : }
1953 12547684 : lr=i+3; rem -= lr;
1954 12547684 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1955 12505395 : else p1 = gc_upto((pari_sp)rem,p1);
1956 12548435 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1957 12548435 : rem[1] = z[-1];
1958 12548435 : rem += 2;
1959 12548435 : gel(rem,i) = p1;
1960 17295846 : for (i--; i>=0; i--)
1961 : {
1962 4747440 : av1=avma; p1 = gel(x,i);
1963 13653672 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1964 4746928 : if (mod && avma==av1) p1 = gmul(p1,mod);
1965 4747180 : gel(rem,i) = avma==av1? gcopy(p1):gc_upto(av1,p1);
1966 : }
1967 12548406 : rem -= 2;
1968 12548406 : if (!sx) (void)normalizepol_lg(rem, lr);
1969 12548660 : if (pr == ONLY_REM) return gc_upto(av,rem);
1970 6761956 : z -= 2; *pr = rem; return gc_all_unsafe(av,avy,2,&z,pr);
1971 : }
1972 :
1973 : GEN
1974 14233929 : RgX_divrem(GEN x, GEN y, GEN *pr)
1975 : {
1976 14233929 : if (pr == ONLY_REM) return RgX_rem(x, y);
1977 14233929 : return RgX_divrem_i(x, y, pr);
1978 : }
1979 :
1980 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1981 : GEN
1982 156871 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1983 : {
1984 156871 : long vx = varn(x), dx = degpol(x), dy = degpol(y), dz, i, j, sx, lr;
1985 : pari_sp av0, av;
1986 : GEN z, p1, rem, lead;
1987 :
1988 156871 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1989 156871 : if (dx < dy)
1990 : {
1991 41644 : if (pr)
1992 : {
1993 41616 : av0 = avma; x = RgXQX_red(x, T);
1994 41616 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1995 41609 : if (pr == ONLY_REM) return x;
1996 0 : *pr = x;
1997 : }
1998 28 : return pol_0(vx);
1999 : }
2000 115227 : lead = leading_coeff(y);
2001 115227 : if (!dy) /* y is constant */
2002 : {
2003 602 : if (pr && pr != ONLY_DIVIDES)
2004 : {
2005 0 : if (pr == ONLY_REM) return pol_0(vx);
2006 0 : *pr = pol_0(vx);
2007 : }
2008 602 : if (gequal1(lead)) return RgX_copy(x);
2009 0 : av0 = avma; x = gmul(x, ginvmod(lead,T));
2010 0 : return gc_upto(av0, RgXQX_red(x,T));
2011 : }
2012 114625 : av0 = avma; dz = dx-dy;
2013 114625 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
2014 114625 : set_avma(av0);
2015 114625 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2016 114625 : x += 2; y += 2; z += 2;
2017 :
2018 114625 : p1 = gel(x,dx); av = avma;
2019 114625 : gel(z,dz) = lead? gc_upto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2020 581855 : for (i=dx-1; i>=dy; i--)
2021 : {
2022 467230 : av = avma; p1 = gel(x,i);
2023 2418735 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2024 467205 : if (lead) p1 = gmul(grem(p1, T), lead);
2025 467206 : gel(z,i-dy) = gc_upto(av, grem(p1, T));
2026 : }
2027 114625 : if (!pr) { guncloneNULL(lead); return z-2; }
2028 :
2029 113351 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2030 191659 : for (sx=0; ; i--)
2031 : {
2032 191659 : p1 = gel(x,i);
2033 663112 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2034 191659 : p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2035 105922 : if (!i) break;
2036 78308 : set_avma(av);
2037 : }
2038 113350 : if (pr == ONLY_DIVIDES)
2039 : {
2040 27142 : guncloneNULL(lead);
2041 27143 : if (sx) return gc_NULL(av0);
2042 24577 : return gc_const((pari_sp)rem, z-2);
2043 : }
2044 86208 : lr=i+3; rem -= lr; av = (pari_sp)rem;
2045 86208 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2046 86208 : rem[1] = z[-1];
2047 86208 : rem += 2; gel(rem,i) = gc_upto(av, p1);
2048 187624 : for (i--; i>=0; i--)
2049 : {
2050 101416 : av = avma; p1 = gel(x,i);
2051 337537 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2052 101416 : gel(rem,i) = gc_upto(av, grem(p1, T));
2053 : }
2054 86208 : rem -= 2;
2055 86208 : guncloneNULL(lead);
2056 86208 : if (!sx) (void)normalizepol_lg(rem, lr);
2057 86208 : if (pr == ONLY_REM) return gc_upto(av0,rem);
2058 168 : *pr = rem; return z-2;
2059 : }
2060 :
2061 : /*******************************************************************/
2062 : /* */
2063 : /* PSEUDO-DIVISION */
2064 : /* */
2065 : /*******************************************************************/
2066 : INLINE GEN
2067 4260634 : rem(GEN c, GEN T)
2068 : {
2069 4260634 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2070 4260636 : return c;
2071 : }
2072 :
2073 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2074 : int
2075 17804 : ZXQX_dvd(GEN x, GEN y, GEN T)
2076 : {
2077 : long dx, dy, i, T_ismonic;
2078 17804 : pari_sp av = avma, av2;
2079 : GEN y_lead;
2080 :
2081 17804 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2082 17804 : dy = degpol(y); y_lead = gel(y,dy+2);
2083 17804 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2084 : /* if monic, no point in using pseudo-division */
2085 17804 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2086 14808 : T_ismonic = gequal1(leading_coeff(T));
2087 14808 : dx = degpol(x);
2088 14808 : if (dx < dy) return !signe(x);
2089 14808 : (void)new_chunk(2);
2090 14808 : x = RgX_recip_i(x)+2;
2091 14808 : y = RgX_recip_i(y)+2;
2092 : /* pay attention to sparse divisors */
2093 30117 : for (i = 1; i <= dy; i++)
2094 15309 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2095 14808 : av2 = avma;
2096 : for (;;)
2097 73023 : {
2098 87831 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2099 87831 : x0 = gneg(x0);
2100 87831 : m = gcdii(cx, y0);
2101 87831 : if (!equali1(m))
2102 : {
2103 85329 : x0 = gdiv(x0, m);
2104 85329 : y0 = diviiexact(y0, m);
2105 85329 : if (equali1(y0)) y0 = NULL;
2106 : }
2107 179435 : for (i=1; i<=dy; i++)
2108 : {
2109 91604 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2110 91604 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2111 91604 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2112 91604 : gel(x,i) = c;
2113 : }
2114 688210 : for ( ; i<=dx; i++)
2115 : {
2116 600379 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2117 600379 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2118 600379 : gel(x,i) = c;
2119 : }
2120 102377 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2121 87831 : if (dx < dy) break;
2122 73023 : if (gc_needed(av2,1))
2123 : {
2124 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2125 28 : gc_slice(av2,x,dx+1);
2126 : }
2127 : }
2128 14808 : return gc_bool(av, dx < 0);
2129 : }
2130 :
2131 : /* T either NULL or a t_POL. */
2132 : GEN
2133 660119 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2134 : {
2135 660119 : long vx = varn(x), dx, dy, dz, i, lx, p;
2136 660119 : pari_sp av = avma, av2;
2137 : GEN y_lead;
2138 :
2139 660119 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2140 660119 : dy = degpol(y); y_lead = gel(y,dy+2);
2141 : /* if monic, no point in using pseudo-division */
2142 660119 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2143 610194 : dx = degpol(x);
2144 610194 : if (dx < dy) return RgX_copy(x);
2145 610187 : (void)new_chunk(2);
2146 610187 : x = RgX_recip_i(x)+2;
2147 610187 : y = RgX_recip_i(y)+2;
2148 : /* pay attention to sparse divisors */
2149 2075254 : for (i = 1; i <= dy; i++)
2150 1465067 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2151 610187 : dz = dx-dy; p = dz+1;
2152 610187 : av2 = avma;
2153 : for (;;)
2154 : {
2155 1308665 : gel(x,0) = gneg(gel(x,0)); p--;
2156 4440168 : for (i=1; i<=dy; i++)
2157 : {
2158 3131516 : GEN c = gmul(y_lead, gel(x,i));
2159 3131428 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2160 3131496 : gel(x,i) = rem(c, T);
2161 : }
2162 2276357 : for ( ; i<=dx; i++)
2163 : {
2164 967700 : GEN c = gmul(y_lead, gel(x,i));
2165 967692 : gel(x,i) = rem(c, T);
2166 : }
2167 1352555 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2168 1308657 : if (dx < dy) break;
2169 698477 : if (gc_needed(av2,1))
2170 : {
2171 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2172 0 : gc_slice(av2,x,dx+1);
2173 : }
2174 : }
2175 610180 : if (dx < 0) return pol_0(vx);
2176 610019 : lx = dx+3; x -= 2;
2177 610019 : x[0] = evaltyp(t_POL) | _evallg(lx);
2178 610019 : x[1] = evalsigne(1) | evalvarn(vx);
2179 610019 : x = RgX_recip_i(x);
2180 610020 : if (p)
2181 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2182 28809 : GEN t = y_lead;
2183 28809 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2184 0 : t = RgXQ_powu(t, p, T);
2185 : else
2186 28809 : t = gpowgs(t, p);
2187 88708 : for (i=2; i<lx; i++)
2188 : {
2189 59899 : GEN c = gmul(gel(x,i), t);
2190 59899 : gel(x,i) = rem(c,T);
2191 : }
2192 28809 : if (!T) return gc_upto(av, x);
2193 : }
2194 581211 : return gc_GEN(av, x);
2195 : }
2196 :
2197 : GEN
2198 660119 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2199 :
2200 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2201 : GEN
2202 12063 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2203 : {
2204 12063 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2205 12063 : pari_sp av = avma, av2;
2206 : GEN z, r, ypow, y_lead;
2207 :
2208 12063 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2209 12063 : dy = degpol(y); y_lead = gel(y,dy+2);
2210 12063 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2211 7561 : dx = degpol(x);
2212 7561 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2213 7561 : if (dx == dy)
2214 : {
2215 98 : GEN x_lead = gel(x,lg(x)-1);
2216 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2217 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2218 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2219 98 : *ptr = gc_upto(av, r); return scalarpol(x_lead, vx);
2220 : }
2221 7463 : (void)new_chunk(2);
2222 7463 : x = RgX_recip_i(x)+2;
2223 7463 : y = RgX_recip_i(y)+2;
2224 : /* pay attention to sparse divisors */
2225 39000 : for (i = 1; i <= dy; i++)
2226 31537 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2227 7463 : dz = dx-dy; p = dz+1;
2228 7463 : lz = dz+3;
2229 7463 : z = cgetg(lz, t_POL);
2230 7463 : z[1] = evalsigne(1) | evalvarn(vx);
2231 28640 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2232 7463 : ypow = new_chunk(dz+1);
2233 7463 : gel(ypow,0) = gen_1;
2234 7463 : gel(ypow,1) = y_lead;
2235 13714 : for (i=2; i<=dz; i++)
2236 : {
2237 6251 : GEN c = gmul(gel(ypow,i-1), y_lead);
2238 6251 : gel(ypow,i) = rem(c,T);
2239 : }
2240 7463 : av2 = avma;
2241 7463 : for (iz=2;;)
2242 : {
2243 15605 : p--;
2244 15605 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2245 69915 : for (i=1; i<=dy; i++)
2246 : {
2247 54310 : GEN c = gmul(y_lead, gel(x,i));
2248 54310 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2249 54310 : gel(x,i) = rem(c, T);
2250 : }
2251 41058 : for ( ; i<=dx; i++)
2252 : {
2253 25453 : GEN c = gmul(y_lead, gel(x,i));
2254 25453 : gel(x,i) = rem(c,T);
2255 : }
2256 15605 : x++; dx--;
2257 21177 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2258 15605 : if (dx < dy) break;
2259 8142 : if (gc_needed(av2,1))
2260 : {
2261 0 : GEN X = x-2;
2262 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2263 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2264 0 : (void)gc_all(av2,2, &X, &z); x = X+2;
2265 : }
2266 : }
2267 14008 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2268 7463 : if (dx < 0)
2269 182 : x = pol_0(vx);
2270 : else
2271 : {
2272 7281 : lx = dx+3; x -= 2;
2273 7281 : x[0] = evaltyp(t_POL) | _evallg(lx);
2274 7281 : x[1] = evalsigne(1) | evalvarn(vx);
2275 7281 : x = RgX_recip_i(x);
2276 : }
2277 7463 : z = RgX_recip_i(z);
2278 7463 : r = x;
2279 7463 : if (p)
2280 : {
2281 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2282 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2283 : }
2284 7463 : *ptr = r; return gc_all(av, 2, &z, ptr);
2285 : }
2286 : GEN
2287 11818 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2288 11818 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2289 :
2290 : GEN
2291 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2292 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2293 : GEN
2294 753034202 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2295 : GEN
2296 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2297 : GEN
2298 26320 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2299 : GEN
2300 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2301 : GEN
2302 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2303 :
2304 : GEN
2305 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2306 :
2307 : GEN
2308 0 : RgXQX_powers(GEN P, long n, GEN T)
2309 : {
2310 0 : GEN v = cgetg(n+2, t_VEC);
2311 : long i;
2312 0 : gel(v, 1) = pol_1(varn(T));
2313 0 : if (n==0) return v;
2314 0 : gel(v, 2) = gcopy(P);
2315 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2316 0 : return v;
2317 : }
2318 :
2319 : static GEN
2320 623387 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2321 : static GEN
2322 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2323 : static GEN
2324 67319 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2325 : static GEN
2326 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2327 : static GEN
2328 276652 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2329 : static GEN
2330 1040650 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2331 : static GEN
2332 1019726 : _one(void *data) { return pol_1(varn((GEN)data)); }
2333 : static GEN
2334 1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2335 : static GEN
2336 724301 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2337 :
2338 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2339 : _mul, _sqr, _one, _zero };
2340 :
2341 : GEN
2342 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2343 : {
2344 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2345 : }
2346 :
2347 : GEN
2348 417521 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2349 : {
2350 417521 : int use_sqr = 2*degpol(x) >= degpol(T);
2351 417519 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2352 : }
2353 :
2354 : /* mod X^n */
2355 : struct modXn {
2356 : long v; /* varn(X) */
2357 : long n;
2358 : } ;
2359 : static GEN
2360 11893 : _sqrXn(void *data, GEN x) {
2361 11893 : struct modXn *S = (struct modXn*)data;
2362 11893 : return RgXn_sqr(x, S->n);
2363 : }
2364 : static GEN
2365 4528 : _mulXn(void *data, GEN x, GEN y) {
2366 4528 : struct modXn *S = (struct modXn*)data;
2367 4528 : return RgXn_mul(x,y, S->n);
2368 : }
2369 : static GEN
2370 1939 : _oneXn(void *data) {
2371 1939 : struct modXn *S = (struct modXn*)data;
2372 1939 : return pol_1(S->v);
2373 : }
2374 : static GEN
2375 0 : _zeroXn(void *data) {
2376 0 : struct modXn *S = (struct modXn*)data;
2377 0 : return pol_0(S->v);
2378 : }
2379 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2380 : _oneXn, _zeroXn };
2381 :
2382 : GEN
2383 483 : RgXn_powers(GEN x, long m, long n)
2384 : {
2385 483 : long d = degpol(x);
2386 483 : int use_sqr = (d<<1) >= n;
2387 : struct modXn S;
2388 483 : S.v = varn(x); S.n = n;
2389 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2390 : }
2391 :
2392 : GEN
2393 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2394 : {
2395 : struct modXn S;
2396 : long v;
2397 2286 : if (n == 0) return x;
2398 2286 : v = RgX_valrem(x, &x);
2399 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2400 2265 : S.v = varn(x); S.n = n;
2401 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2402 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2403 2265 : return x;
2404 : }
2405 : GEN
2406 0 : RgXn_powu(GEN x, ulong m, long n)
2407 : {
2408 : pari_sp av;
2409 0 : if (n == 0) return gcopy(x);
2410 0 : av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
2411 : }
2412 :
2413 : GEN
2414 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2415 : {
2416 : struct modXn S;
2417 966 : S.v = varn(gel(x,2)); S.n = n;
2418 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2419 : }
2420 :
2421 : GEN
2422 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2423 : {
2424 0 : int use_sqr = 2*degpol(x) >= n;
2425 : struct modXn S;
2426 0 : S.v = varn(x); S.n = n;
2427 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2428 : }
2429 :
2430 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2431 : GEN
2432 5187 : RgXn_eval(GEN Q, GEN x, long n)
2433 : {
2434 5187 : long d = degpol(x);
2435 : int use_sqr;
2436 : struct modXn S;
2437 5187 : if (d == 1 && isrationalzero(gel(x,2)))
2438 : {
2439 5180 : GEN y = RgX_unscale(Q, gel(x,3));
2440 5180 : setvarn(y, varn(x)); return y;
2441 : }
2442 7 : S.v = varn(x);
2443 7 : S.n = n;
2444 7 : use_sqr = (d<<1) >= n;
2445 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2446 : }
2447 :
2448 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2449 : static GEN
2450 2614799 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2451 : {
2452 2614799 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2453 2614799 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2454 : }
2455 :
2456 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2457 : static GEN
2458 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2459 : {
2460 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2461 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2462 : }
2463 :
2464 : static GEN
2465 2188989 : RgXn_div_gen(GEN g, GEN f, long e)
2466 : {
2467 : pari_sp av;
2468 : ulong mask;
2469 : GEN W, a;
2470 2188989 : long v = varn(f), n = 1;
2471 :
2472 2188989 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2473 2188982 : a = ginv(gel(f,2));
2474 2188981 : if (e == 1 && !g) return scalarpol(a, v);
2475 2184067 : else if (e == 2 && !g)
2476 : {
2477 : GEN b;
2478 836936 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2479 246255 : b = gneg(b);
2480 246255 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2481 246256 : return deg1pol(b, a, v);
2482 : }
2483 1347131 : av = avma;
2484 1347131 : W = scalarpol_shallow(a,v);
2485 1347131 : mask = quadratic_prec_mask(e);
2486 3955973 : while (mask > 1)
2487 : {
2488 : GEN u, fr;
2489 2608842 : long n2 = n;
2490 2608842 : n<<=1; if (mask & 1) n--;
2491 2608842 : mask >>= 1;
2492 2608842 : fr = RgXn_red_shallow(f, n);
2493 2608842 : if (mask>1 || !g)
2494 : {
2495 1324855 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2496 1324855 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2497 : }
2498 : else
2499 : {
2500 1283987 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2501 1283987 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2502 1283987 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2503 : }
2504 2608842 : if (gc_needed(av,2))
2505 : {
2506 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2507 0 : W = gc_upto(av, W);
2508 : }
2509 : }
2510 1347131 : return W;
2511 : }
2512 :
2513 : static GEN
2514 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2515 : {
2516 : GEN r;
2517 119 : if (lgefint(p) == 3)
2518 : {
2519 119 : ulong pp = uel(p, 2);
2520 119 : if (pp == 2)
2521 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2522 : else
2523 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2524 : }
2525 : else
2526 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2527 119 : return FpX_to_mod(r, p);
2528 : }
2529 :
2530 : static GEN
2531 7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2532 : {
2533 7 : GEN r, T = RgX_to_FpX(pol, p);
2534 7 : if (signe(T) == 0) pari_err_OP("/", x, y);
2535 7 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2536 7 : return FpXQX_to_mod(r, T, p);
2537 : }
2538 :
2539 : static GEN
2540 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2541 : {
2542 : GEN r;
2543 91 : if (lgefint(p) == 3)
2544 : {
2545 91 : ulong pp = uel(p, 2);
2546 91 : if (pp == 2)
2547 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2548 : else
2549 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2550 : }
2551 : else
2552 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2553 91 : return FpX_to_mod(r, p);
2554 : }
2555 :
2556 : static GEN
2557 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2558 : {
2559 0 : GEN r, T = RgX_to_FpX(pol, p);
2560 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2561 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2562 0 : return FpXQX_to_mod(r, T, p);
2563 : }
2564 :
2565 : static GEN
2566 905091 : RgXn_inv_fast(GEN x, long e)
2567 : {
2568 : GEN p, pol;
2569 : long pa;
2570 905091 : long t = RgX_type(x,&p,&pol,&pa);
2571 905094 : switch(t)
2572 : {
2573 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2574 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2575 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2576 905003 : default: return NULL;
2577 : }
2578 : }
2579 :
2580 : static GEN
2581 1284113 : RgXn_div_fast(GEN x, GEN y, long e)
2582 : {
2583 : GEN p, pol;
2584 : long pa;
2585 1284113 : long t = RgX_type2(x,y,&p,&pol,&pa);
2586 1284113 : switch(t)
2587 : {
2588 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2589 7 : case RgX_type_code(t_POLMOD, t_INTMOD):
2590 7 : return RgXn_div_FpXQX(x, y, e, pol, p);
2591 1283987 : default: return NULL;
2592 : }
2593 : }
2594 :
2595 : GEN
2596 1284113 : RgXn_div_i(GEN g, GEN f, long e)
2597 : {
2598 1284113 : GEN h = RgXn_div_fast(g, f, e);
2599 1284113 : if (h) return h;
2600 1283987 : return RgXn_div_gen(g, f, e);
2601 : }
2602 :
2603 : GEN
2604 584450 : RgXn_div(GEN g, GEN f, long e)
2605 : {
2606 584450 : pari_sp av = avma;
2607 584450 : return gc_upto(av, RgXn_div_i(g, f, e));
2608 : }
2609 :
2610 : GEN
2611 905091 : RgXn_inv_i(GEN f, long e)
2612 : {
2613 905091 : GEN h = RgXn_inv_fast(f, e);
2614 905093 : if (h) return h;
2615 905002 : return RgXn_div_gen(NULL, f, e);
2616 : }
2617 :
2618 : GEN
2619 807480 : RgXn_inv(GEN f, long e)
2620 : {
2621 807480 : pari_sp av = avma;
2622 807480 : return gc_upto(av, RgXn_inv_i(f, e));
2623 : }
2624 :
2625 : /* intformal(x^n*S) / x^(n+1) */
2626 : static GEN
2627 56374 : RgX_integXn(GEN x, long n)
2628 114351 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2629 :
2630 : GEN
2631 52916 : RgXn_expint(GEN h, long e)
2632 : {
2633 52916 : pari_sp av = avma, av2;
2634 52916 : long v = varn(h), n;
2635 52916 : GEN f = pol_1(v), g;
2636 : ulong mask;
2637 :
2638 52916 : if (!signe(h)) return f;
2639 50424 : g = pol_1(v);
2640 50425 : n = 1; mask = quadratic_prec_mask(e);
2641 50426 : av2 = avma;
2642 56376 : for (;mask>1;)
2643 : {
2644 : GEN u, w;
2645 56376 : long n2 = n;
2646 56376 : n<<=1; if (mask & 1) n--;
2647 56376 : mask >>= 1;
2648 56376 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2649 56375 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2650 56375 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2651 56376 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2652 56374 : if (mask<=1) break;
2653 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2654 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2655 5950 : if (gc_needed(av2,2))
2656 : {
2657 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2658 0 : (void)gc_all(av2, 2, &f, &g);
2659 : }
2660 : }
2661 50424 : return gc_upto(av, f);
2662 : }
2663 :
2664 : GEN
2665 0 : RgXn_exp(GEN h, long e)
2666 : {
2667 0 : long d = degpol(h);
2668 0 : if (d < 0) return pol_1(varn(h));
2669 0 : if (!d || !gequal0(gel(h,2)))
2670 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2671 0 : return RgXn_expint(RgX_deriv(h), e);
2672 : }
2673 :
2674 : GEN
2675 154 : RgXn_reverse(GEN f, long e)
2676 : {
2677 154 : pari_sp av = avma, av2;
2678 : ulong mask;
2679 : GEN fi, a, df, W, an;
2680 154 : long v = varn(f), n=1;
2681 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2682 0 : pari_err_INV("serreverse",f);
2683 154 : fi = ginv(gel(f,3));
2684 154 : a = deg1pol_shallow(fi,gen_0,v);
2685 154 : if (e <= 2) return gc_GEN(av, a);
2686 133 : W = scalarpol(fi,v);
2687 133 : df = RgX_deriv(f);
2688 133 : mask = quadratic_prec_mask(e);
2689 133 : av2 = avma;
2690 616 : for (;mask>1;)
2691 : {
2692 : GEN u, fa, fr;
2693 483 : long n2 = n, rt;
2694 483 : n<<=1; if (mask & 1) n--;
2695 483 : mask >>= 1;
2696 483 : fr = RgXn_red_shallow(f, n);
2697 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2698 483 : an = RgXn_powers(a, rt, n);
2699 483 : if (n>1)
2700 : {
2701 483 : long n4 = (n2+1)>>1;
2702 483 : GEN dfr = RgXn_red_shallow(df, n2);
2703 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2704 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2705 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2706 : }
2707 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2708 483 : fa = RgX_shift(fa, -n2);
2709 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2710 483 : if (gc_needed(av2,2))
2711 : {
2712 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2713 0 : (void)gc_all(av2, 2, &a, &W);
2714 : }
2715 : }
2716 133 : return gc_upto(av, a);
2717 : }
2718 :
2719 : GEN
2720 7 : RgXn_sqrt(GEN h, long e)
2721 : {
2722 7 : pari_sp av = avma, av2;
2723 7 : long v = varn(h), n = 1;
2724 7 : GEN f = scalarpol(gen_1, v), df = f;
2725 7 : ulong mask = quadratic_prec_mask(e);
2726 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2727 0 : pari_err_SQRTN("RgXn_sqrt",h);
2728 7 : av2 = avma;
2729 : while(1)
2730 7 : {
2731 14 : long n2 = n, m;
2732 : GEN g;
2733 14 : n<<=1; if (mask & 1) n--;
2734 14 : mask >>= 1;
2735 14 : m = n-n2;
2736 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2737 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2738 14 : if (mask==1) return gc_upto(av, f);
2739 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2740 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2741 7 : if (gc_needed(av2,2))
2742 : {
2743 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2744 0 : (void)gc_all(av2, 2, &f, &df);
2745 : }
2746 : }
2747 : }
2748 :
2749 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2750 : GEN
2751 116731 : RgXQ_powu(GEN x, ulong n, GEN T)
2752 : {
2753 116731 : pari_sp av = avma;
2754 :
2755 116731 : if (!n) return pol_1(varn(x));
2756 68470 : if (n == 1) return RgX_copy(x);
2757 19824 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2758 19824 : return gc_GEN(av, x);
2759 : }
2760 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2761 : GEN
2762 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2763 : {
2764 : pari_sp av;
2765 102160 : long s = signe(n);
2766 :
2767 102160 : if (!s) return pol_1(varn(x));
2768 102160 : if (is_pm1(n) == 1)
2769 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2770 13853 : av = avma;
2771 13853 : if (s < 0) x = RgXQ_inv(x, T);
2772 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2773 13853 : return gc_GEN(av, x);
2774 : }
2775 : static GEN
2776 200602 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2777 : static GEN
2778 111194 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2779 :
2780 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2781 : GEN
2782 12145 : ZXQ_powers(GEN x, long l, GEN T)
2783 : {
2784 12145 : int use_sqr = 2*degpol(x) >= degpol(T);
2785 12145 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2786 : }
2787 :
2788 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2789 : GEN
2790 168507 : ZXQ_powu(GEN x, ulong n, GEN T)
2791 : {
2792 168507 : pari_sp av = avma;
2793 :
2794 168507 : if (!n) return pol_1(varn(x));
2795 168507 : if (n == 1) return ZX_copy(x);
2796 113873 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2797 113874 : return gc_GEN(av, x);
2798 : }
2799 :
2800 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2801 : GEN
2802 8596 : RgXQ_powers(GEN x, long l, GEN T)
2803 : {
2804 8596 : int use_sqr = 2*degpol(x) >= degpol(T);
2805 8596 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2806 : }
2807 :
2808 : GEN
2809 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2810 : {
2811 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2812 : }
2813 :
2814 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2815 : GEN
2816 9590 : QXQ_powers(GEN a, long n, GEN T)
2817 : {
2818 : GEN den, v;
2819 9590 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2820 9562 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2821 : /* den*a integral; v[i+1] = (den*a)^i in K */
2822 9562 : if (den)
2823 : { /* restore denominators */
2824 6012 : GEN d = den;
2825 : long i;
2826 6012 : gel(v,2) = a;
2827 27749 : for (i=3; i<=n+1; i++) {
2828 21737 : d = mulii(d,den);
2829 21737 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2830 : }
2831 : }
2832 9562 : return v;
2833 : }
2834 :
2835 : static GEN
2836 3612 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2837 : {
2838 3612 : long l, i, m = 0;
2839 : GEN dz, z;
2840 3612 : GEN V = cgetg_copy(v, &l);
2841 12796 : for (i = imin; i < l; i++)
2842 : {
2843 9184 : GEN c = gel(v, i);
2844 9184 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2845 : }
2846 3612 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2847 3983 : for (i = 1; i < imin; i++) V[i] = v[i];
2848 12796 : for (i = imin; i < l; i++)
2849 : {
2850 9184 : GEN c = gel(v,i);
2851 9184 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2852 9184 : gel(V,i) = c;
2853 : }
2854 3612 : return V;
2855 : }
2856 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2857 : GEN
2858 3241 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2859 3241 : { return do_QXQ_eval(v, 1, a, T); }
2860 :
2861 : GEN
2862 371 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2863 371 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2864 :
2865 : GEN
2866 2940 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2867 : {
2868 2940 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2869 : }
2870 :
2871 : GEN
2872 2680 : RgXQ_norm(GEN x, GEN T)
2873 : {
2874 : pari_sp av;
2875 2680 : long dx = degpol(x);
2876 : GEN L, y;
2877 2680 : if (degpol(T)==0) return gpowgs(x,0);
2878 2673 : av = avma; y = resultant(T, x);
2879 2673 : L = leading_coeff(T);
2880 2673 : if (gequal1(L) || !signe(x)) return y;
2881 0 : return gc_upto(av, gdiv(y, gpowgs(L, dx)));
2882 : }
2883 :
2884 : GEN
2885 476 : RgXQ_trace(GEN x, GEN T)
2886 : {
2887 476 : pari_sp av = avma;
2888 : GEN dT, z;
2889 : long n;
2890 476 : if (degpol(T)==0) return gmulgs(x,0);
2891 469 : dT = RgX_deriv(T); n = degpol(dT);
2892 469 : z = RgXQ_mul(x, dT, T);
2893 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2894 420 : return gc_upto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2895 : }
2896 :
2897 : GEN
2898 2788484 : RgX_blocks(GEN P, long n, long m)
2899 : {
2900 2788484 : GEN z = cgetg(m+1,t_VEC);
2901 2788484 : long i,j, k=2, l = lg(P);
2902 8559030 : for(i=1; i<=m; i++)
2903 : {
2904 5770546 : GEN zi = cgetg(n+2,t_POL);
2905 5770546 : zi[1] = P[1];
2906 5770546 : gel(z,i) = zi;
2907 17024680 : for(j=2; j<n+2; j++)
2908 11254134 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2909 5770546 : zi = RgX_renormalize_lg(zi, n+2);
2910 : }
2911 2788484 : return z;
2912 : }
2913 :
2914 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2915 : void
2916 49757556 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2917 : {
2918 49757556 : long n = degpol(p), v = varn(p), n0, n1, i;
2919 : GEN p0, p1;
2920 :
2921 49759157 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2922 :
2923 49665873 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2924 49665873 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2925 49668663 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2926 146784882 : for (i=0; i<n1; i++)
2927 : {
2928 97116129 : p0[2+i] = p[2+(i<<1)];
2929 97116129 : p1[2+i] = p[3+(i<<1)];
2930 : }
2931 49668753 : if (n1 != n0)
2932 18024996 : p0[2+i] = p[2+(i<<1)];
2933 49668753 : *pe = normalizepol(p0);
2934 49673593 : *po = normalizepol(p1);
2935 : }
2936 :
2937 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2938 : GEN
2939 43631 : RgX_splitting(GEN p, long k)
2940 : {
2941 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
2942 : GEN r;
2943 :
2944 43631 : m = n/k;
2945 43631 : r = cgetg(k+1,t_VEC);
2946 234059 : for(i=1; i<=k; i++)
2947 : {
2948 190428 : gel(r,i) = cgetg(m+3, t_POL);
2949 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2950 : }
2951 571991 : for (j=1, i=0, l=2; i<=n; i++)
2952 : {
2953 528360 : gmael(r,j,l) = gel(p,2+i);
2954 528360 : if (j==k) { j=1; l++; } else j++;
2955 : }
2956 234059 : for(i=1; i<=k; i++)
2957 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2958 43631 : return r;
2959 : }
2960 :
2961 : /*******************************************************************/
2962 : /* */
2963 : /* Kronecker form */
2964 : /* */
2965 : /*******************************************************************/
2966 :
2967 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2968 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2969 : * normalized coefficients */
2970 : GEN
2971 186210 : Kronecker_to_mod(GEN z, GEN T)
2972 : {
2973 186210 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2974 186210 : GEN x, t = cgetg(N,t_POL);
2975 186210 : t[1] = T[1];
2976 186210 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2977 186210 : x[1] = z[1];
2978 186210 : T = RgX_copy(T);
2979 1399501 : for (i=2; i<lx+2; i++, z+= N-2)
2980 : {
2981 5761438 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2982 1213291 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2983 : }
2984 186210 : N = (l-2) % (N-2) + 2;
2985 472263 : for (j=2; j<N; j++) t[j] = z[j];
2986 186210 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2987 186210 : return normalizepol_lg(x, i+1);
2988 : }
2989 :
2990 : /*******************************************************************/
2991 : /* */
2992 : /* Domain detection */
2993 : /* */
2994 : /*******************************************************************/
2995 :
2996 : static GEN
2997 609812 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2998 : {
2999 609812 : pari_sp av = avma;
3000 : GEN r;
3001 609812 : if (lgefint(p) == 3)
3002 : {
3003 558456 : ulong pp = uel(p, 2);
3004 558456 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
3005 : RgX_to_Flx(y, pp), pp));
3006 : }
3007 : else
3008 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3009 609812 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3010 544614 : return gc_upto(av, FpX_to_mod(r, p));
3011 : }
3012 :
3013 : static GEN
3014 2079 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3015 : {
3016 2079 : pari_sp av = avma;
3017 : long dT;
3018 : GEN kx, ky, r;
3019 2079 : GEN T = RgX_to_FpX(pol, p);
3020 2079 : if (signe(T)==0) pari_err_OP("*", x, y);
3021 2072 : dT = degpol(T);
3022 2072 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3023 2072 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3024 2072 : r = FpX_mul(kx, ky, p);
3025 2072 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3026 1463 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3027 : }
3028 :
3029 : static GEN
3030 455393 : RgX_liftred(GEN x, GEN T)
3031 455393 : { return RgXQX_red(liftpol_shallow(x), T); }
3032 :
3033 : static GEN
3034 171133 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3035 : {
3036 171133 : pari_sp av = avma;
3037 171133 : long dT = degpol(T);
3038 171133 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3039 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3040 171133 : return gc_upto(av, Kronecker_to_mod(r, T));
3041 : }
3042 :
3043 : static GEN
3044 1407 : RgX_sqr_FpX(GEN x, GEN p)
3045 : {
3046 1407 : pari_sp av = avma;
3047 : GEN r;
3048 1407 : if (lgefint(p) == 3)
3049 : {
3050 1203 : ulong pp = uel(p, 2);
3051 1203 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3052 : }
3053 : else
3054 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3055 1407 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3056 1274 : return gc_upto(av, FpX_to_mod(r, p));
3057 : }
3058 :
3059 : static GEN
3060 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3061 : {
3062 196 : pari_sp av = avma;
3063 : long dT;
3064 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3065 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3066 189 : dT = degpol(T);
3067 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3068 189 : r = FpX_sqr(kx, p);
3069 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3070 189 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3071 : }
3072 :
3073 : static GEN
3074 13425 : RgX_sqr_QXQX(GEN x, GEN T)
3075 : {
3076 13425 : pari_sp av = avma;
3077 13425 : long dT = degpol(T);
3078 13425 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3079 13425 : return gc_upto(av, Kronecker_to_mod(r, T));
3080 : }
3081 :
3082 : static GEN
3083 68327 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3084 : {
3085 68327 : pari_sp av = avma;
3086 : GEN r;
3087 68327 : if (lgefint(p) == 3)
3088 : {
3089 51691 : ulong pp = uel(p, 2);
3090 51691 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3091 : RgX_to_Flx(y, pp), pp));
3092 : }
3093 : else
3094 16636 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3095 68327 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3096 31031 : return gc_upto(av, FpX_to_mod(r, p));
3097 : }
3098 :
3099 : static GEN
3100 49851 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3101 : {
3102 49851 : pari_sp av = avma;
3103 : GEN r;
3104 49851 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3105 49851 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3106 : }
3107 : static GEN
3108 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3109 : {
3110 70 : pari_sp av = avma;
3111 : GEN r;
3112 70 : GEN T = RgX_to_FpX(pol, p);
3113 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3114 63 : if (lgefint(p) == 3)
3115 : {
3116 55 : ulong pp = uel(p, 2);
3117 55 : GEN Tp = ZX_to_Flx(T, pp);
3118 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3119 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3120 : }
3121 : else
3122 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3123 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3124 56 : return gc_upto(av, FpXQX_to_mod(r, T, p));
3125 : }
3126 :
3127 : static GEN
3128 121568271 : RgX_mul_fast(GEN x, GEN y)
3129 : {
3130 : GEN p, pol;
3131 : long pa;
3132 121568271 : long t = RgX_type2(x,y, &p,&pol,&pa);
3133 121568553 : switch(t)
3134 : {
3135 67620441 : case t_INT: return ZX_mul(x,y);
3136 3357958 : case t_FRAC: return QX_mul(x,y);
3137 105127 : case t_FFELT: return FFX_mul(x, y, pol);
3138 609812 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3139 171133 : case RgX_type_code(t_POLMOD, t_INT):
3140 : case RgX_type_code(t_POLMOD, t_FRAC):
3141 171133 : return RgX_mul_QXQX(x, y, pol);
3142 2042 : case RgX_type_code(t_POLMOD, t_INTMOD):
3143 2042 : return RgX_mul_FpXQX(x, y, pol, p);
3144 49702040 : default: return NULL;
3145 : }
3146 : }
3147 : static GEN
3148 1385650 : RgX_sqr_fast(GEN x)
3149 : {
3150 : GEN p, pol;
3151 : long pa;
3152 1385650 : long t = RgX_type(x,&p,&pol,&pa);
3153 1385665 : switch(t)
3154 : {
3155 1246912 : case t_INT: return ZX_sqr(x);
3156 83269 : case t_FRAC: return QX_sqr(x);
3157 2499 : case t_FFELT: return FFX_sqr(x, pol);
3158 1407 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3159 13425 : case RgX_type_code(t_POLMOD, t_INT):
3160 : case RgX_type_code(t_POLMOD, t_FRAC):
3161 13425 : return RgX_sqr_QXQX(x, pol);
3162 192 : case RgX_type_code(t_POLMOD, t_INTMOD):
3163 192 : return RgX_sqr_FpXQX(x, pol, p);
3164 37961 : default: return NULL;
3165 : }
3166 : }
3167 :
3168 : static GEN
3169 14716746 : RgX_rem_fast(GEN x, GEN y)
3170 : {
3171 : GEN p, pol;
3172 : long pa;
3173 14716746 : long t = RgX_type2(x,y, &p,&pol,&pa);
3174 14716843 : switch(t)
3175 : {
3176 3492001 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3177 1837192 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3178 84 : case t_FFELT: return FFX_rem(x, y, pol);
3179 68327 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3180 49851 : case RgX_type_code(t_POLMOD, t_INT):
3181 : case RgX_type_code(t_POLMOD, t_FRAC):
3182 49851 : return RgX_rem_QXQX(x, y, pol);
3183 66 : case RgX_type_code(t_POLMOD, t_INTMOD):
3184 66 : return RgX_rem_FpXQX(x, y, pol, p);
3185 9269322 : default: return NULL;
3186 : }
3187 : }
3188 :
3189 : GEN
3190 107835097 : RgX_mul(GEN x, GEN y)
3191 : {
3192 107835097 : GEN z = RgX_mul_fast(x,y);
3193 107835198 : if (!z) z = RgX_mul_i(x,y);
3194 107834838 : return z;
3195 : }
3196 :
3197 : GEN
3198 1373498 : RgX_sqr(GEN x)
3199 : {
3200 1373498 : GEN z = RgX_sqr_fast(x);
3201 1373509 : if (!z) z = RgX_sqr_i(x);
3202 1373509 : return z;
3203 : }
3204 :
3205 : GEN
3206 14716764 : RgX_rem(GEN x, GEN y)
3207 : {
3208 14716764 : GEN z = RgX_rem_fast(x, y);
3209 14716831 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3210 14716801 : return z;
3211 : }
3212 :
3213 : static GEN
3214 0 : _RgX_mul(void* E, GEN x, GEN y)
3215 0 : { (void) E; return RgX_mul(x, y); }
3216 :
3217 : GEN
3218 0 : RgXV_prod(GEN V)
3219 0 : { return gen_product(V, NULL, &_RgX_mul); }
3220 :
3221 : GEN
3222 11061993 : RgXn_mul(GEN f, GEN g, long n)
3223 : {
3224 11061993 : pari_sp av = avma;
3225 11061993 : GEN h = RgX_mul_fast(f,g);
3226 11061995 : if (!h) return RgXn_mul2(f,g,n);
3227 1783591 : if (degpol(h) < n) return h;
3228 383046 : return gc_GEN(av, RgXn_red_shallow(h, n));
3229 : }
3230 :
3231 : GEN
3232 12152 : RgXn_sqr(GEN f, long n)
3233 : {
3234 12152 : pari_sp av = avma;
3235 12152 : GEN g = RgX_sqr_fast(f);
3236 12152 : if (!g) return RgXn_sqr2(f,n);
3237 11837 : if (degpol(g) < n) return g;
3238 10640 : return gc_GEN(av, RgXn_red_shallow(g, n));
3239 : }
3240 :
3241 : /* (f*g) \/ x^n */
3242 : GEN
3243 2671189 : RgX_mulhigh_i(GEN f, GEN g, long n)
3244 : {
3245 2671189 : GEN h = RgX_mul_fast(f,g);
3246 2671189 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3247 : }
3248 :
3249 : /* (f*g) \/ x^n */
3250 : GEN
3251 0 : RgX_sqrhigh_i(GEN f, long n)
3252 : {
3253 0 : GEN h = RgX_sqr_fast(f);
3254 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3255 : }
|
