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 19437467 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19437467 : long pold=1, rold=n*(d-1);
35 91143886 : for(p=2; p<=d; p++)
36 : {
37 71706419 : r = m*(p-1) + n*((d-1)/p);
38 71706419 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19437467 : return pold;
41 : }
42 :
43 : static GEN
44 14560542 : 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 14560542 : pari_sp av = avma;
48 : long i;
49 14560542 : GEN z = cmul(E,P,a,ff->one(E));
50 14537960 : if (!z) z = gen_0;
51 76396967 : for (i=1; i<=n; i++)
52 : {
53 61860616 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 61879897 : if (t) {
55 46946895 : z = ff->add(E, z, t);
56 46919479 : if (gc_needed(av,2)) z = gerepileupto(av, z);
57 : }
58 : }
59 14536351 : 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 12467139 : 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 12467139 : pari_sp av = avma;
72 12467139 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12467139 : if (d < 0) return ff->zero(E);
76 11859385 : if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1157461 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1157461 : 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 1157461 : d -= l;
84 1157461 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2700277 : while (d >= l-1)
86 : {
87 1541507 : d -= l-1;
88 1541507 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1541418 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1541482 : if (gc_needed(av,2))
91 91 : z = gerepileupto(av, z);
92 : }
93 1158770 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1158773 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1158785 : return gerepileupto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 873956 : 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 873956 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 873956 : if (d < 0) return ff->zero(E);
106 872661 : rtd = (long) sqrt((double)d);
107 872661 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 872664 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 872657 : return gerepileupto(av, z);
110 : }
111 :
112 : static GEN
113 2027653 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19241858 : _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 1839783 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 591680 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2068754 : _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 668729 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 668729 : pari_sp av = avma, btop;
163 668729 : long i, d = degpol(P), o;
164 : GEN s;
165 668726 : if (signe(P)==0) return pol_0(varn(P));
166 668725 : s = gel(P, d+2);
167 668725 : if (d == 0) return gcopy(s);
168 665449 : o = RgX_deflate_order(P);
169 665471 : if (o > 1) A = gpowgs(A, o);
170 665483 : btop = avma;
171 2313693 : for (i = d-o; i >= 0; i-=o)
172 : {
173 1648214 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 1648210 : if (gc_needed(btop,1))
175 : {
176 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
177 13 : s = gerepileupto(btop, s);
178 : }
179 : }
180 665479 : return gerepileupto(av, s);
181 : }
182 :
183 : GEN
184 1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
185 : {
186 1652 : pari_sp av = avma;
187 1652 : long i, d = degpol(P), v = varn(A);
188 : GEN s;
189 1652 : if (signe(P)==0) return pol_0(v);
190 1652 : if (d == 0) return scalarpol(gel(P, d+2), v);
191 1232 : s = scalarpol_shallow(gel(P, d+2), v);
192 4963 : for (i = d-1; i >= 0; i--)
193 : {
194 3731 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
195 3731 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
196 3731 : if (gc_needed(av,1))
197 : {
198 0 : if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
199 0 : s = gerepileupto(av, s);
200 : }
201 : }
202 1232 : return gerepileupto(av, s);
203 : }
204 :
205 : const struct bb_algebra *
206 281051 : get_Rg_algebra(void)
207 : {
208 281051 : 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 24838187 : RgX_equal(GEN x, GEN y)
236 : {
237 24838187 : long i = lg(x);
238 :
239 24838187 : if (i != lg(y)) return 0;
240 106915584 : for (i--; i > 1; i--)
241 82405327 : if (!gequal(gel(x,i),gel(y,i))) return 0;
242 24510257 : 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 156946164 : Rg_get_1(GEN x)
249 : {
250 : GEN p, T;
251 156946164 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
252 156946164 : if (RgX_type_is_composite(tx))
253 11657560 : RgX_type_decode(tx, &i /*junk*/, &tx);
254 156946164 : 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 156511863 : 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 6942912 : Rg_get_0(GEN x)
266 : {
267 : GEN p, T;
268 6942912 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
269 6942912 : if (RgX_type_is_composite(tx))
270 61334 : RgX_type_decode(tx, &i /*junk*/, &tx);
271 6942912 : 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 6942163 : default: return gen_0;
277 : }
278 : }
279 :
280 : GEN
281 6874 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
282 : {
283 6874 : long i, n = degpol(P);
284 : GEN z, dz, dP;
285 6874 : if (n < 0) return gen_0;
286 6874 : P = Q_remove_denom(P, &dP);
287 6874 : z = gel(P,2); if (n == 0) return icopy(z);
288 3878 : if (dV) z = mulii(dV, z); /* V[1] = dV */
289 3878 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
290 7483 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
291 3878 : dz = mul_denom(dP, dV);
292 3878 : return dz? RgX_Rg_div(z, dz): z;
293 : }
294 :
295 : /* Return P(h * x), not memory clean */
296 : GEN
297 51360 : RgX_unscale(GEN P, GEN h)
298 : {
299 51360 : long i, l = lg(P);
300 51360 : GEN hi = gen_1, Q = cgetg(l, t_POL);
301 51360 : Q[1] = P[1];
302 51360 : if (l == 2) return Q;
303 40118 : gel(Q,2) = gcopy(gel(P,2));
304 94473 : for (i=3; i<l; i++)
305 : {
306 54355 : hi = gmul(hi,h);
307 54355 : gel(Q,i) = gmul(gel(P,i), hi);
308 : }
309 40118 : return Q;
310 : }
311 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
312 : GEN
313 1467268 : ZX_z_unscale(GEN P, long h)
314 : {
315 1467268 : long i, l = lg(P);
316 1467268 : GEN Q = cgetg(l, t_POL);
317 1467268 : Q[1] = P[1];
318 1467268 : if (l == 2) return Q;
319 1378722 : gel(Q,2) = gel(P,2);
320 1378722 : if (l == 3) return Q;
321 1352416 : if (h == -1)
322 236959 : for (i = 3; i < l; i++)
323 : {
324 196378 : gel(Q,i) = negi(gel(P,i));
325 196379 : if (++i == l) break;
326 145375 : gel(Q,i) = gel(P,i);
327 : }
328 : else
329 : {
330 : GEN hi;
331 1260832 : gel(Q,3) = mulis(gel(P,3), h);
332 1260832 : hi = sqrs(h);
333 5415797 : for (i = 4; i < l; i++)
334 : {
335 4154965 : gel(Q,i) = mulii(gel(P,i), hi);
336 4154964 : if (i != l-1) hi = mulis(hi,h);
337 : }
338 : }
339 1352417 : 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 924397 : ZX_unscale2n(GEN P, long n)
365 : {
366 924397 : long i, ni = n, l = lg(P);
367 924397 : GEN Q = cgetg(l, t_POL);
368 924396 : Q[1] = P[1];
369 924396 : if (l == 2) return Q;
370 924396 : gel(Q,2) = gel(P,2);
371 924396 : if (l == 3) return Q;
372 924396 : gel(Q,3) = shifti(gel(P,3), ni);
373 3166750 : for (i=4; i<l; i++)
374 : {
375 2242388 : ni += n;
376 2242388 : gel(Q,i) = shifti(gel(P,i), ni);
377 : }
378 924362 : 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 6223 : RgXV_unscale(GEN x, GEN h)
419 : {
420 6223 : if (isint1(h)) return gcopy(x);
421 18176 : 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 4342208 : RgX_rescale(GEN P, GEN h)
427 : {
428 4342208 : long i, l = lg(P);
429 4342208 : GEN Q = cgetg(l,t_POL), hi = h;
430 4342201 : gel(Q,l-1) = gel(P,l-1);
431 11427138 : for (i=l-2; i>=2; i--)
432 : {
433 11424529 : gel(Q,i) = gmul(gel(P,i), hi);
434 11424379 : if (i == 2) break;
435 7084611 : hi = gmul(hi,h);
436 : }
437 4342377 : 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 1175173 : RgX_deflate(GEN x0, long d)
450 : {
451 : GEN z, y, x;
452 1175173 : long i,id, dy, dx = degpol(x0);
453 1175170 : if (d == 1 || dx <= 0) return leafcopy(x0);
454 459089 : dy = dx/d;
455 459089 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
456 459087 : z = y + 2;
457 459087 : x = x0+ 2;
458 1762689 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
459 459087 : 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 126 : rfrac_deflate_order(GEN F)
485 : {
486 126 : GEN N = gel(F,1), D = gel(F,2);
487 126 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
488 126 : if (m == 1) return 1;
489 42 : if (typ(N) == t_POL && varn(N) == varn(D))
490 28 : m = cgcd(m, RgX_deflate_order(N));
491 42 : return m;
492 : }
493 : /* F a t_RFRAC */
494 : GEN
495 126 : rfrac_deflate_max(GEN F, long *m)
496 : {
497 126 : *m = rfrac_deflate_order(F);
498 126 : return rfrac_deflate(F, *m);
499 : }
500 : /* F a t_RFRAC */
501 : GEN
502 126 : rfrac_deflate(GEN F, long m)
503 : {
504 126 : GEN N = gel(F,1), D = gel(F,2);
505 126 : if (m == 1) return F;
506 42 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
507 42 : D = RgX_deflate(D, m); return mkrfrac(N, D);
508 : }
509 :
510 : /* return x0(X^d) */
511 : GEN
512 881758 : RgX_inflate(GEN x0, long d)
513 : {
514 881758 : long i, id, dy, dx = degpol(x0);
515 881758 : GEN x = x0 + 2, z, y;
516 881758 : if (dx <= 0) return leafcopy(x0);
517 804468 : dy = dx*d;
518 804468 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
519 804468 : z = y + 2;
520 28075875 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
521 11486543 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
522 804468 : return y;
523 : }
524 :
525 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
526 : static GEN
527 5701655 : RgX_translate_basecase(GEN P, GEN c)
528 : {
529 5701655 : pari_sp av = avma;
530 : GEN Q, R;
531 : long i, k, n;
532 :
533 5701655 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
534 5699725 : Q = leafcopy(P);
535 5699725 : R = Q+2; n = degpol(P);
536 5699725 : if (isint1(c))
537 : {
538 7595 : for (i=1; i<=n; i++)
539 : {
540 20265 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
541 5390 : if (gc_needed(av,2))
542 : {
543 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
544 0 : Q = gc_GEN(av, Q); R = Q+2;
545 : }
546 : }
547 : }
548 5697517 : else if (isintm1(c))
549 : {
550 15974 : for (i=1; i<=n; i++)
551 : {
552 49630 : for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
553 12299 : if (gc_needed(av,2))
554 : {
555 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
556 0 : Q = gc_GEN(av, Q); R = Q+2;
557 : }
558 : }
559 : }
560 : else
561 : {
562 19424578 : for (i=1; i<=n; i++)
563 : {
564 45449872 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
565 13730738 : if (gc_needed(av,2))
566 : {
567 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
568 0 : Q = gc_GEN(av, Q); R = Q+2;
569 : }
570 : }
571 : }
572 5699389 : return gc_GEN(av, Q);
573 : }
574 :
575 : static GEN
576 102627 : zero_FpX_mod(GEN p, long v)
577 : {
578 102627 : GEN r = cgetg(3,t_POL);
579 102627 : r[1] = evalvarn(v);
580 102627 : gel(r,2) = mkintmod(gen_0, icopy(p));
581 102627 : return r;
582 : }
583 :
584 : static GEN
585 0 : RgX_translate_FpX(GEN P, GEN c, GEN p)
586 : {
587 0 : pari_sp av = avma;
588 : GEN r;
589 : #if 0
590 : if (lgefint(p) == 3)
591 : {
592 : ulong pp = uel(p, 2);
593 : r = Flx_to_ZX_inplace(Flx_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
594 : }
595 : else
596 : #endif
597 0 : r = FpX_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
598 0 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
599 0 : return gerepileupto(av, FpX_to_mod(r, p));
600 : }
601 :
602 : static GEN
603 6063357 : RgX_translate_fast(GEN P, GEN c)
604 : {
605 : GEN p, pol;
606 : long pa;
607 6063357 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
608 6063343 : switch(t)
609 : {
610 361875 : case t_INT: return ZX_translate(P, c);
611 0 : case t_INTMOD: return RgX_translate_FpX(P, c, p);
612 5701468 : default: return NULL;
613 : }
614 : }
615 :
616 : static GEN
617 5701843 : RgX_translate_i(GEN P, GEN c)
618 : {
619 5701843 : pari_sp av = avma;
620 : long n;
621 5701843 : n = degpol(P);
622 5701844 : if (n < 40)
623 5701655 : return RgX_translate_basecase(P, c);
624 : else
625 : {
626 189 : long d = n >> 1;
627 189 : GEN Q = RgX_translate_i(RgX_shift_shallow(P, -d), c);
628 189 : GEN R = RgX_translate_i(RgXn_red_shallow(P, d), c);
629 189 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
630 189 : return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
631 : }
632 : }
633 :
634 : GEN
635 6063362 : RgX_translate(GEN P, GEN c)
636 : {
637 6063362 : GEN R = RgX_translate_fast(P, c);
638 6063344 : return R ? R: RgX_translate_i(P,c);
639 : }
640 : /* P(ax + b) */
641 : GEN
642 30212 : RgX_affine(GEN P, GEN a, GEN b)
643 : {
644 30212 : if (!gequal0(b)) P = RgX_translate(P, b);
645 30212 : return RgX_unscale(P, a);
646 : }
647 :
648 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
649 : GEN
650 32464 : RgXQX_translate(GEN P, GEN c, GEN T)
651 : {
652 32464 : pari_sp av = avma;
653 : GEN Q, R;
654 : long i, k, n;
655 :
656 32464 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
657 32121 : Q = leafcopy(P);
658 32121 : R = Q+2; n = degpol(P);
659 103762 : for (i=1; i<=n; i++)
660 : {
661 301670 : for (k=n-i; k<n; k++)
662 : {
663 230029 : pari_sp av2 = avma;
664 230029 : gel(R,k) = gerepileupto(av2,
665 230029 : RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
666 : }
667 71641 : if (gc_needed(av,2))
668 : {
669 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
670 0 : Q = gc_GEN(av, Q); R = Q+2;
671 : }
672 : }
673 32121 : return gc_GEN(av, Q);
674 : }
675 :
676 : /********************************************************************/
677 : /** **/
678 : /** CONVERSIONS **/
679 : /** (not memory clean) **/
680 : /** **/
681 : /********************************************************************/
682 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
683 : * but everything else is */
684 : static GEN
685 166663 : QXQ_to_mod(GEN x, GEN T)
686 : {
687 : long d;
688 166663 : switch(typ(x))
689 : {
690 62395 : case t_INT: return icopy(x);
691 2079 : case t_FRAC: return gcopy(x);
692 102189 : case t_POL:
693 102189 : d = degpol(x);
694 102189 : if (d < 0) return gen_0;
695 99928 : if (d == 0) return gcopy(gel(x,2));
696 97817 : return mkpolmod(RgX_copy(x), T);
697 0 : default: pari_err_TYPE("QXQ_to_mod",x);
698 : return NULL;/* LCOV_EXCL_LINE */
699 : }
700 : }
701 : /* pure shallow version */
702 : GEN
703 839309 : QXQ_to_mod_shallow(GEN x, GEN T)
704 : {
705 : long d;
706 839309 : switch(typ(x))
707 : {
708 544308 : case t_INT:
709 544308 : case t_FRAC: return x;
710 295001 : case t_POL:
711 295001 : d = degpol(x);
712 295001 : if (d < 0) return gen_0;
713 248750 : if (d == 0) return gel(x,2);
714 232150 : return mkpolmod(x, T);
715 0 : default: pari_err_TYPE("QXQ_to_mod",x);
716 : return NULL;/* LCOV_EXCL_LINE */
717 : }
718 : }
719 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
720 : * Not memory clean because T not copied, but everything else is */
721 : static GEN
722 37184 : QXQX_to_mod(GEN z, GEN T)
723 : {
724 37184 : long i,l = lg(z);
725 37184 : GEN x = cgetg(l,t_POL);
726 185598 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
727 37184 : x[1] = z[1]; return normalizepol_lg(x,l);
728 : }
729 : /* pure shallow version */
730 : GEN
731 203487 : QXQX_to_mod_shallow(GEN z, GEN T)
732 : {
733 203487 : long i,l = lg(z);
734 203487 : GEN x = cgetg(l,t_POL);
735 959608 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
736 203487 : x[1] = z[1]; return normalizepol_lg(x,l);
737 : }
738 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
739 : GEN
740 12516 : QXQXV_to_mod(GEN V, GEN T)
741 : {
742 12516 : long i, l = lg(V);
743 12516 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
744 49700 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
745 12516 : return z;
746 : }
747 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
748 : GEN
749 18358 : QXQV_to_mod(GEN V, GEN T)
750 : {
751 18358 : long i, l = lg(V);
752 18358 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
753 36607 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
754 18358 : return z;
755 : }
756 :
757 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
758 : GEN
759 14644 : QXQC_to_mod_shallow(GEN V, GEN T)
760 : {
761 14644 : long i, l = lg(V);
762 14644 : GEN z = cgetg(l, t_COL);
763 97832 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
764 14644 : return z;
765 : }
766 :
767 : GEN
768 6622 : QXQM_to_mod_shallow(GEN V, GEN T)
769 : {
770 6622 : long i, l = lg(V);
771 6622 : GEN z = cgetg(l, t_MAT);
772 21266 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
773 6622 : return z;
774 : }
775 :
776 : GEN
777 7885673 : RgX_renormalize_lg(GEN x, long lx)
778 : {
779 : long i;
780 10927274 : for (i = lx-1; i>1; i--)
781 10474334 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
782 7885673 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
783 7885673 : setlg(x, i+1); setsigne(x, i != 1); return x;
784 : }
785 :
786 : GEN
787 1492276 : RgV_to_RgX(GEN x, long v)
788 : {
789 1492276 : long i, k = lg(x);
790 : GEN p;
791 :
792 4152522 : while (--k && gequal0(gel(x,k)));
793 1492276 : if (!k) return pol_0(v);
794 1485026 : i = k+2; p = cgetg(i,t_POL);
795 1485025 : p[1] = evalsigne(1) | evalvarn(v);
796 10228782 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
797 1485025 : return p;
798 : }
799 : GEN
800 188283 : RgV_to_RgX_reverse(GEN x, long v)
801 : {
802 188283 : long j, k, l = lg(x);
803 : GEN p;
804 :
805 189851 : for (k = 1; k < l; k++)
806 189851 : if (!gequal0(gel(x,k))) break;
807 188283 : if (k == l) return pol_0(v);
808 188283 : k -= 1;
809 188283 : l -= k;
810 188283 : x += k;
811 188283 : p = cgetg(l+1,t_POL);
812 188283 : p[1] = evalsigne(1) | evalvarn(v);
813 985038 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
814 188283 : return p;
815 : }
816 :
817 : /* return the (N-dimensional) vector of coeffs of p */
818 : GEN
819 13407670 : RgX_to_RgC(GEN x, long N)
820 : {
821 : long i, l;
822 : GEN z;
823 13407670 : l = lg(x)-1; x++;
824 13407670 : if (l > N+1) l = N+1; /* truncate higher degree terms */
825 13407670 : z = cgetg(N+1,t_COL);
826 80850896 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
827 24469386 : for ( ; i<=N; i++) gel(z,i) = gen_0;
828 13407775 : return z;
829 : }
830 : GEN
831 1355167 : Rg_to_RgC(GEN x, long N)
832 : {
833 1355167 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
834 : }
835 :
836 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
837 : GEN
838 301943 : RgM_to_RgXV(GEN x, long v)
839 1284495 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
840 : GEN
841 7202 : RgM_to_RgXV_reverse(GEN x, long v)
842 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
843 :
844 : /* matrix whose entries are given by the coeffs of the polynomials in
845 : * vector v (considered as degree n-1 polynomials) */
846 : GEN
847 335583 : RgV_to_RgM(GEN x, long n)
848 1686221 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
849 :
850 : GEN
851 78528 : RgXV_to_RgM(GEN x, long n)
852 399385 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
853 :
854 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
855 : GEN
856 23469 : RgM_to_RgXX(GEN x, long v,long w)
857 : {
858 23469 : long j, lx = lg(x);
859 23469 : GEN y = cgetg(lx+1, t_POL);
860 23469 : y[1] = evalsigne(1) | evalvarn(v);
861 23469 : y++;
862 130426 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
863 23469 : return normalizepol_lg(--y, lx+1);
864 : }
865 :
866 : /* matrix whose entries are given by the coeffs of the polynomial v in
867 : * two variables (considered as degree n-1 polynomials) */
868 : GEN
869 322 : RgXX_to_RgM(GEN v, long n)
870 : {
871 322 : long j, N = lg(v)-1;
872 322 : GEN y = cgetg(N, t_MAT);
873 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
874 322 : return y;
875 : }
876 :
877 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
878 : GEN
879 32767 : RgXY_swapspec(GEN x, long n, long w, long nx)
880 : {
881 32767 : long j, ly = n+3;
882 32767 : GEN y = cgetg(ly, t_POL);
883 32767 : y[1] = evalsigne(1);
884 401015 : for (j=2; j<ly; j++)
885 : {
886 : long k;
887 368248 : GEN a = cgetg(nx+2,t_POL);
888 368248 : a[1] = evalsigne(1) | evalvarn(w);
889 2028272 : for (k=0; k<nx; k++)
890 : {
891 1660024 : GEN xk = gel(x,k);
892 1660024 : if (typ(xk)==t_POL && varn(xk)==w)
893 1565649 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
894 : else
895 94375 : gel(a,k+2) = j==2 ? xk: gen_0;
896 : }
897 368248 : gel(y,j) = normalizepol_lg(a, nx+2);
898 : }
899 32767 : return normalizepol_lg(y,ly);
900 : }
901 :
902 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
903 : GEN
904 2051 : RgXY_swap(GEN x, long n, long w)
905 : {
906 2051 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
907 2051 : setvarn(z, varn(x)); return z;
908 : }
909 :
910 : long
911 1387 : RgXY_degreex(GEN b)
912 : {
913 1387 : long deg = 0, i;
914 1387 : if (!signe(b)) return -1;
915 15592 : for (i = 2; i < lg(b); ++i)
916 : {
917 14205 : GEN bi = gel(b, i);
918 14205 : if (typ(bi) == t_POL)
919 13742 : deg = maxss(deg, degpol(bi));
920 : }
921 1387 : return deg;
922 : }
923 :
924 : GEN
925 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
926 :
927 : /* return (x % X^n). Shallow */
928 : GEN
929 7943799 : RgXn_red_shallow(GEN a, long n)
930 : {
931 7943799 : long i, L = n+2, l = lg(a);
932 : GEN b;
933 7943799 : if (L >= l) return a; /* deg(x) < n */
934 5751549 : b = cgetg(L, t_POL); b[1] = a[1];
935 37323242 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
936 5751550 : return normalizepol_lg(b,L);
937 : }
938 :
939 : GEN
940 483 : RgXnV_red_shallow(GEN x, long n)
941 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
942 :
943 : /* return (x * X^n). Shallow */
944 : GEN
945 175575247 : RgX_shift_shallow(GEN a, long n)
946 : {
947 175575247 : long i, l = lg(a);
948 : GEN b;
949 175575247 : if (l == 2 || !n) return a;
950 109159181 : l += n;
951 109159181 : if (n < 0)
952 : {
953 54203047 : if (l <= 2) return pol_0(varn(a));
954 52728950 : b = cgetg(l, t_POL); b[1] = a[1];
955 52729696 : a -= n;
956 163991726 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
957 : } else {
958 54956134 : b = cgetg(l, t_POL); b[1] = a[1];
959 54959824 : a -= n; n += 2;
960 119453934 : for (i=2; i<n; i++) gel(b,i) = gen_0;
961 211966402 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
962 : }
963 107689520 : return b;
964 : }
965 : /* return (x * X^n). */
966 : GEN
967 1578029 : RgX_shift(GEN a, long n)
968 : {
969 1578029 : long i, l = lg(a);
970 : GEN b;
971 1578029 : if (l == 2 || !n) return RgX_copy(a);
972 1577077 : l += n;
973 1577077 : if (n < 0)
974 : {
975 1442 : if (l <= 2) return pol_0(varn(a));
976 1372 : b = cgetg(l, t_POL); b[1] = a[1];
977 1372 : a -= n;
978 9534 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
979 : } else {
980 1575635 : b = cgetg(l, t_POL); b[1] = a[1];
981 1575635 : a -= n; n += 2;
982 3979076 : for (i=2; i<n; i++) gel(b,i) = gen_0;
983 4346090 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
984 : }
985 1577007 : return b;
986 : }
987 :
988 : GEN
989 317037 : RgX_rotate_shallow(GEN P, long k, long p)
990 : {
991 317037 : long i, l = lgpol(P);
992 : GEN r;
993 317037 : if (signe(P)==0)
994 1365 : return pol_0(varn(P));
995 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
996 2100644 : for(i=0; i<p; i++)
997 : {
998 1784972 : long s = 2+(i+k)%p;
999 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
1000 : }
1001 315672 : return RgX_renormalize(r);
1002 : }
1003 :
1004 : GEN
1005 2997279 : RgX_mulXn(GEN x, long d)
1006 : {
1007 : pari_sp av;
1008 : GEN z;
1009 : long v;
1010 2997279 : if (d >= 0) return RgX_shift(x, d);
1011 1472517 : d = -d;
1012 1472517 : v = RgX_val(x);
1013 1472517 : if (v >= d) return RgX_shift(x, -d);
1014 1472503 : av = avma;
1015 1472503 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1016 1472503 : return gerepileupto(av, z);
1017 : }
1018 :
1019 : long
1020 588 : RgXV_maxdegree(GEN x)
1021 : {
1022 588 : long d = -1, i, l = lg(x);
1023 4494 : for (i = 1; i < l; i++)
1024 3906 : d = maxss(d, degpol(gel(x,i)));
1025 588 : return d;
1026 : }
1027 :
1028 : long
1029 3667765 : RgX_val(GEN x)
1030 : {
1031 3667765 : long i, lx = lg(x);
1032 3667765 : if (lx == 2) return LONG_MAX;
1033 4619541 : for (i = 2; i < lx; i++)
1034 4619485 : if (!isexactzero(gel(x,i))) break;
1035 3656838 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1036 3656782 : return i - 2;
1037 : }
1038 : long
1039 80770910 : RgX_valrem(GEN x, GEN *Z)
1040 : {
1041 80770910 : long v, i, lx = lg(x);
1042 80770910 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1043 125181214 : for (i = 2; i < lx; i++)
1044 125183053 : if (!isexactzero(gel(x,i))) break;
1045 : /* possible with nonrational zeros */
1046 80771176 : if (i == lx)
1047 : {
1048 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1049 14 : return LONG_MAX;
1050 : }
1051 80771162 : v = i - 2;
1052 80771162 : *Z = RgX_shift_shallow(x, -v);
1053 80784021 : return v;
1054 : }
1055 : long
1056 851993 : RgX_valrem_inexact(GEN x, GEN *Z)
1057 : {
1058 : long v;
1059 851993 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1060 868980 : for (v = 0;; v++)
1061 868980 : if (!gequal0(gel(x,2+v))) break;
1062 851986 : if (Z) *Z = RgX_shift_shallow(x, -v);
1063 851986 : return v;
1064 : }
1065 :
1066 : GEN
1067 67683 : RgXQC_red(GEN x, GEN T)
1068 412979 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1069 :
1070 : GEN
1071 1211 : RgXQV_red(GEN x, GEN T)
1072 27832 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1073 :
1074 : GEN
1075 12999 : RgXQM_red(GEN x, GEN T)
1076 80682 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1077 :
1078 : GEN
1079 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1080 : {
1081 322 : return RgXQM_red(RgM_mul(P, Q), T);
1082 : }
1083 :
1084 : GEN
1085 508195 : RgXQX_red(GEN P, GEN T)
1086 : {
1087 508195 : long i, l = lg(P);
1088 508195 : GEN Q = cgetg(l, t_POL);
1089 508195 : Q[1] = P[1];
1090 2626824 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1091 508195 : return normalizepol_lg(Q, l);
1092 : }
1093 :
1094 : GEN
1095 825254 : RgX_deriv(GEN x)
1096 : {
1097 825254 : long i,lx = lg(x)-1;
1098 : GEN y;
1099 :
1100 825254 : if (lx<3) return pol_0(varn(x));
1101 822153 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1102 3639453 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1103 822146 : y[1] = x[1]; return normalizepol_lg(y,i);
1104 : }
1105 :
1106 : GEN
1107 2583662 : RgX_recipspec_shallow(GEN x, long l, long n)
1108 : {
1109 : long i;
1110 2583662 : GEN z = cgetg(n+2,t_POL);
1111 2583673 : z[1] = 0; z += 2;
1112 144302472 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1113 2833408 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1114 2583673 : return normalizepol_lg(z-2,n+2);
1115 : }
1116 :
1117 : GEN
1118 642832 : RgXn_recip_shallow(GEN P, long n)
1119 : {
1120 642832 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1121 642844 : setvarn(Q, varn(P));
1122 642844 : return Q;
1123 : }
1124 :
1125 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1126 : GEN
1127 31367 : RgX_recip(GEN x)
1128 : {
1129 : long lx, i, j;
1130 31367 : GEN y = cgetg_copy(x, &lx);
1131 274442 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1132 31367 : return normalizepol_lg(y,lx);
1133 : }
1134 : /* shallow version */
1135 : GEN
1136 59388 : RgX_recip_shallow(GEN x)
1137 : {
1138 : long lx, i, j;
1139 59388 : GEN y = cgetg_copy(x, &lx);
1140 356118 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1141 59388 : return normalizepol_lg(y,lx);
1142 : }
1143 :
1144 : GEN
1145 5112458 : RgX_recip_i(GEN x)
1146 : {
1147 : long lx, i, j;
1148 5112458 : GEN y = cgetg_copy(x, &lx);
1149 25942409 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1150 5112450 : return y;
1151 : }
1152 : /*******************************************************************/
1153 : /* */
1154 : /* ADDITION / SUBTRACTION */
1155 : /* */
1156 : /*******************************************************************/
1157 : /* same variable */
1158 : GEN
1159 144725281 : RgX_add(GEN x, GEN y)
1160 : {
1161 144725281 : long i, lx = lg(x), ly = lg(y);
1162 : GEN z;
1163 144725281 : if (ly <= lx) {
1164 130475040 : z = cgetg(lx,t_POL); z[1] = x[1];
1165 508052723 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1166 190222182 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1167 130452916 : z = normalizepol_lg(z, lx);
1168 : } else {
1169 14250241 : z = cgetg(ly,t_POL); z[1] = y[1];
1170 54481878 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1171 41633669 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1172 14251697 : z = normalizepol_lg(z, ly);
1173 : }
1174 144718142 : return z;
1175 : }
1176 : GEN
1177 72768852 : RgX_sub(GEN x, GEN y)
1178 : {
1179 72768852 : long i, lx = lg(x), ly = lg(y);
1180 : GEN z;
1181 72768852 : if (ly <= lx) {
1182 38938241 : z = cgetg(lx,t_POL); z[1] = x[1];
1183 188182991 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1184 66415136 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1185 38917770 : z = normalizepol_lg(z, lx);
1186 : } else {
1187 33830611 : z = cgetg(ly,t_POL); z[1] = y[1];
1188 125718380 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1189 74936478 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1190 33798402 : z = normalizepol_lg(z, ly);
1191 : }
1192 72761250 : return z;
1193 : }
1194 : GEN
1195 8820843 : RgX_neg(GEN x)
1196 53144957 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1197 :
1198 : GEN
1199 44377720 : RgX_Rg_add(GEN y, GEN x)
1200 : {
1201 : GEN z;
1202 44377720 : long lz = lg(y), i;
1203 44377720 : if (lz == 2) return scalarpol(x,varn(y));
1204 22633853 : z = cgetg(lz,t_POL); z[1] = y[1];
1205 22633809 : gel(z,2) = gadd(gel(y,2),x);
1206 78671734 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1207 : /* probably useless unless lz = 3, but cannot be skipped if y is
1208 : * an inexact 0 */
1209 22633728 : return normalizepol_lg(z,lz);
1210 : }
1211 : GEN
1212 65071 : RgX_Rg_add_shallow(GEN y, GEN x)
1213 : {
1214 : GEN z;
1215 65071 : long lz = lg(y), i;
1216 65071 : if (lz == 2) return scalarpol(x,varn(y));
1217 65071 : z = cgetg(lz,t_POL); z[1] = y[1];
1218 65071 : gel(z,2) = gadd(gel(y,2),x);
1219 130289 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1220 65070 : return normalizepol_lg(z,lz);
1221 : }
1222 : GEN
1223 1222886 : RgX_Rg_sub(GEN y, GEN x)
1224 : {
1225 : GEN z;
1226 1222886 : long lz = lg(y), i;
1227 1222886 : if (lz == 2)
1228 : { /* scalarpol(gneg(x),varn(y)) optimized */
1229 922867 : long v = varn(y);
1230 922867 : if (isrationalzero(x)) return pol_0(v);
1231 2235 : z = cgetg(3,t_POL);
1232 2235 : z[1] = gequal0(x)? evalvarn(v)
1233 2235 : : evalvarn(v) | evalsigne(1);
1234 2235 : gel(z,2) = gneg(x); return z;
1235 : }
1236 300019 : z = cgetg(lz,t_POL); z[1] = y[1];
1237 300019 : gel(z,2) = gsub(gel(y,2),x);
1238 686831 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1239 300020 : return normalizepol_lg(z,lz);
1240 : }
1241 : GEN
1242 4951499 : Rg_RgX_sub(GEN x, GEN y)
1243 : {
1244 : GEN z;
1245 4951499 : long lz = lg(y), i;
1246 4951499 : if (lz == 2) return scalarpol(x,varn(y));
1247 4881105 : z = cgetg(lz,t_POL); z[1] = y[1];
1248 4881112 : gel(z,2) = gsub(x, gel(y,2));
1249 7332772 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1250 4881046 : return normalizepol_lg(z,lz);
1251 : }
1252 : /*******************************************************************/
1253 : /* */
1254 : /* KARATSUBA MULTIPLICATION */
1255 : /* */
1256 : /*******************************************************************/
1257 : #if 0
1258 : /* to debug Karatsuba-like routines */
1259 : GEN
1260 : zx_debug_spec(GEN x, long nx)
1261 : {
1262 : GEN z = cgetg(nx+2,t_POL);
1263 : long i;
1264 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1265 : z[1] = evalsigne(1); return z;
1266 : }
1267 :
1268 : GEN
1269 : RgX_debug_spec(GEN x, long nx)
1270 : {
1271 : GEN z = cgetg(nx+2,t_POL);
1272 : long i;
1273 : for (i=0; i<nx; i++) z[i+2] = x[i];
1274 : z[1] = evalsigne(1); return z;
1275 : }
1276 : #endif
1277 :
1278 : /* generic multiplication */
1279 : GEN
1280 8907612 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1281 : {
1282 : GEN z, t;
1283 : long i;
1284 8907612 : if (nx == ny) {
1285 1536117 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1286 4953086 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1287 1536108 : return normalizepol_lg(z, nx+2);
1288 : }
1289 7371495 : if (ny < nx) {
1290 7186316 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1291 26452544 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1292 17172130 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1293 7185819 : return normalizepol_lg(z, nx+2);
1294 : } else {
1295 185179 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1296 3671822 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1297 447113 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1298 185192 : return normalizepol_lg(z, ny+2);
1299 : }
1300 : }
1301 : GEN
1302 222287 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1303 : {
1304 : GEN z, t;
1305 : long i;
1306 222287 : if (nx == ny) {
1307 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1308 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1309 12824 : return normalizepol_lg(z, nx+2);
1310 : }
1311 209463 : if (ny < nx) {
1312 207685 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1313 3724613 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1314 2371419 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1315 207685 : return normalizepol_lg(z, nx+2);
1316 : } else {
1317 1778 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1318 330904 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1319 12222 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1320 1778 : return normalizepol_lg(z, ny+2);
1321 : }
1322 : }
1323 :
1324 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1325 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1326 : * t_VECSMALL, to hold this, which can be left on stack: gerepile
1327 : * will not crash on it. The returned vector itself is not a proper GEN,
1328 : * we access the coefficients as x[i], i = 0..deg(x) */
1329 : static GEN
1330 81693405 : RgXspec_kill0(GEN x, long lx)
1331 : {
1332 81693405 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1333 : long i;
1334 246566013 : for (i=0; i <lx; i++)
1335 : {
1336 164872592 : GEN c = gel(x,i);
1337 164872592 : z[i] = (long)(isrationalzero(c)? NULL: c);
1338 : }
1339 81693421 : return z;
1340 : }
1341 :
1342 : INLINE GEN
1343 111052089 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1344 : {
1345 111052089 : pari_sp av = avma;
1346 111052089 : GEN s = NULL;
1347 : long i;
1348 :
1349 391647123 : for (i=a; i<b; i++)
1350 280604418 : if (gel(y,i) && gel(x,-i))
1351 : {
1352 202279837 : GEN t = gmul(gel(y,i), gel(x,-i));
1353 202274729 : s = s? gadd(s, t): t;
1354 : }
1355 111042705 : return s? gerepileupto(av, s): gen_0;
1356 : }
1357 :
1358 : /* assume nx >= ny > 0, return x * y * t^v */
1359 : static GEN
1360 34199818 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1361 : {
1362 : long i, lz, nz;
1363 : GEN z;
1364 :
1365 34199818 : x = RgXspec_kill0(x,nx);
1366 34199817 : y = RgXspec_kill0(y,ny);
1367 34199826 : lz = nx + ny + 1; nz = lz-2;
1368 34199826 : lz += v;
1369 34199826 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1370 72776981 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1371 83191328 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1372 56628052 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1373 48991583 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1374 34199189 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1375 : }
1376 :
1377 : /* return (x * X^d) + y. Assume d > 0 */
1378 : GEN
1379 10066003 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1380 : {
1381 : GEN x, y, xd, yd, zd;
1382 : long a, lz, nx, ny;
1383 :
1384 10066003 : if (!signe(x0)) return y0;
1385 9431179 : ny = lgpol(y0);
1386 9431179 : nx = lgpol(x0);
1387 9431299 : zd = (GEN)avma;
1388 9431299 : x = x0 + 2; y = y0 + 2; a = ny-d;
1389 9431299 : if (a <= 0)
1390 : {
1391 1479426 : lz = nx+d+2;
1392 1479426 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1393 3397392 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1394 1479428 : x = zd + a;
1395 1497125 : while (zd > x) gel(--zd,0) = gen_0;
1396 : }
1397 : else
1398 : {
1399 7951873 : xd = new_chunk(d); yd = y+d;
1400 7951870 : x = RgX_addspec_shallow(x,yd, nx,a);
1401 7951761 : lz = (a>nx)? ny+2: lg(x)+d;
1402 39624802 : x += 2; while (xd > x) *--zd = *--xd;
1403 : }
1404 22973665 : while (yd > y) *--zd = *--yd;
1405 9431189 : *--zd = x0[1];
1406 9431189 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1407 : }
1408 : GEN
1409 514757 : RgX_addmulXn(GEN x0, GEN y0, long d)
1410 : {
1411 : GEN x, y, xd, yd, zd;
1412 : long a, lz, nx, ny;
1413 :
1414 514757 : if (!signe(x0)) return RgX_copy(y0);
1415 513973 : nx = lgpol(x0);
1416 513973 : ny = lgpol(y0);
1417 513973 : zd = (GEN)avma;
1418 513973 : x = x0 + 2; y = y0 + 2; a = ny-d;
1419 513973 : if (a <= 0)
1420 : {
1421 291686 : lz = nx+d+2;
1422 291686 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1423 4293251 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1424 291686 : x = zd + a;
1425 757796 : while (zd > x) gel(--zd,0) = gen_0;
1426 : }
1427 : else
1428 : {
1429 222287 : xd = new_chunk(d); yd = y+d;
1430 222287 : x = RgX_addspec(x,yd, nx,a);
1431 222287 : lz = (a>nx)? ny+2: lg(x)+d;
1432 8415473 : x += 2; while (xd > x) *--zd = *--xd;
1433 : }
1434 2614898 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1435 513973 : *--zd = x0[1];
1436 513973 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1437 : }
1438 :
1439 : /* return x * y mod t^n */
1440 : static GEN
1441 6627959 : RgXn_mul_basecase(GEN x, GEN y, long n)
1442 : {
1443 6627959 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1444 : GEN z;
1445 6627959 : if (lx < 0) return pol_0(varn(x));
1446 6627959 : if (ly < 0) return pol_0(varn(x));
1447 6627959 : z = cgetg(lz, t_POL) + 2;
1448 6627959 : x+=2; if (lx > n) lx = n;
1449 6627959 : y+=2; if (ly > n) ly = n;
1450 6627959 : z[-1] = x[-1];
1451 6627959 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1452 6627959 : x = RgXspec_kill0(x, lx);
1453 6627959 : y = RgXspec_kill0(y, ly);
1454 : /* x:y:z [i] = term of degree i */
1455 26223885 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1456 11826184 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1457 6674031 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1458 6627959 : return normalizepol_lg(z - 2, lz);
1459 : }
1460 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1461 : static GEN
1462 9414318 : RgXn_mul2(GEN f, GEN g, long n)
1463 : {
1464 9414318 : pari_sp av = avma;
1465 : GEN fe,fo, ge,go, l,h,m;
1466 : long n0, n1;
1467 9414318 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1468 6663449 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1469 35490 : n0 = n>>1; n1 = n-n0;
1470 35490 : RgX_even_odd(f, &fe, &fo);
1471 35490 : RgX_even_odd(g, &ge, &go);
1472 35490 : l = RgXn_mul2(fe,ge,n1);
1473 35490 : h = RgXn_mul2(fo,go,n0);
1474 35490 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1475 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1476 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1477 35490 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1478 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1479 35490 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1480 35490 : m = RgX_inflate(m,2);
1481 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1482 35490 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1483 35490 : h = RgX_inflate(h,2);
1484 35490 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1485 35490 : return gerepileupto(av, h);
1486 : }
1487 : /* (f*g) \/ x^n */
1488 : static GEN
1489 1589447 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1490 : {
1491 1589447 : long d = degpol(f)+degpol(g) + 1 - n;
1492 : GEN h;
1493 1589447 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1494 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1495 : RgX_recip_i(g), d));
1496 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1497 : }
1498 :
1499 : /* (f*g) \/ x^n */
1500 : static GEN
1501 0 : RgX_sqrhigh_i2(GEN f, long n)
1502 : {
1503 0 : long d = 2*degpol(f)+ 1 - n;
1504 : GEN h;
1505 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1506 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1507 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1508 : }
1509 :
1510 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1511 : * b+2 were sent instead. na, nb = number of terms of a, b.
1512 : * Only c, c0, c1, c2 are genuine GEN.
1513 : */
1514 : GEN
1515 40265987 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1516 : {
1517 : GEN a0, c, c0;
1518 40265987 : long n0, n0a, i, v = 0;
1519 : pari_sp av;
1520 :
1521 65264295 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1522 61956664 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1523 40265963 : if (na < nb) swapspec(a,b, na,nb);
1524 40265963 : if (!nb) return pol_0(0);
1525 :
1526 34679022 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1527 479214 : RgX_shift_inplace_init(v);
1528 479228 : i = (na>>1); n0 = na-i; na = i;
1529 479228 : av = avma; a0 = a+n0; n0a = n0;
1530 1334087 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1531 :
1532 479228 : if (nb > n0)
1533 : {
1534 : GEN b0,c1,c2;
1535 : long n0b;
1536 :
1537 477869 : nb -= n0; b0 = b+n0; n0b = n0;
1538 1444999 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1539 477869 : c = RgX_mulspec(a,b,n0a,n0b);
1540 477869 : c0 = RgX_mulspec(a0,b0, na,nb);
1541 :
1542 477869 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1543 477869 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1544 :
1545 477869 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1546 477869 : c2 = RgX_sub(c1, RgX_add(c0,c));
1547 477869 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1548 : }
1549 : else
1550 : {
1551 1359 : c = RgX_mulspec(a,b,n0a,nb);
1552 1359 : c0 = RgX_mulspec(a0,b,na,nb);
1553 : }
1554 479228 : c0 = RgX_addmulXn(c0,c,n0);
1555 479228 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1556 : }
1557 :
1558 : INLINE GEN
1559 100082 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1560 : {
1561 100082 : pari_sp av = avma;
1562 100082 : GEN s = NULL;
1563 100082 : long j, l = (i+1)>>1;
1564 203924 : for (j=a; j<l; j++)
1565 : {
1566 103842 : GEN xj = gel(x,j), xx = gel(x,i-j);
1567 103842 : if (xj && xx)
1568 : {
1569 94707 : GEN t = gmul(xj, xx);
1570 94707 : s = s? gadd(s, t): t;
1571 : }
1572 : }
1573 100082 : if (s) s = gshift(s,1);
1574 100082 : if ((i&1) == 0)
1575 : {
1576 69041 : GEN t = gel(x, i>>1);
1577 69041 : if (t) {
1578 65982 : t = gsqr(t);
1579 65982 : s = s? gadd(s, t): t;
1580 : }
1581 : }
1582 100082 : return s? gerepileupto(av,s): gen_0;
1583 : }
1584 : static GEN
1585 38000 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1586 : {
1587 : long i, lz, nz;
1588 : GEN z;
1589 :
1590 38000 : if (!nx) return pol_0(0);
1591 38000 : x = RgXspec_kill0(x,nx);
1592 38000 : lz = (nx << 1) + 1, nz = lz-2;
1593 38000 : lz += v;
1594 38000 : z = cgetg(lz,t_POL) + 2;
1595 94602 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1596 107041 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1597 69041 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1598 38000 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1599 : }
1600 : /* return x^2 mod t^n */
1601 : static GEN
1602 0 : RgXn_sqr_basecase(GEN x, long n)
1603 : {
1604 0 : long i, lz = n+2, lx = lgpol(x);
1605 : GEN z;
1606 0 : if (lx < 0) return pol_0(varn(x));
1607 0 : z = cgetg(lz, t_POL);
1608 0 : z[1] = x[1];
1609 0 : x+=2; if (lx > n) lx = n;
1610 0 : x = RgXspec_kill0(x,lx);
1611 0 : z+=2;/* x:z [i] = term of degree i */
1612 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1613 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1614 0 : z -= 2; return normalizepol_lg(z, lz);
1615 : }
1616 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1617 : static GEN
1618 315 : RgXn_sqr2(GEN f, long n)
1619 : {
1620 315 : pari_sp av = avma;
1621 : GEN fe,fo, l,h,m;
1622 : long n0, n1;
1623 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1624 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1625 0 : n0 = n>>1; n1 = n-n0;
1626 0 : RgX_even_odd(f, &fe, &fo);
1627 0 : l = RgXn_sqr(fe,n1);
1628 0 : h = RgXn_sqr(fo,n0);
1629 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1630 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1631 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1632 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1633 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1634 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1635 0 : m = RgX_inflate(m,2);
1636 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1637 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1638 0 : h = RgX_inflate(h,2);
1639 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1640 0 : return gerepileupto(av, h);
1641 : }
1642 : GEN
1643 38039 : RgX_sqrspec(GEN a, long na)
1644 : {
1645 : GEN a0, c, c0, c1;
1646 38039 : long n0, n0a, i, v = 0;
1647 : pari_sp av;
1648 :
1649 66340 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1650 38039 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1651 39 : RgX_shift_inplace_init(v);
1652 39 : i = (na>>1); n0 = na-i; na = i;
1653 39 : av = avma; a0 = a+n0; n0a = n0;
1654 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1655 :
1656 39 : c = RgX_sqrspec(a,n0a);
1657 39 : c0 = RgX_sqrspec(a0,na);
1658 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1659 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1660 39 : c0 = RgX_addmulXn(c0,c,n0);
1661 39 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1662 : }
1663 :
1664 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1665 : GEN
1666 1719019 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1667 : {
1668 1719019 : GEN z = RgX_mul(A, B);
1669 1719005 : if (a < b)
1670 10376 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1671 1708629 : else if (a > b)
1672 1104945 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1673 : else
1674 603684 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1675 1719006 : return z;
1676 : }
1677 :
1678 : GEN
1679 38829700 : RgX_mul_i(GEN x, GEN y)
1680 : {
1681 38829700 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1682 38829289 : setvarn(z, varn(x)); return z;
1683 : }
1684 :
1685 : GEN
1686 37961 : RgX_sqr_i(GEN x)
1687 : {
1688 37961 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1689 37961 : setvarn(z,varn(x)); return z;
1690 : }
1691 :
1692 : /*******************************************************************/
1693 : /* */
1694 : /* DIVISION */
1695 : /* */
1696 : /*******************************************************************/
1697 : GEN
1698 7177428 : RgX_Rg_divexact(GEN x, GEN y) {
1699 7177428 : long i, lx = lg(x);
1700 : GEN z;
1701 7177428 : if (lx == 2) return gcopy(x);
1702 7126465 : switch(typ(y))
1703 : {
1704 6965197 : case t_INT:
1705 6965197 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1706 6589776 : break;
1707 5243 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1708 : }
1709 6745801 : z = cgetg(lx, t_POL); z[1] = x[1];
1710 34443321 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1711 6744851 : return z;
1712 : }
1713 : GEN
1714 30320883 : RgX_Rg_div(GEN x, GEN y) {
1715 30320883 : long i, lx = lg(x);
1716 : GEN z;
1717 30320883 : if (lx == 2) return gcopy(x);
1718 30048006 : switch(typ(y))
1719 : {
1720 20815471 : case t_INT:
1721 20815471 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1722 4192013 : break;
1723 6055 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1724 : }
1725 13418493 : z = cgetg(lx, t_POL); z[1] = x[1];
1726 50019255 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1727 13418077 : return normalizepol_lg(z, lx);
1728 : }
1729 : GEN
1730 39130 : RgX_normalize(GEN x)
1731 : {
1732 39130 : GEN z, d = NULL;
1733 39130 : long i, n = lg(x)-1;
1734 39130 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1735 39130 : if (i == 1) return pol_0(varn(x));
1736 39130 : if (i == n && isint1(d)) return x;
1737 17612 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1738 31612 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1739 17612 : gel(z,n) = Rg_get_1(d); return z;
1740 : }
1741 : GEN
1742 10458 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1743 : GEN
1744 250999 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1745 : {
1746 250999 : long l = lg(a), i;
1747 : GEN a0, z0, z;
1748 :
1749 250999 : if (l <= 3)
1750 : {
1751 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1752 0 : return pol_0(varn(a));
1753 : }
1754 250999 : z = cgetg(l-1, t_POL);
1755 251003 : z[1] = a[1];
1756 251003 : a0 = a + l-1;
1757 251003 : z0 = z + l-2; *z0 = *a0--;
1758 2954751 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1759 : {
1760 2703747 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1761 2703748 : gel(z0,0) = t;
1762 : }
1763 251004 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1764 251004 : return z;
1765 : }
1766 : /* Polynomial division x / y:
1767 : * if pr = ONLY_REM return remainder, otherwise return quotient
1768 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1769 : * if pr != NULL set *pr to remainder, as the last object on stack */
1770 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1771 : static GEN
1772 23511934 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1773 : {
1774 : pari_sp avy, av, av1;
1775 : long dx,dy,dz,i,j,sx,lr;
1776 : GEN z,p1,p2,rem,y_lead,mod,p;
1777 : GEN (*f)(GEN,GEN);
1778 :
1779 23511934 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1780 :
1781 23511934 : dy = degpol(y);
1782 23511929 : y_lead = gel(y,dy+2);
1783 23511929 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1784 : {
1785 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1786 0 : for (dy--; dy>=0; dy--)
1787 : {
1788 0 : y_lead = gel(y,dy+2);
1789 0 : if (!gequal0(y_lead)) break;
1790 : }
1791 : }
1792 23511915 : if (!dy) /* y is constant */
1793 : {
1794 6844 : if (pr == ONLY_REM) return pol_0(varn(x));
1795 6844 : z = RgX_Rg_div(x, y_lead);
1796 6844 : if (pr == ONLY_DIVIDES) return z;
1797 2630 : if (pr) *pr = pol_0(varn(x));
1798 2630 : return z;
1799 : }
1800 23505071 : dx = degpol(x);
1801 23505040 : if (dx < dy)
1802 : {
1803 3847509 : if (pr == ONLY_REM) return RgX_copy(x);
1804 348542 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1805 348521 : z = pol_0(varn(x));
1806 348521 : if (pr) *pr = RgX_copy(x);
1807 348521 : return z;
1808 : }
1809 :
1810 : /* x,y in R[X], y non constant */
1811 19657531 : av = avma;
1812 19657531 : p = NULL;
1813 19657531 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1814 : {
1815 130277 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1816 130277 : if (!z) return gc_NULL(av);
1817 130277 : z = FpX_to_mod(z, p);
1818 130277 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1819 71134 : return gerepileupto(av, z);
1820 59143 : *pr = FpX_to_mod(*pr, p);
1821 59143 : return gc_all(av, 2, &z, pr);
1822 : }
1823 19527372 : switch(typ(y_lead))
1824 : {
1825 0 : case t_REAL:
1826 0 : y_lead = ginv(y_lead);
1827 0 : f = gmul; mod = NULL;
1828 0 : break;
1829 1878 : case t_INTMOD:
1830 1878 : case t_POLMOD: y_lead = ginv(y_lead);
1831 1878 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1832 1878 : break;
1833 19525494 : default: if (gequal1(y_lead)) y_lead = NULL;
1834 19525485 : f = gdiv; mod = NULL;
1835 : }
1836 :
1837 19527363 : if (y_lead == NULL)
1838 17515507 : p2 = gel(x,dx+2);
1839 : else {
1840 : for(;;) {
1841 2011856 : p2 = f(gel(x,dx+2),y_lead);
1842 2011856 : p2 = simplify_shallow(p2);
1843 2011856 : if (!isexactzero(p2) || (--dx < 0)) break;
1844 : }
1845 2011856 : if (dx < dy) /* leading coeff of x was in fact zero */
1846 : {
1847 0 : if (pr == ONLY_DIVIDES) {
1848 0 : set_avma(av);
1849 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1850 : }
1851 0 : if (pr == ONLY_REM)
1852 : {
1853 0 : if (dx < 0)
1854 0 : return gc_GEN(av, scalarpol(p2, varn(x)));
1855 : else
1856 : {
1857 : GEN t;
1858 0 : set_avma(av);
1859 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1860 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1861 0 : return t;
1862 : }
1863 : }
1864 0 : if (pr) /* cf ONLY_REM above */
1865 : {
1866 0 : if (dx < 0)
1867 : {
1868 0 : p2 = gclone(p2);
1869 0 : set_avma(av);
1870 0 : z = pol_0(varn(x));
1871 0 : x = scalarpol(p2, varn(x));
1872 0 : gunclone(p2);
1873 : }
1874 : else
1875 : {
1876 : GEN t;
1877 0 : set_avma(av);
1878 0 : z = pol_0(varn(x));
1879 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1880 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1881 0 : x = t;
1882 : }
1883 0 : *pr = x;
1884 : }
1885 : else
1886 : {
1887 0 : set_avma(av);
1888 0 : z = pol_0(varn(x));
1889 : }
1890 0 : return z;
1891 : }
1892 : }
1893 : /* dx >= dy */
1894 19527363 : avy = avma;
1895 19527363 : dz = dx-dy;
1896 19527363 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1897 19527346 : x += 2;
1898 19527346 : z += 2;
1899 19527346 : y += 2;
1900 19527346 : gel(z,dz) = gcopy(p2);
1901 :
1902 54084384 : for (i=dx-1; i>=dy; i--)
1903 : {
1904 34557191 : av1=avma; p1=gel(x,i);
1905 1140247646 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1906 34551982 : if (y_lead) p1 = simplify(f(p1,y_lead));
1907 :
1908 34551983 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1909 : else
1910 24562875 : p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1911 34556635 : gel(z,i-dy) = p1;
1912 : }
1913 19527193 : if (!pr) return gerepileupto(av,z-2);
1914 :
1915 12549997 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1916 12907028 : for (sx=0; ; i--)
1917 : {
1918 12907028 : p1 = gel(x,i);
1919 : /* we always enter this loop at least once */
1920 32217753 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1921 12905795 : if (mod && avma==av1) p1 = gmul(p1,mod);
1922 12905795 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1923 3576642 : if (!isexactzero(p1)) break;
1924 3426506 : if (!i) break;
1925 356989 : set_avma(av1);
1926 : }
1927 12549143 : if (pr == ONLY_DIVIDES)
1928 : {
1929 7931 : if (sx) return gc_NULL(av);
1930 7910 : set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1931 : }
1932 12541212 : lr=i+3; rem -= lr;
1933 12541212 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1934 12498923 : else p1 = gerepileupto((pari_sp)rem,p1);
1935 12542000 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1936 12542000 : rem[1] = z[-1];
1937 12542000 : rem += 2;
1938 12542000 : gel(rem,i) = p1;
1939 17281781 : for (i--; i>=0; i--)
1940 : {
1941 4739783 : av1=avma; p1 = gel(x,i);
1942 13634007 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1943 4739263 : if (mod && avma==av1) p1 = gmul(p1,mod);
1944 4739530 : gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1945 : }
1946 12541998 : rem -= 2;
1947 12541998 : if (!sx) (void)normalizepol_lg(rem, lr);
1948 12542203 : if (pr == ONLY_REM) return gerepileupto(av,rem);
1949 6760644 : z -= 2;
1950 : {
1951 6760644 : GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1952 6760644 : gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1953 : }
1954 : }
1955 :
1956 : GEN
1957 14231463 : RgX_divrem(GEN x, GEN y, GEN *pr)
1958 : {
1959 14231463 : if (pr == ONLY_REM) return RgX_rem(x, y);
1960 14231463 : return RgX_divrem_i(x, y, pr);
1961 : }
1962 :
1963 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1964 : GEN
1965 155548 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1966 : {
1967 : long vx, dx, dy, dz, i, j, sx, lr;
1968 : pari_sp av0, av, tetpil;
1969 : GEN z,p1,rem,lead;
1970 :
1971 155548 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1972 155548 : vx = varn(x);
1973 155548 : dx = degpol(x);
1974 155548 : dy = degpol(y);
1975 155548 : if (dx < dy)
1976 : {
1977 41112 : if (pr)
1978 : {
1979 41084 : av0 = avma; x = RgXQX_red(x, T);
1980 41084 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1981 41077 : if (pr == ONLY_REM) return x;
1982 0 : *pr = x;
1983 : }
1984 28 : return pol_0(vx);
1985 : }
1986 114436 : lead = leading_coeff(y);
1987 114436 : if (!dy) /* y is constant */
1988 : {
1989 602 : if (pr && pr != ONLY_DIVIDES)
1990 : {
1991 0 : if (pr == ONLY_REM) return pol_0(vx);
1992 0 : *pr = pol_0(vx);
1993 : }
1994 602 : if (gequal1(lead)) return RgX_copy(x);
1995 0 : av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
1996 0 : return gerepile(av0,tetpil,RgXQX_red(x,T));
1997 : }
1998 113834 : av0 = avma; dz = dx-dy;
1999 113834 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
2000 113834 : set_avma(av0);
2001 113834 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2002 113834 : x += 2; y += 2; z += 2;
2003 :
2004 113834 : p1 = gel(x,dx); av = avma;
2005 113834 : gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2006 579254 : for (i=dx-1; i>=dy; i--)
2007 : {
2008 465420 : av=avma; p1=gel(x,i);
2009 2412457 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2010 465383 : if (lead) p1 = gmul(grem(p1, T), lead);
2011 465386 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
2012 : }
2013 113834 : if (!pr) { guncloneNULL(lead); return z-2; }
2014 :
2015 112560 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2016 189216 : for (sx=0; ; i--)
2017 : {
2018 189216 : p1 = gel(x,i);
2019 654709 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2020 189213 : tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2021 103738 : if (!i) break;
2022 76656 : set_avma(av);
2023 : }
2024 112559 : if (pr == ONLY_DIVIDES)
2025 : {
2026 27065 : guncloneNULL(lead);
2027 27065 : if (sx) return gc_NULL(av0);
2028 24506 : return gc_const((pari_sp)rem, z-2);
2029 : }
2030 85494 : lr=i+3; rem -= lr;
2031 85494 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2032 85494 : rem[1] = z[-1];
2033 85494 : p1 = gerepile((pari_sp)rem,tetpil,p1);
2034 85494 : rem += 2; gel(rem,i) = p1;
2035 186658 : for (i--; i>=0; i--)
2036 : {
2037 101164 : av=avma; p1 = gel(x,i);
2038 337033 : for (j=0; j<=i && j<=dz; j++)
2039 235869 : p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2040 101164 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
2041 : }
2042 85494 : rem -= 2;
2043 85494 : guncloneNULL(lead);
2044 85494 : if (!sx) (void)normalizepol_lg(rem, lr);
2045 85494 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
2046 161 : *pr = rem; return z-2;
2047 : }
2048 :
2049 : /*******************************************************************/
2050 : /* */
2051 : /* PSEUDO-DIVISION */
2052 : /* */
2053 : /*******************************************************************/
2054 : INLINE GEN
2055 4042863 : rem(GEN c, GEN T)
2056 : {
2057 4042863 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2058 4042865 : return c;
2059 : }
2060 :
2061 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2062 : int
2063 17335 : ZXQX_dvd(GEN x, GEN y, GEN T)
2064 : {
2065 : long dx, dy, i, T_ismonic;
2066 17335 : pari_sp av = avma, av2;
2067 : GEN y_lead;
2068 :
2069 17335 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2070 17335 : dy = degpol(y); y_lead = gel(y,dy+2);
2071 17335 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2072 : /* if monic, no point in using pseudo-division */
2073 17335 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2074 14794 : T_ismonic = gequal1(leading_coeff(T));
2075 14794 : dx = degpol(x);
2076 14794 : if (dx < dy) return !signe(x);
2077 14794 : (void)new_chunk(2);
2078 14794 : x = RgX_recip_i(x)+2;
2079 14794 : y = RgX_recip_i(y)+2;
2080 : /* pay attention to sparse divisors */
2081 30089 : for (i = 1; i <= dy; i++)
2082 15295 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2083 14794 : av2 = avma;
2084 : for (;;)
2085 73008 : {
2086 87802 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2087 87802 : x0 = gneg(x0);
2088 87803 : m = gcdii(cx, y0);
2089 87802 : if (!equali1(m))
2090 : {
2091 85300 : x0 = gdiv(x0, m);
2092 85301 : y0 = diviiexact(y0, m);
2093 85300 : if (equali1(y0)) y0 = NULL;
2094 : }
2095 179377 : for (i=1; i<=dy; i++)
2096 : {
2097 91575 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2098 91575 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2099 91576 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2100 91575 : gel(x,i) = c;
2101 : }
2102 688163 : for ( ; i<=dx; i++)
2103 : {
2104 600361 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2105 600361 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2106 600361 : gel(x,i) = c;
2107 : }
2108 102348 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2109 87802 : if (dx < dy) break;
2110 73008 : if (gc_needed(av2,1))
2111 : {
2112 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2113 28 : gerepilecoeffs(av2,x,dx+1);
2114 : }
2115 : }
2116 14794 : return gc_bool(av, dx < 0);
2117 : }
2118 :
2119 : /* T either NULL or a t_POL. */
2120 : GEN
2121 627563 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2122 : {
2123 627563 : long vx = varn(x), dx, dy, dz, i, lx, p;
2124 627563 : pari_sp av = avma, av2;
2125 : GEN y_lead;
2126 :
2127 627563 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2128 627563 : dy = degpol(y); y_lead = gel(y,dy+2);
2129 : /* if monic, no point in using pseudo-division */
2130 627563 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2131 579251 : dx = degpol(x);
2132 579251 : if (dx < dy) return RgX_copy(x);
2133 579244 : (void)new_chunk(2);
2134 579244 : x = RgX_recip_i(x)+2;
2135 579244 : y = RgX_recip_i(y)+2;
2136 : /* pay attention to sparse divisors */
2137 1951547 : for (i = 1; i <= dy; i++)
2138 1372304 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2139 579243 : dz = dx-dy; p = dz+1;
2140 579243 : av2 = avma;
2141 : for (;;)
2142 : {
2143 1246575 : gel(x,0) = gneg(gel(x,0)); p--;
2144 4192209 : for (i=1; i<=dy; i++)
2145 : {
2146 2945663 : GEN c = gmul(y_lead, gel(x,i));
2147 2945596 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2148 2945624 : gel(x,i) = rem(c, T);
2149 : }
2150 2182560 : for ( ; i<=dx; i++)
2151 : {
2152 935988 : GEN c = gmul(y_lead, gel(x,i));
2153 935981 : gel(x,i) = rem(c, T);
2154 : }
2155 1289399 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2156 1246572 : if (dx < dy) break;
2157 667332 : if (gc_needed(av2,1))
2158 : {
2159 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2160 0 : gerepilecoeffs(av2,x,dx+1);
2161 : }
2162 : }
2163 579240 : if (dx < 0) return pol_0(vx);
2164 579079 : lx = dx+3; x -= 2;
2165 579079 : x[0] = evaltyp(t_POL) | _evallg(lx);
2166 579079 : x[1] = evalsigne(1) | evalvarn(vx);
2167 579079 : x = RgX_recip_i(x);
2168 579081 : if (p)
2169 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2170 28674 : GEN t = y_lead;
2171 28674 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2172 0 : t = RgXQ_powu(t, p, T);
2173 : else
2174 28674 : t = gpowgs(t, p);
2175 88374 : for (i=2; i<lx; i++)
2176 : {
2177 59700 : GEN c = gmul(gel(x,i), t);
2178 59700 : gel(x,i) = rem(c,T);
2179 : }
2180 28674 : if (!T) return gerepileupto(av, x);
2181 : }
2182 550407 : return gc_GEN(av, x);
2183 : }
2184 :
2185 : GEN
2186 627563 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2187 :
2188 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2189 : GEN
2190 12056 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2191 : {
2192 12056 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2193 12056 : pari_sp av = avma, av2;
2194 : GEN z, r, ypow, y_lead;
2195 :
2196 12056 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2197 12056 : dy = degpol(y); y_lead = gel(y,dy+2);
2198 12056 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2199 7561 : dx = degpol(x);
2200 7561 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2201 7561 : if (dx == dy)
2202 : {
2203 98 : GEN x_lead = gel(x,lg(x)-1);
2204 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2205 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2206 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2207 98 : *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
2208 : }
2209 7463 : (void)new_chunk(2);
2210 7463 : x = RgX_recip_i(x)+2;
2211 7463 : y = RgX_recip_i(y)+2;
2212 : /* pay attention to sparse divisors */
2213 39000 : for (i = 1; i <= dy; i++)
2214 31537 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2215 7463 : dz = dx-dy; p = dz+1;
2216 7463 : lz = dz+3;
2217 7463 : z = cgetg(lz, t_POL);
2218 7463 : z[1] = evalsigne(1) | evalvarn(vx);
2219 28640 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2220 7463 : ypow = new_chunk(dz+1);
2221 7463 : gel(ypow,0) = gen_1;
2222 7463 : gel(ypow,1) = y_lead;
2223 13714 : for (i=2; i<=dz; i++)
2224 : {
2225 6251 : GEN c = gmul(gel(ypow,i-1), y_lead);
2226 6251 : gel(ypow,i) = rem(c,T);
2227 : }
2228 7463 : av2 = avma;
2229 7463 : for (iz=2;;)
2230 : {
2231 15605 : p--;
2232 15605 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2233 69915 : for (i=1; i<=dy; i++)
2234 : {
2235 54310 : GEN c = gmul(y_lead, gel(x,i));
2236 54310 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2237 54310 : gel(x,i) = rem(c, T);
2238 : }
2239 41058 : for ( ; i<=dx; i++)
2240 : {
2241 25453 : GEN c = gmul(y_lead, gel(x,i));
2242 25453 : gel(x,i) = rem(c,T);
2243 : }
2244 15605 : x++; dx--;
2245 21177 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2246 15605 : if (dx < dy) break;
2247 8142 : if (gc_needed(av2,1))
2248 : {
2249 0 : GEN X = x-2;
2250 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2251 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2252 0 : gerepileall(av2,2, &X, &z); x = X+2;
2253 : }
2254 : }
2255 14008 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2256 7463 : if (dx < 0)
2257 182 : x = pol_0(vx);
2258 : else
2259 : {
2260 7281 : lx = dx+3; x -= 2;
2261 7281 : x[0] = evaltyp(t_POL) | _evallg(lx);
2262 7281 : x[1] = evalsigne(1) | evalvarn(vx);
2263 7281 : x = RgX_recip_i(x);
2264 : }
2265 7463 : z = RgX_recip_i(z);
2266 7463 : r = x;
2267 7463 : if (p)
2268 : {
2269 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2270 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2271 : }
2272 7463 : *ptr = r; return gc_all(av, 2, &z, ptr);
2273 : }
2274 : GEN
2275 11818 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2276 11818 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2277 :
2278 : GEN
2279 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2280 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2281 : GEN
2282 748049253 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2283 : GEN
2284 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2285 : GEN
2286 26320 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2287 : GEN
2288 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2289 : GEN
2290 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2291 :
2292 : GEN
2293 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2294 :
2295 : GEN
2296 0 : RgXQX_powers(GEN P, long n, GEN T)
2297 : {
2298 0 : GEN v = cgetg(n+2, t_VEC);
2299 : long i;
2300 0 : gel(v, 1) = pol_1(varn(T));
2301 0 : if (n==0) return v;
2302 0 : gel(v, 2) = gcopy(P);
2303 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2304 0 : return v;
2305 : }
2306 :
2307 : static GEN
2308 623021 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2309 : static GEN
2310 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2311 : static GEN
2312 67310 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2313 : static GEN
2314 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2315 : static GEN
2316 276612 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2317 : static GEN
2318 1039955 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2319 : static GEN
2320 1018592 : _one(void *data) { return pol_1(varn((GEN)data)); }
2321 : static GEN
2322 1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2323 : static GEN
2324 723900 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2325 :
2326 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2327 : _mul, _sqr, _one, _zero };
2328 :
2329 : GEN
2330 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2331 : {
2332 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2333 : }
2334 :
2335 : GEN
2336 417188 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2337 : {
2338 417188 : int use_sqr = 2*degpol(x) >= degpol(T);
2339 417189 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2340 : }
2341 :
2342 : /* mod X^n */
2343 : struct modXn {
2344 : long v; /* varn(X) */
2345 : long n;
2346 : } ;
2347 : static GEN
2348 11893 : _sqrXn(void *data, GEN x) {
2349 11893 : struct modXn *S = (struct modXn*)data;
2350 11893 : return RgXn_sqr(x, S->n);
2351 : }
2352 : static GEN
2353 4528 : _mulXn(void *data, GEN x, GEN y) {
2354 4528 : struct modXn *S = (struct modXn*)data;
2355 4528 : return RgXn_mul(x,y, S->n);
2356 : }
2357 : static GEN
2358 1939 : _oneXn(void *data) {
2359 1939 : struct modXn *S = (struct modXn*)data;
2360 1939 : return pol_1(S->v);
2361 : }
2362 : static GEN
2363 0 : _zeroXn(void *data) {
2364 0 : struct modXn *S = (struct modXn*)data;
2365 0 : return pol_0(S->v);
2366 : }
2367 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2368 : _oneXn, _zeroXn };
2369 :
2370 : GEN
2371 483 : RgXn_powers(GEN x, long m, long n)
2372 : {
2373 483 : long d = degpol(x);
2374 483 : int use_sqr = (d<<1) >= n;
2375 : struct modXn S;
2376 483 : S.v = varn(x); S.n = n;
2377 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2378 : }
2379 :
2380 : GEN
2381 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2382 : {
2383 : struct modXn S;
2384 : long v;
2385 2286 : if (n == 0) return x;
2386 2286 : v = RgX_valrem(x, &x);
2387 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2388 2265 : S.v = varn(x); S.n = n;
2389 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2390 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2391 2265 : return x;
2392 : }
2393 : GEN
2394 0 : RgXn_powu(GEN x, ulong m, long n)
2395 : {
2396 : pari_sp av;
2397 0 : if (n == 0) return gcopy(x);
2398 0 : av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
2399 : }
2400 :
2401 : GEN
2402 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2403 : {
2404 : struct modXn S;
2405 966 : S.v = varn(gel(x,2)); S.n = n;
2406 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2407 : }
2408 :
2409 : GEN
2410 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2411 : {
2412 0 : int use_sqr = 2*degpol(x) >= n;
2413 : struct modXn S;
2414 0 : S.v = varn(x); S.n = n;
2415 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2416 : }
2417 :
2418 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2419 : GEN
2420 5159 : RgXn_eval(GEN Q, GEN x, long n)
2421 : {
2422 5159 : long d = degpol(x);
2423 : int use_sqr;
2424 : struct modXn S;
2425 5159 : if (d == 1 && isrationalzero(gel(x,2)))
2426 : {
2427 5152 : GEN y = RgX_unscale(Q, gel(x,3));
2428 5152 : setvarn(y, varn(x)); return y;
2429 : }
2430 7 : S.v = varn(x);
2431 7 : S.n = n;
2432 7 : use_sqr = (d<<1) >= n;
2433 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2434 : }
2435 :
2436 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2437 : static GEN
2438 2614673 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2439 : {
2440 2614673 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2441 2614673 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2442 : }
2443 :
2444 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2445 : static GEN
2446 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2447 : {
2448 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2449 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2450 : }
2451 :
2452 : static GEN
2453 2188941 : RgXn_div_gen(GEN g, GEN f, long e)
2454 : {
2455 : pari_sp av;
2456 : ulong mask;
2457 : GEN W, a;
2458 2188941 : long v = varn(f), n = 1;
2459 :
2460 2188941 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2461 2188934 : a = ginv(gel(f,2));
2462 2188932 : if (e == 1 && !g) return scalarpol(a, v);
2463 2184046 : else if (e == 2 && !g)
2464 : {
2465 : GEN b;
2466 836929 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2467 246248 : b = gneg(b);
2468 246251 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2469 246251 : return deg1pol(b, a, v);
2470 : }
2471 1347117 : av = avma;
2472 1347117 : W = scalarpol_shallow(a,v);
2473 1347117 : mask = quadratic_prec_mask(e);
2474 3955833 : while (mask > 1)
2475 : {
2476 : GEN u, fr;
2477 2608716 : long n2 = n;
2478 2608716 : n<<=1; if (mask & 1) n--;
2479 2608716 : mask >>= 1;
2480 2608716 : fr = RgXn_red_shallow(f, n);
2481 2608716 : if (mask>1 || !g)
2482 : {
2483 1324736 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2484 1324736 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2485 : }
2486 : else
2487 : {
2488 1283980 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2489 1283980 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2490 1283980 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2491 : }
2492 2608716 : if (gc_needed(av,2))
2493 : {
2494 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2495 0 : W = gerepileupto(av, W);
2496 : }
2497 : }
2498 1347117 : return W;
2499 : }
2500 :
2501 : static GEN
2502 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2503 : {
2504 : GEN r;
2505 119 : if (lgefint(p) == 3)
2506 : {
2507 119 : ulong pp = uel(p, 2);
2508 119 : if (pp == 2)
2509 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2510 : else
2511 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2512 : }
2513 : else
2514 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2515 119 : return FpX_to_mod(r, p);
2516 : }
2517 :
2518 : static GEN
2519 7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2520 : {
2521 7 : GEN r, T = RgX_to_FpX(pol, p);
2522 7 : if (signe(T) == 0) pari_err_OP("/", x, y);
2523 7 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2524 7 : return FpXQX_to_mod(r, T, p);
2525 : }
2526 :
2527 : static GEN
2528 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2529 : {
2530 : GEN r;
2531 91 : if (lgefint(p) == 3)
2532 : {
2533 91 : ulong pp = uel(p, 2);
2534 91 : if (pp == 2)
2535 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2536 : else
2537 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2538 : }
2539 : else
2540 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2541 91 : return FpX_to_mod(r, p);
2542 : }
2543 :
2544 : static GEN
2545 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2546 : {
2547 0 : GEN r, T = RgX_to_FpX(pol, p);
2548 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2549 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2550 0 : return FpXQX_to_mod(r, T, p);
2551 : }
2552 :
2553 : static GEN
2554 905054 : RgXn_inv_fast(GEN x, long e)
2555 : {
2556 : GEN p, pol;
2557 : long pa;
2558 905054 : long t = RgX_type(x,&p,&pol,&pa);
2559 905052 : switch(t)
2560 : {
2561 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2562 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2563 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2564 904961 : default: return NULL;
2565 : }
2566 : }
2567 :
2568 : static GEN
2569 1284106 : RgXn_div_fast(GEN x, GEN y, long e)
2570 : {
2571 : GEN p, pol;
2572 : long pa;
2573 1284106 : long t = RgX_type2(x,y,&p,&pol,&pa);
2574 1284106 : switch(t)
2575 : {
2576 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2577 7 : case RgX_type_code(t_POLMOD, t_INTMOD):
2578 7 : return RgXn_div_FpXQX(x, y, e, pol, p);
2579 1283980 : default: return NULL;
2580 : }
2581 : }
2582 :
2583 : GEN
2584 1284106 : RgXn_div_i(GEN g, GEN f, long e)
2585 : {
2586 1284106 : GEN h = RgXn_div_fast(g, f, e);
2587 1284106 : if (h) return h;
2588 1283980 : return RgXn_div_gen(g, f, e);
2589 : }
2590 :
2591 : GEN
2592 584443 : RgXn_div(GEN g, GEN f, long e)
2593 : {
2594 584443 : pari_sp av = avma;
2595 584443 : return gerepileupto(av, RgXn_div_i(g, f, e));
2596 : }
2597 :
2598 : GEN
2599 905054 : RgXn_inv_i(GEN f, long e)
2600 : {
2601 905054 : GEN h = RgXn_inv_fast(f, e);
2602 905052 : if (h) return h;
2603 904961 : return RgXn_div_gen(NULL, f, e);
2604 : }
2605 :
2606 : GEN
2607 807474 : RgXn_inv(GEN f, long e)
2608 : {
2609 807474 : pari_sp av = avma;
2610 807474 : return gerepileupto(av, RgXn_inv_i(f, e));
2611 : }
2612 :
2613 : /* intformal(x^n*S) / x^(n+1) */
2614 : static GEN
2615 56374 : RgX_integXn(GEN x, long n)
2616 114351 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2617 :
2618 : GEN
2619 52916 : RgXn_expint(GEN h, long e)
2620 : {
2621 52916 : pari_sp av = avma, av2;
2622 52916 : long v = varn(h), n;
2623 52916 : GEN f = pol_1(v), g;
2624 : ulong mask;
2625 :
2626 52916 : if (!signe(h)) return f;
2627 50424 : g = pol_1(v);
2628 50423 : n = 1; mask = quadratic_prec_mask(e);
2629 50423 : av2 = avma;
2630 56373 : for (;mask>1;)
2631 : {
2632 : GEN u, w;
2633 56373 : long n2 = n;
2634 56373 : n<<=1; if (mask & 1) n--;
2635 56373 : mask >>= 1;
2636 56373 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2637 56374 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2638 56374 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2639 56374 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2640 56373 : if (mask<=1) break;
2641 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2642 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2643 5950 : if (gc_needed(av2,2))
2644 : {
2645 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2646 0 : gerepileall(av2, 2, &f, &g);
2647 : }
2648 : }
2649 50423 : return gerepileupto(av, f);
2650 : }
2651 :
2652 : GEN
2653 0 : RgXn_exp(GEN h, long e)
2654 : {
2655 0 : long d = degpol(h);
2656 0 : if (d < 0) return pol_1(varn(h));
2657 0 : if (!d || !gequal0(gel(h,2)))
2658 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2659 0 : return RgXn_expint(RgX_deriv(h), e);
2660 : }
2661 :
2662 : GEN
2663 154 : RgXn_reverse(GEN f, long e)
2664 : {
2665 154 : pari_sp av = avma, av2;
2666 : ulong mask;
2667 : GEN fi, a, df, W, an;
2668 154 : long v = varn(f), n=1;
2669 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2670 0 : pari_err_INV("serreverse",f);
2671 154 : fi = ginv(gel(f,3));
2672 154 : a = deg1pol_shallow(fi,gen_0,v);
2673 154 : if (e <= 2) return gc_GEN(av, a);
2674 133 : W = scalarpol(fi,v);
2675 133 : df = RgX_deriv(f);
2676 133 : mask = quadratic_prec_mask(e);
2677 133 : av2 = avma;
2678 616 : for (;mask>1;)
2679 : {
2680 : GEN u, fa, fr;
2681 483 : long n2 = n, rt;
2682 483 : n<<=1; if (mask & 1) n--;
2683 483 : mask >>= 1;
2684 483 : fr = RgXn_red_shallow(f, n);
2685 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2686 483 : an = RgXn_powers(a, rt, n);
2687 483 : if (n>1)
2688 : {
2689 483 : long n4 = (n2+1)>>1;
2690 483 : GEN dfr = RgXn_red_shallow(df, n2);
2691 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2692 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2693 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2694 : }
2695 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2696 483 : fa = RgX_shift(fa, -n2);
2697 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2698 483 : if (gc_needed(av2,2))
2699 : {
2700 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2701 0 : gerepileall(av2, 2, &a, &W);
2702 : }
2703 : }
2704 133 : return gerepileupto(av, a);
2705 : }
2706 :
2707 : GEN
2708 7 : RgXn_sqrt(GEN h, long e)
2709 : {
2710 7 : pari_sp av = avma, av2;
2711 7 : long v = varn(h), n = 1;
2712 7 : GEN f = scalarpol(gen_1, v), df = f;
2713 7 : ulong mask = quadratic_prec_mask(e);
2714 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2715 0 : pari_err_SQRTN("RgXn_sqrt",h);
2716 7 : av2 = avma;
2717 : while(1)
2718 7 : {
2719 14 : long n2 = n, m;
2720 : GEN g;
2721 14 : n<<=1; if (mask & 1) n--;
2722 14 : mask >>= 1;
2723 14 : m = n-n2;
2724 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2725 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2726 14 : if (mask==1) return gerepileupto(av, f);
2727 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2728 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2729 7 : if (gc_needed(av2,2))
2730 : {
2731 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2732 0 : gerepileall(av2, 2, &f, &df);
2733 : }
2734 : }
2735 : }
2736 :
2737 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2738 : GEN
2739 116675 : RgXQ_powu(GEN x, ulong n, GEN T)
2740 : {
2741 116675 : pari_sp av = avma;
2742 :
2743 116675 : if (!n) return pol_1(varn(x));
2744 68442 : if (n == 1) return RgX_copy(x);
2745 19824 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2746 19824 : return gc_GEN(av, x);
2747 : }
2748 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2749 : GEN
2750 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2751 : {
2752 : pari_sp av;
2753 102160 : long s = signe(n);
2754 :
2755 102160 : if (!s) return pol_1(varn(x));
2756 102160 : if (is_pm1(n) == 1)
2757 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2758 13853 : av = avma;
2759 13853 : if (s < 0) x = RgXQ_inv(x, T);
2760 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2761 13853 : return gc_GEN(av, x);
2762 : }
2763 : static GEN
2764 200566 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2765 : static GEN
2766 111152 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2767 :
2768 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2769 : GEN
2770 11746 : ZXQ_powers(GEN x, long l, GEN T)
2771 : {
2772 11746 : int use_sqr = 2*degpol(x) >= degpol(T);
2773 11746 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2774 : }
2775 :
2776 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2777 : GEN
2778 168506 : ZXQ_powu(GEN x, ulong n, GEN T)
2779 : {
2780 168506 : pari_sp av = avma;
2781 :
2782 168506 : if (!n) return pol_1(varn(x));
2783 168506 : if (n == 1) return ZX_copy(x);
2784 113872 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2785 113879 : return gc_GEN(av, x);
2786 : }
2787 :
2788 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2789 : GEN
2790 8554 : RgXQ_powers(GEN x, long l, GEN T)
2791 : {
2792 8554 : int use_sqr = 2*degpol(x) >= degpol(T);
2793 8554 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2794 : }
2795 :
2796 : GEN
2797 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2798 : {
2799 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2800 : }
2801 :
2802 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2803 : GEN
2804 9191 : QXQ_powers(GEN a, long n, GEN T)
2805 : {
2806 : GEN den, v;
2807 9191 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2808 9163 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2809 : /* den*a integral; v[i+1] = (den*a)^i in K */
2810 9163 : if (den)
2811 : { /* restore denominators */
2812 5620 : GEN d = den;
2813 : long i;
2814 5620 : gel(v,2) = a;
2815 27280 : for (i=3; i<=n+1; i++) {
2816 21660 : d = mulii(d,den);
2817 21660 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2818 : }
2819 : }
2820 9163 : return v;
2821 : }
2822 :
2823 : static GEN
2824 3241 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2825 : {
2826 3241 : long l, i, m = 0;
2827 : GEN dz, z;
2828 3241 : GEN V = cgetg_copy(v, &l);
2829 10794 : for (i = imin; i < l; i++)
2830 : {
2831 7553 : GEN c = gel(v, i);
2832 7553 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2833 : }
2834 3241 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2835 3479 : for (i = 1; i < imin; i++) V[i] = v[i];
2836 10794 : for (i = imin; i < l; i++)
2837 : {
2838 7553 : GEN c = gel(v,i);
2839 7553 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2840 7553 : gel(V,i) = c;
2841 : }
2842 3241 : return V;
2843 : }
2844 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2845 : GEN
2846 3003 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2847 3003 : { return do_QXQ_eval(v, 1, a, T); }
2848 :
2849 : GEN
2850 238 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2851 238 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2852 :
2853 : GEN
2854 2940 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2855 : {
2856 2940 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2857 : }
2858 :
2859 : GEN
2860 5144 : RgXQ_norm(GEN x, GEN T)
2861 : {
2862 : pari_sp av;
2863 5144 : long dx = degpol(x);
2864 : GEN L, y;
2865 5144 : if (degpol(T)==0) return gpowgs(x,0);
2866 5137 : av = avma; y = resultant(T, x);
2867 5137 : L = leading_coeff(T);
2868 5137 : if (gequal1(L) || !signe(x)) return y;
2869 0 : return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2870 : }
2871 :
2872 : GEN
2873 476 : RgXQ_trace(GEN x, GEN T)
2874 : {
2875 476 : pari_sp av = avma;
2876 : GEN dT, z;
2877 : long n;
2878 476 : if (degpol(T)==0) return gmulgs(x,0);
2879 469 : dT = RgX_deriv(T); n = degpol(dT);
2880 469 : z = RgXQ_mul(x, dT, T);
2881 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2882 420 : return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2883 : }
2884 :
2885 : GEN
2886 2788358 : RgX_blocks(GEN P, long n, long m)
2887 : {
2888 2788358 : GEN z = cgetg(m+1,t_VEC);
2889 2788358 : long i,j, k=2, l = lg(P);
2890 8558652 : for(i=1; i<=m; i++)
2891 : {
2892 5770294 : GEN zi = cgetg(n+2,t_POL);
2893 5770294 : zi[1] = P[1];
2894 5770294 : gel(z,i) = zi;
2895 17014348 : for(j=2; j<n+2; j++)
2896 11244054 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2897 5770294 : zi = RgX_renormalize_lg(zi, n+2);
2898 : }
2899 2788358 : return z;
2900 : }
2901 :
2902 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2903 : void
2904 49732077 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2905 : {
2906 49732077 : long n = degpol(p), v = varn(p), n0, n1, i;
2907 : GEN p0, p1;
2908 :
2909 49733029 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2910 :
2911 49639752 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2912 49639752 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2913 49642119 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2914 146714448 : for (i=0; i<n1; i++)
2915 : {
2916 97071251 : p0[2+i] = p[2+(i<<1)];
2917 97071251 : p1[2+i] = p[3+(i<<1)];
2918 : }
2919 49643197 : if (n1 != n0)
2920 18008197 : p0[2+i] = p[2+(i<<1)];
2921 49643197 : *pe = normalizepol(p0);
2922 49646687 : *po = normalizepol(p1);
2923 : }
2924 :
2925 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2926 : GEN
2927 43631 : RgX_splitting(GEN p, long k)
2928 : {
2929 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
2930 : GEN r;
2931 :
2932 43631 : m = n/k;
2933 43631 : r = cgetg(k+1,t_VEC);
2934 234059 : for(i=1; i<=k; i++)
2935 : {
2936 190428 : gel(r,i) = cgetg(m+3, t_POL);
2937 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2938 : }
2939 571991 : for (j=1, i=0, l=2; i<=n; i++)
2940 : {
2941 528360 : gmael(r,j,l) = gel(p,2+i);
2942 528360 : if (j==k) { j=1; l++; } else j++;
2943 : }
2944 234059 : for(i=1; i<=k; i++)
2945 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2946 43631 : return r;
2947 : }
2948 :
2949 : /*******************************************************************/
2950 : /* */
2951 : /* Kronecker form */
2952 : /* */
2953 : /*******************************************************************/
2954 :
2955 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2956 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2957 : * normalized coefficients */
2958 : GEN
2959 186168 : Kronecker_to_mod(GEN z, GEN T)
2960 : {
2961 186168 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2962 186168 : GEN x, t = cgetg(N,t_POL);
2963 186168 : t[1] = T[1];
2964 186168 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2965 186168 : x[1] = z[1];
2966 186168 : T = RgX_copy(T);
2967 1399395 : for (i=2; i<lx+2; i++, z+= N-2)
2968 : {
2969 5761180 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2970 1213227 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2971 : }
2972 186168 : N = (l-2) % (N-2) + 2;
2973 472160 : for (j=2; j<N; j++) t[j] = z[j];
2974 186168 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2975 186168 : return normalizepol_lg(x, i+1);
2976 : }
2977 :
2978 : /*******************************************************************/
2979 : /* */
2980 : /* Domain detection */
2981 : /* */
2982 : /*******************************************************************/
2983 :
2984 : static GEN
2985 609812 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2986 : {
2987 609812 : pari_sp av = avma;
2988 : GEN r;
2989 609812 : if (lgefint(p) == 3)
2990 : {
2991 558456 : ulong pp = uel(p, 2);
2992 558456 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2993 : RgX_to_Flx(y, pp), pp));
2994 : }
2995 : else
2996 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2997 609812 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2998 544614 : return gerepileupto(av, FpX_to_mod(r, p));
2999 : }
3000 :
3001 : static GEN
3002 616 : zero_FpXQX_mod(GEN pol, GEN p, long v)
3003 : {
3004 616 : GEN r = cgetg(3,t_POL);
3005 616 : r[1] = evalvarn(v);
3006 616 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
3007 616 : return r;
3008 : }
3009 :
3010 : static GEN
3011 2079 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3012 : {
3013 2079 : pari_sp av = avma;
3014 : long dT;
3015 : GEN kx, ky, r;
3016 2079 : GEN T = RgX_to_FpX(pol, p);
3017 2079 : if (signe(T)==0) pari_err_OP("*", x, y);
3018 2072 : dT = degpol(T);
3019 2072 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3020 2072 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3021 2072 : r = FpX_mul(kx, ky, p);
3022 2072 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3023 1463 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3024 : }
3025 :
3026 : static GEN
3027 455141 : RgX_liftred(GEN x, GEN T)
3028 455141 : { return RgXQX_red(liftpol_shallow(x), T); }
3029 :
3030 : static GEN
3031 171091 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3032 : {
3033 171091 : pari_sp av = avma;
3034 171091 : long dT = degpol(T);
3035 171091 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3036 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3037 171091 : return gerepileupto(av, Kronecker_to_mod(r, T));
3038 : }
3039 :
3040 : static GEN
3041 1407 : RgX_sqr_FpX(GEN x, GEN p)
3042 : {
3043 1407 : pari_sp av = avma;
3044 : GEN r;
3045 1407 : if (lgefint(p) == 3)
3046 : {
3047 1203 : ulong pp = uel(p, 2);
3048 1203 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3049 : }
3050 : else
3051 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3052 1407 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3053 1274 : return gerepileupto(av, FpX_to_mod(r, p));
3054 : }
3055 :
3056 : static GEN
3057 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3058 : {
3059 196 : pari_sp av = avma;
3060 : long dT;
3061 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3062 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3063 189 : dT = degpol(T);
3064 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3065 189 : r = FpX_sqr(kx, p);
3066 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3067 189 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3068 : }
3069 :
3070 : static GEN
3071 13425 : RgX_sqr_QXQX(GEN x, GEN T)
3072 : {
3073 13425 : pari_sp av = avma;
3074 13425 : long dT = degpol(T);
3075 13425 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3076 13425 : return gerepileupto(av, Kronecker_to_mod(r, T));
3077 : }
3078 :
3079 : static GEN
3080 68320 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3081 : {
3082 68320 : pari_sp av = avma;
3083 : GEN r;
3084 68320 : if (lgefint(p) == 3)
3085 : {
3086 51684 : ulong pp = uel(p, 2);
3087 51684 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3088 : RgX_to_Flx(y, pp), pp));
3089 : }
3090 : else
3091 16636 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3092 68320 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3093 31024 : return gerepileupto(av, FpX_to_mod(r, p));
3094 : }
3095 :
3096 : static GEN
3097 49767 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3098 : {
3099 49767 : pari_sp av = avma;
3100 : GEN r;
3101 49767 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3102 49767 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3103 : }
3104 : static GEN
3105 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3106 : {
3107 70 : pari_sp av = avma;
3108 : GEN r;
3109 70 : GEN T = RgX_to_FpX(pol, p);
3110 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3111 63 : if (lgefint(p) == 3)
3112 : {
3113 55 : ulong pp = uel(p, 2);
3114 55 : GEN Tp = ZX_to_Flx(T, pp);
3115 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3116 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3117 : }
3118 : else
3119 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3120 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3121 56 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3122 : }
3123 :
3124 : static GEN
3125 121130765 : RgX_mul_fast(GEN x, GEN y)
3126 : {
3127 : GEN p, pol;
3128 : long pa;
3129 121130765 : long t = RgX_type2(x,y, &p,&pol,&pa);
3130 121131101 : switch(t)
3131 : {
3132 67189445 : case t_INT: return ZX_mul(x,y);
3133 3356236 : case t_FRAC: return QX_mul(x,y);
3134 105127 : case t_FFELT: return FFX_mul(x, y, pol);
3135 609812 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3136 171091 : case RgX_type_code(t_POLMOD, t_INT):
3137 : case RgX_type_code(t_POLMOD, t_FRAC):
3138 171091 : return RgX_mul_QXQX(x, y, pol);
3139 2058 : case RgX_type_code(t_POLMOD, t_INTMOD):
3140 2058 : return RgX_mul_FpXQX(x, y, pol, p);
3141 49697332 : default: return NULL;
3142 : }
3143 : }
3144 : static GEN
3145 1384641 : RgX_sqr_fast(GEN x)
3146 : {
3147 : GEN p, pol;
3148 : long pa;
3149 1384641 : long t = RgX_type(x,&p,&pol,&pa);
3150 1384647 : switch(t)
3151 : {
3152 1245852 : case t_INT: return ZX_sqr(x);
3153 83309 : case t_FRAC: return QX_sqr(x);
3154 2499 : case t_FFELT: return FFX_sqr(x, pol);
3155 1407 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3156 13425 : case RgX_type_code(t_POLMOD, t_INT):
3157 : case RgX_type_code(t_POLMOD, t_FRAC):
3158 13425 : return RgX_sqr_QXQX(x, pol);
3159 194 : case RgX_type_code(t_POLMOD, t_INTMOD):
3160 194 : return RgX_sqr_FpXQX(x, pol, p);
3161 37961 : default: return NULL;
3162 : }
3163 : }
3164 :
3165 : static GEN
3166 14689295 : RgX_rem_fast(GEN x, GEN y)
3167 : {
3168 : GEN p, pol;
3169 : long pa;
3170 14689295 : long t = RgX_type2(x,y, &p,&pol,&pa);
3171 14689462 : switch(t)
3172 : {
3173 3471718 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3174 1835456 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3175 84 : case t_FFELT: return FFX_rem(x, y, pol);
3176 68320 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3177 49767 : case RgX_type_code(t_POLMOD, t_INT):
3178 : case RgX_type_code(t_POLMOD, t_FRAC):
3179 49767 : return RgX_rem_QXQX(x, y, pol);
3180 63 : case RgX_type_code(t_POLMOD, t_INTMOD):
3181 63 : return RgX_rem_FpXQX(x, y, pol, p);
3182 9264054 : default: return NULL;
3183 : }
3184 : }
3185 :
3186 : GEN
3187 107398079 : RgX_mul(GEN x, GEN y)
3188 : {
3189 107398079 : GEN z = RgX_mul_fast(x,y);
3190 107398242 : if (!z) z = RgX_mul_i(x,y);
3191 107397808 : return z;
3192 : }
3193 :
3194 : GEN
3195 1372489 : RgX_sqr(GEN x)
3196 : {
3197 1372489 : GEN z = RgX_sqr_fast(x);
3198 1372494 : if (!z) z = RgX_sqr_i(x);
3199 1372494 : return z;
3200 : }
3201 :
3202 : GEN
3203 14689302 : RgX_rem(GEN x, GEN y)
3204 : {
3205 14689302 : GEN z = RgX_rem_fast(x, y);
3206 14689447 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3207 14689355 : return z;
3208 : }
3209 :
3210 : static GEN
3211 0 : _RgX_mul(void* E, GEN x, GEN y)
3212 0 : { (void) E; return RgX_mul(x, y); }
3213 :
3214 : GEN
3215 0 : RgXV_prod(GEN V)
3216 0 : { return gen_product(V, NULL, &_RgX_mul); }
3217 :
3218 : GEN
3219 11061650 : RgXn_mul(GEN f, GEN g, long n)
3220 : {
3221 11061650 : pari_sp av = avma;
3222 11061650 : GEN h = RgX_mul_fast(f,g);
3223 11061644 : if (!h) return RgXn_mul2(f,g,n);
3224 1783450 : if (degpol(h) < n) return h;
3225 382975 : return gc_GEN(av, RgXn_red_shallow(h, n));
3226 : }
3227 :
3228 : GEN
3229 12152 : RgXn_sqr(GEN f, long n)
3230 : {
3231 12152 : pari_sp av = avma;
3232 12152 : GEN g = RgX_sqr_fast(f);
3233 12152 : if (!g) return RgXn_sqr2(f,n);
3234 11837 : if (degpol(g) < n) return g;
3235 10640 : return gc_GEN(av, RgXn_red_shallow(g, n));
3236 : }
3237 :
3238 : /* (f*g) \/ x^n */
3239 : GEN
3240 2671060 : RgX_mulhigh_i(GEN f, GEN g, long n)
3241 : {
3242 2671060 : GEN h = RgX_mul_fast(f,g);
3243 2671062 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3244 : }
3245 :
3246 : /* (f*g) \/ x^n */
3247 : GEN
3248 0 : RgX_sqrhigh_i(GEN f, long n)
3249 : {
3250 0 : GEN h = RgX_sqr_fast(f);
3251 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3252 : }
|
