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 19446966 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19446966 : long pold=1, rold=n*(d-1);
35 90973453 : for(p=2; p<=d; p++)
36 : {
37 71526487 : r = m*(p-1) + n*((d-1)/p);
38 71526487 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19446966 : return pold;
41 : }
42 :
43 : static GEN
44 14494218 : 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 14494218 : pari_sp av = avma;
48 : long i;
49 14494218 : GEN z = cmul(E,P,a,ff->one(E));
50 14469613 : if (!z) z = gen_0;
51 75590478 : for (i=1; i<=n; i++)
52 : {
53 61120650 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 61144832 : if (t) {
55 46196621 : z = ff->add(E, z, t);
56 46164385 : if (gc_needed(av,2)) z = gerepileupto(av, z);
57 : }
58 : }
59 14469828 : 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 12445155 : 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 12445155 : pari_sp av = avma;
72 12445155 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12445155 : if (d < 0) return ff->zero(E);
76 11838480 : if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1148512 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1148512 : 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 1148512 : d -= l;
84 1148512 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2654728 : while (d >= l-1)
86 : {
87 1505008 : d -= l-1;
88 1505008 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1505064 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1504990 : if (gc_needed(av,2))
91 91 : z = gerepileupto(av, z);
92 : }
93 1149720 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1149721 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1149726 : return gerepileupto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 868903 : 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 868903 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 868903 : if (d < 0) return ff->zero(E);
106 867685 : rtd = (long) sqrt((double)d);
107 867685 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 867693 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 867694 : return gerepileupto(av, z);
110 : }
111 :
112 : static GEN
113 2027929 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19244729 : _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 1840019 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 591655 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2068893 : _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 667289 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 667289 : pari_sp av = avma, btop;
163 667289 : long i, d = degpol(P), o;
164 : GEN s;
165 667288 : if (signe(P)==0) return pol_0(varn(P));
166 667288 : s = gel(P, d+2);
167 667288 : if (d == 0) return gcopy(s);
168 664012 : o = RgX_deflate_order(P);
169 664022 : if (o > 1) A = gpowgs(A, o);
170 664022 : btop = avma;
171 2310783 : for (i = d-o; i >= 0; i-=o)
172 : {
173 1646762 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 1646761 : 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 664021 : 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 281046 : get_Rg_algebra(void)
207 : {
208 281046 : 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 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
225 : }
226 :
227 : /*******************************************************************/
228 : /* */
229 : /* RgX */
230 : /* */
231 : /*******************************************************************/
232 :
233 : long
234 24804298 : RgX_equal(GEN x, GEN y)
235 : {
236 24804298 : long i = lg(x);
237 :
238 24804298 : if (i != lg(y)) return 0;
239 106783915 : for (i--; i > 1; i--)
240 82300858 : if (!gequal(gel(x,i),gel(y,i))) return 0;
241 24483057 : return 1;
242 : }
243 :
244 : /* Returns 1 in the base ring over which x is defined */
245 : /* HACK: this also works for t_SER */
246 : GEN
247 156370317 : Rg_get_1(GEN x)
248 : {
249 : GEN p, T;
250 156370317 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
251 156370318 : if (RgX_type_is_composite(tx))
252 11634460 : RgX_type_decode(tx, &i /*junk*/, &tx);
253 156370318 : switch(tx)
254 : {
255 427217 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
256 945 : case t_PADIC: return cvtop(gen_1, p, lx);
257 4998 : case t_FFELT: return FF_1(T);
258 155937158 : default: return gen_1;
259 : }
260 : }
261 : /* Returns 0 in the base ring over which x is defined */
262 : /* HACK: this also works for t_SER */
263 : GEN
264 6979168 : Rg_get_0(GEN x)
265 : {
266 : GEN p, T;
267 6979168 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
268 6979168 : if (RgX_type_is_composite(tx))
269 56504 : RgX_type_decode(tx, &i /*junk*/, &tx);
270 6979168 : switch(tx)
271 : {
272 497 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
273 42 : case t_PADIC: return zeropadic(p, lx);
274 231 : case t_FFELT: return FF_zero(T);
275 6978398 : default: return gen_0;
276 : }
277 : }
278 :
279 : GEN
280 6489 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
281 : {
282 6489 : long i, n = degpol(P);
283 : GEN z, dz, dP;
284 6489 : if (n < 0) return gen_0;
285 6489 : P = Q_remove_denom(P, &dP);
286 6489 : z = gel(P,2); if (n == 0) return icopy(z);
287 3717 : if (dV) z = mulii(dV, z); /* V[1] = dV */
288 3717 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
289 7308 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
290 3717 : dz = mul_denom(dP, dV);
291 3717 : return dz? RgX_Rg_div(z, dz): z;
292 : }
293 :
294 : /* Return P(h * x), not memory clean */
295 : GEN
296 48261 : RgX_unscale(GEN P, GEN h)
297 : {
298 48261 : long i, l = lg(P);
299 48261 : GEN hi = gen_1, Q = cgetg(l, t_POL);
300 48261 : Q[1] = P[1];
301 48261 : if (l == 2) return Q;
302 38503 : gel(Q,2) = gcopy(gel(P,2));
303 91686 : for (i=3; i<l; i++)
304 : {
305 53182 : hi = gmul(hi,h);
306 53181 : gel(Q,i) = gmul(gel(P,i), hi);
307 : }
308 38504 : return Q;
309 : }
310 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
311 : GEN
312 1473125 : ZX_z_unscale(GEN P, long h)
313 : {
314 1473125 : long i, l = lg(P);
315 1473125 : GEN Q = cgetg(l, t_POL);
316 1473125 : Q[1] = P[1];
317 1473125 : if (l == 2) return Q;
318 1384579 : gel(Q,2) = gel(P,2);
319 1384579 : if (l == 3) return Q;
320 1358273 : if (h == -1)
321 236837 : for (i = 3; i < l; i++)
322 : {
323 196286 : gel(Q,i) = negi(gel(P,i));
324 196286 : if (++i == l) break;
325 145324 : gel(Q,i) = gel(P,i);
326 : }
327 : else
328 : {
329 : GEN hi;
330 1266760 : gel(Q,3) = mulis(gel(P,3), h);
331 1266760 : hi = sqrs(h);
332 5450358 : for (i = 4; i < l; i++)
333 : {
334 4183596 : gel(Q,i) = mulii(gel(P,i), hi);
335 4183601 : if (i != l-1) hi = mulis(hi,h);
336 : }
337 : }
338 1358275 : return Q;
339 : }
340 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
341 : GEN
342 1019154 : ZX_unscale(GEN P, GEN h)
343 : {
344 : long i, l;
345 : GEN Q, hi;
346 1019154 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
347 881 : l = lg(P); Q = cgetg(l, t_POL);
348 881 : Q[1] = P[1];
349 881 : if (l == 2) return Q;
350 881 : gel(Q,2) = gel(P,2);
351 881 : if (l == 3) return Q;
352 881 : hi = h;
353 881 : gel(Q,3) = mulii(gel(P,3), hi);
354 2741 : for (i = 4; i < l; i++)
355 : {
356 1860 : hi = mulii(hi,h);
357 1860 : gel(Q,i) = mulii(gel(P,i), hi);
358 : }
359 881 : return Q;
360 : }
361 : /* P a ZX. Return P(x << n), not memory clean */
362 : GEN
363 928273 : ZX_unscale2n(GEN P, long n)
364 : {
365 928273 : long i, ni = n, l = lg(P);
366 928273 : GEN Q = cgetg(l, t_POL);
367 928273 : Q[1] = P[1];
368 928273 : if (l == 2) return Q;
369 928273 : gel(Q,2) = gel(P,2);
370 928273 : if (l == 3) return Q;
371 928273 : gel(Q,3) = shifti(gel(P,3), ni);
372 3172380 : for (i=4; i<l; i++)
373 : {
374 2244149 : ni += n;
375 2244149 : gel(Q,i) = shifti(gel(P,i), ni);
376 : }
377 928231 : return Q;
378 : }
379 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
380 : GEN
381 12997 : ZX_unscale_div(GEN P, GEN h)
382 : {
383 12997 : long i, l = lg(P);
384 12997 : GEN hi, Q = cgetg(l, t_POL);
385 12997 : Q[1] = P[1];
386 12997 : if (l == 2) return Q;
387 12997 : gel(Q,2) = diviiexact(gel(P,2), h);
388 12997 : if (l == 3) return Q;
389 12997 : gel(Q,3) = gel(P,3);
390 12997 : if (l == 4) return Q;
391 12997 : hi = h;
392 12997 : gel(Q,4) = mulii(gel(P,4), hi);
393 64121 : for (i=5; i<l; i++)
394 : {
395 51124 : hi = mulii(hi,h);
396 51124 : gel(Q,i) = mulii(gel(P,i), hi);
397 : }
398 12997 : return Q;
399 : }
400 : /* P(h*X) / h^k, assuming the result is a ZX */
401 : GEN
402 1393 : ZX_unscale_divpow(GEN P, GEN h, long k)
403 : {
404 1393 : long i, j, l = lg(P);
405 1393 : GEN H, Q = cgetg(l, t_POL);
406 1393 : Q[1] = P[1]; if (l == 2) return Q;
407 1393 : H = gpowers(h, maxss(k, l - 3 - k));
408 5572 : for (i = 2, j = k+1; j > 1 && i < l; i++)
409 4179 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
410 1393 : if (i == l) return Q;
411 1393 : gel(Q, i) = gel(P, i); i++;
412 5082 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
413 1393 : return Q;
414 : }
415 :
416 : GEN
417 6237 : RgXV_unscale(GEN x, GEN h)
418 : {
419 6237 : if (isint1(h)) return gcopy(x);
420 18253 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
421 : }
422 :
423 : /* Return h^degpol(P) P(x / h), not memory clean */
424 : GEN
425 4340655 : RgX_rescale(GEN P, GEN h)
426 : {
427 4340655 : long i, l = lg(P);
428 4340655 : GEN Q = cgetg(l,t_POL), hi = h;
429 4340654 : gel(Q,l-1) = gel(P,l-1);
430 11420200 : for (i=l-2; i>=2; i--)
431 : {
432 11417567 : gel(Q,i) = gmul(gel(P,i), hi);
433 11417440 : if (i == 2) break;
434 7079264 : hi = gmul(hi,h);
435 : }
436 4340809 : Q[1] = P[1]; return Q;
437 : }
438 :
439 : GEN
440 2401 : RgXV_rescale(GEN x, GEN h)
441 : {
442 2401 : if (isint1(h)) return RgX_copy(x);
443 16086 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
444 : }
445 :
446 : /* A(X^d) --> A(X) */
447 : GEN
448 985976 : RgX_deflate(GEN x0, long d)
449 : {
450 : GEN z, y, x;
451 985976 : long i,id, dy, dx = degpol(x0);
452 985972 : if (d == 1 || dx <= 0) return leafcopy(x0);
453 375531 : dy = dx/d;
454 375531 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
455 375531 : z = y + 2;
456 375531 : x = x0+ 2;
457 1511403 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
458 375531 : return y;
459 : }
460 :
461 : GEN
462 1260 : RgX_homogenize_deg(GEN P, long d, long v)
463 : {
464 : long i, l;
465 1260 : GEN Q = cgetg_copy(P, &l);
466 1260 : Q[1] = P[1];
467 4018 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
468 1260 : return Q;
469 : }
470 :
471 : GEN
472 18732 : RgX_homogenize(GEN P, long v)
473 : {
474 : long i, l, d;
475 18732 : GEN Q = cgetg_copy(P, &l);
476 18732 : Q[1] = P[1]; d = l-3;
477 165998 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
478 18732 : return Q;
479 : }
480 :
481 : /* F a t_RFRAC */
482 : long
483 105 : rfrac_deflate_order(GEN F)
484 : {
485 105 : GEN N = gel(F,1), D = gel(F,2);
486 105 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
487 105 : if (m == 1) return 1;
488 42 : if (typ(N) == t_POL && varn(N) == varn(D))
489 28 : m = cgcd(m, RgX_deflate_order(N));
490 42 : return m;
491 : }
492 : /* F a t_RFRAC */
493 : GEN
494 105 : rfrac_deflate_max(GEN F, long *m)
495 : {
496 105 : *m = rfrac_deflate_order(F);
497 105 : return rfrac_deflate(F, *m);
498 : }
499 : /* F a t_RFRAC */
500 : GEN
501 105 : rfrac_deflate(GEN F, long m)
502 : {
503 105 : GEN N = gel(F,1), D = gel(F,2);
504 105 : if (m == 1) return F;
505 42 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
506 42 : D = RgX_deflate(D, m); return mkrfrac(N, D);
507 : }
508 :
509 : /* return x0(X^d) */
510 : GEN
511 798657 : RgX_inflate(GEN x0, long d)
512 : {
513 798657 : long i, id, dy, dx = degpol(x0);
514 798657 : GEN x = x0 + 2, z, y;
515 798657 : if (dx <= 0) return leafcopy(x0);
516 791041 : dy = dx*d;
517 791041 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
518 791041 : z = y + 2;
519 27878129 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
520 11422007 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
521 791041 : return y;
522 : }
523 :
524 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
525 : static GEN
526 5697125 : RgX_translate_basecase(GEN P, GEN c)
527 : {
528 5697125 : pari_sp av = avma;
529 : GEN Q, R;
530 : long i, k, n;
531 :
532 5697125 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
533 5697104 : Q = leafcopy(P);
534 5697105 : R = Q+2; n = degpol(P);
535 5697103 : if (isint1(c))
536 : {
537 7595 : for (i=1; i<=n; i++)
538 : {
539 20265 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
540 5390 : if (gc_needed(av,2))
541 : {
542 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
543 0 : Q = gerepilecopy(av, Q); R = Q+2;
544 : }
545 : }
546 : }
547 5694897 : else if (isintm1(c))
548 : {
549 15834 : for (i=1; i<=n; i++)
550 : {
551 49532 : for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
552 12250 : if (gc_needed(av,2))
553 : {
554 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
555 0 : Q = gerepilecopy(av, Q); R = Q+2;
556 : }
557 : }
558 : }
559 : else
560 : {
561 19413935 : for (i=1; i<=n; i++)
562 : {
563 45414616 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
564 13722620 : if (gc_needed(av,2))
565 : {
566 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
567 0 : Q = gerepilecopy(av, Q); R = Q+2;
568 : }
569 : }
570 : }
571 5696780 : return gerepilecopy(av, Q);
572 : }
573 :
574 : static GEN
575 142269 : zero_FpX_mod(GEN p, long v)
576 : {
577 142269 : GEN r = cgetg(3,t_POL);
578 142269 : r[1] = evalvarn(v);
579 142269 : gel(r,2) = mkintmod(gen_0, icopy(p));
580 142269 : return r;
581 : }
582 :
583 : static GEN
584 0 : RgX_translate_FpX(GEN P, GEN c, GEN p)
585 : {
586 0 : pari_sp av = avma;
587 : GEN r;
588 : #if 0
589 : if (lgefint(p) == 3)
590 : {
591 : ulong pp = uel(p, 2);
592 : r = Flx_to_ZX_inplace(Flx_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
593 : }
594 : else
595 : #endif
596 0 : r = FpX_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
597 0 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
598 0 : return gerepileupto(av, FpX_to_mod(r, p));
599 : }
600 :
601 : static GEN
602 5978591 : RgX_translate_fast(GEN P, GEN c)
603 : {
604 : GEN p, pol;
605 : long pa;
606 5978591 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
607 5978590 : switch(t)
608 : {
609 281648 : case t_INT: return ZX_translate(P, c);
610 0 : case t_INTMOD: return RgX_translate_FpX(P, c, p);
611 5696943 : default: return NULL;
612 : }
613 : }
614 :
615 : static GEN
616 5697319 : RgX_translate_i(GEN P, GEN c)
617 : {
618 5697319 : pari_sp av = avma;
619 : long n;
620 5697319 : n = degpol(P);
621 5697316 : if (n < 40)
622 5697127 : return RgX_translate_basecase(P, c);
623 : else
624 : {
625 189 : long d = n >> 1;
626 189 : GEN Q = RgX_translate_i(RgX_shift_shallow(P, -d), c);
627 189 : GEN R = RgX_translate_i(RgXn_red_shallow(P, d), c);
628 189 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
629 189 : return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
630 : }
631 : }
632 :
633 : GEN
634 5978591 : RgX_translate(GEN P, GEN c)
635 : {
636 5978591 : GEN R = RgX_translate_fast(P, c);
637 5978590 : return R ? R: RgX_translate_i(P,c);
638 : }
639 : /* P(ax + b) */
640 : GEN
641 27062 : RgX_affine(GEN P, GEN a, GEN b)
642 : {
643 27062 : if (!gequal0(b)) P = RgX_translate(P, b);
644 27062 : return RgX_unscale(P, a);
645 : }
646 :
647 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
648 : GEN
649 29685 : RgXQX_translate(GEN P, GEN c, GEN T)
650 : {
651 29685 : pari_sp av = avma;
652 : GEN Q, R;
653 : long i, k, n;
654 :
655 29685 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
656 29342 : Q = leafcopy(P);
657 29342 : R = Q+2; n = degpol(P);
658 95079 : for (i=1; i<=n; i++)
659 : {
660 279782 : for (k=n-i; k<n; k++)
661 : {
662 214045 : pari_sp av2 = avma;
663 214045 : gel(R,k) = gerepileupto(av2,
664 214045 : RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
665 : }
666 65737 : if (gc_needed(av,2))
667 : {
668 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
669 0 : Q = gerepilecopy(av, Q); R = Q+2;
670 : }
671 : }
672 29342 : return gerepilecopy(av, Q);
673 : }
674 :
675 : /********************************************************************/
676 : /** **/
677 : /** CONVERSIONS **/
678 : /** (not memory clean) **/
679 : /** **/
680 : /********************************************************************/
681 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
682 : * but everything else is */
683 : static GEN
684 167587 : QXQ_to_mod(GEN x, GEN T)
685 : {
686 : long d;
687 167587 : switch(typ(x))
688 : {
689 62696 : case t_INT: return icopy(x);
690 2079 : case t_FRAC: return gcopy(x);
691 102812 : case t_POL:
692 102812 : d = degpol(x);
693 102812 : if (d < 0) return gen_0;
694 100544 : if (d == 0) return gcopy(gel(x,2));
695 98398 : return mkpolmod(RgX_copy(x), T);
696 0 : default: pari_err_TYPE("QXQ_to_mod",x);
697 : return NULL;/* LCOV_EXCL_LINE */
698 : }
699 : }
700 : /* pure shallow version */
701 : GEN
702 807615 : QXQ_to_mod_shallow(GEN x, GEN T)
703 : {
704 : long d;
705 807615 : switch(typ(x))
706 : {
707 519246 : case t_INT:
708 519246 : case t_FRAC: return x;
709 288369 : case t_POL:
710 288369 : d = degpol(x);
711 288369 : if (d < 0) return gen_0;
712 242116 : if (d == 0) return gel(x,2);
713 225489 : return mkpolmod(x, T);
714 0 : default: pari_err_TYPE("QXQ_to_mod",x);
715 : return NULL;/* LCOV_EXCL_LINE */
716 : }
717 : }
718 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
719 : * Not memory clean because T not copied, but everything else is */
720 : static GEN
721 37415 : QXQX_to_mod(GEN z, GEN T)
722 : {
723 37415 : long i,l = lg(z);
724 37415 : GEN x = cgetg(l,t_POL);
725 186753 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
726 37415 : x[1] = z[1]; return normalizepol_lg(x,l);
727 : }
728 : /* pure shallow version */
729 : GEN
730 192308 : QXQX_to_mod_shallow(GEN z, GEN T)
731 : {
732 192308 : long i,l = lg(z);
733 192308 : GEN x = cgetg(l,t_POL);
734 917099 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
735 192308 : x[1] = z[1]; return normalizepol_lg(x,l);
736 : }
737 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
738 : GEN
739 12614 : QXQXV_to_mod(GEN V, GEN T)
740 : {
741 12614 : long i, l = lg(V);
742 12614 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
743 50029 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
744 12614 : return z;
745 : }
746 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
747 : GEN
748 18267 : QXQV_to_mod(GEN V, GEN T)
749 : {
750 18267 : long i, l = lg(V);
751 18267 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
752 36516 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
753 18267 : return z;
754 : }
755 :
756 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
757 : GEN
758 14448 : QXQC_to_mod_shallow(GEN V, GEN T)
759 : {
760 14448 : long i, l = lg(V);
761 14448 : GEN z = cgetg(l, t_COL);
762 97272 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
763 14448 : return z;
764 : }
765 :
766 : GEN
767 6587 : QXQM_to_mod_shallow(GEN V, GEN T)
768 : {
769 6587 : long i, l = lg(V);
770 6587 : GEN z = cgetg(l, t_MAT);
771 21035 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
772 6587 : return z;
773 : }
774 :
775 : GEN
776 6166084 : RgX_renormalize_lg(GEN x, long lx)
777 : {
778 : long i;
779 9190416 : for (i = lx-1; i>1; i--)
780 8738204 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
781 6166084 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
782 6166084 : setlg(x, i+1); setsigne(x, i != 1); return x;
783 : }
784 :
785 : GEN
786 1491527 : RgV_to_RgX(GEN x, long v)
787 : {
788 1491527 : long i, k = lg(x);
789 : GEN p;
790 :
791 4158281 : while (--k && gequal0(gel(x,k)));
792 1491525 : if (!k) return pol_0(v);
793 1484247 : i = k+2; p = cgetg(i,t_POL);
794 1484248 : p[1] = evalsigne(1) | evalvarn(v);
795 10235166 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
796 1484248 : return p;
797 : }
798 : GEN
799 188283 : RgV_to_RgX_reverse(GEN x, long v)
800 : {
801 188283 : long j, k, l = lg(x);
802 : GEN p;
803 :
804 189851 : for (k = 1; k < l; k++)
805 189851 : if (!gequal0(gel(x,k))) break;
806 188283 : if (k == l) return pol_0(v);
807 188283 : k -= 1;
808 188283 : l -= k;
809 188283 : x += k;
810 188283 : p = cgetg(l+1,t_POL);
811 188283 : p[1] = evalsigne(1) | evalvarn(v);
812 985038 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
813 188283 : return p;
814 : }
815 :
816 : /* return the (N-dimensional) vector of coeffs of p */
817 : GEN
818 13369469 : RgX_to_RgC(GEN x, long N)
819 : {
820 : long i, l;
821 : GEN z;
822 13369469 : l = lg(x)-1; x++;
823 13369469 : if (l > N+1) l = N+1; /* truncate higher degree terms */
824 13369469 : z = cgetg(N+1,t_COL);
825 80708816 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
826 24440535 : for ( ; i<=N; i++) gel(z,i) = gen_0;
827 13369582 : return z;
828 : }
829 : GEN
830 1353734 : Rg_to_RgC(GEN x, long N)
831 : {
832 1353734 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
833 : }
834 :
835 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
836 : GEN
837 301239 : RgM_to_RgXV(GEN x, long v)
838 1283141 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
839 : GEN
840 7202 : RgM_to_RgXV_reverse(GEN x, long v)
841 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
842 :
843 : /* matrix whose entries are given by the coeffs of the polynomials in
844 : * vector v (considered as degree n-1 polynomials) */
845 : GEN
846 335268 : RgV_to_RgM(GEN x, long n)
847 1684495 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
848 :
849 : GEN
850 78408 : RgXV_to_RgM(GEN x, long n)
851 398766 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
852 :
853 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
854 : GEN
855 23487 : RgM_to_RgXX(GEN x, long v,long w)
856 : {
857 23487 : long j, lx = lg(x);
858 23487 : GEN y = cgetg(lx+1, t_POL);
859 23487 : y[1] = evalsigne(1) | evalvarn(v);
860 23487 : y++;
861 130516 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
862 23487 : return normalizepol_lg(--y, lx+1);
863 : }
864 :
865 : /* matrix whose entries are given by the coeffs of the polynomial v in
866 : * two variables (considered as degree n-1 polynomials) */
867 : GEN
868 322 : RgXX_to_RgM(GEN v, long n)
869 : {
870 322 : long j, N = lg(v)-1;
871 322 : GEN y = cgetg(N, t_MAT);
872 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
873 322 : return y;
874 : }
875 :
876 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
877 : GEN
878 32209 : RgXY_swapspec(GEN x, long n, long w, long nx)
879 : {
880 32209 : long j, ly = n+3;
881 32209 : GEN y = cgetg(ly, t_POL);
882 32209 : y[1] = evalsigne(1);
883 397386 : for (j=2; j<ly; j++)
884 : {
885 : long k;
886 365177 : GEN a = cgetg(nx+2,t_POL);
887 365177 : a[1] = evalsigne(1) | evalvarn(w);
888 2015615 : for (k=0; k<nx; k++)
889 : {
890 1650438 : GEN xk = gel(x,k);
891 1650438 : if (typ(xk)==t_POL && varn(xk)==w)
892 1556063 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
893 : else
894 94375 : gel(a,k+2) = j==2 ? xk: gen_0;
895 : }
896 365177 : gel(y,j) = normalizepol_lg(a, nx+2);
897 : }
898 32209 : return normalizepol_lg(y,ly);
899 : }
900 :
901 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
902 : GEN
903 2051 : RgXY_swap(GEN x, long n, long w)
904 : {
905 2051 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
906 2051 : setvarn(z, varn(x)); return z;
907 : }
908 :
909 : long
910 1387 : RgXY_degreex(GEN b)
911 : {
912 1387 : long deg = 0, i;
913 1387 : if (!signe(b)) return -1;
914 15592 : for (i = 2; i < lg(b); ++i)
915 : {
916 14205 : GEN bi = gel(b, i);
917 14205 : if (typ(bi) == t_POL)
918 13714 : deg = maxss(deg, degpol(bi));
919 : }
920 1387 : return deg;
921 : }
922 :
923 : GEN
924 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
925 :
926 : /* return (x % X^n). Shallow */
927 : GEN
928 7924182 : RgXn_red_shallow(GEN a, long n)
929 : {
930 7924182 : long i, L = n+2, l = lg(a);
931 : GEN b;
932 7924182 : if (L >= l) return a; /* deg(x) < n */
933 5735344 : b = cgetg(L, t_POL); b[1] = a[1];
934 36584803 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
935 5735346 : return normalizepol_lg(b,L);
936 : }
937 :
938 : GEN
939 483 : RgXnV_red_shallow(GEN x, long n)
940 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
941 :
942 : /* return (x * X^n). Shallow */
943 : GEN
944 175443327 : RgX_shift_shallow(GEN a, long n)
945 : {
946 175443327 : long i, l = lg(a);
947 : GEN b;
948 175443327 : if (l == 2 || !n) return a;
949 110453825 : l += n;
950 110453825 : if (n < 0)
951 : {
952 55521237 : if (l <= 2) return pol_0(varn(a));
953 54049695 : b = cgetg(l, t_POL); b[1] = a[1];
954 54050503 : a -= n;
955 165419561 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
956 : } else {
957 54932588 : b = cgetg(l, t_POL); b[1] = a[1];
958 54936849 : a -= n; n += 2;
959 119219368 : for (i=2; i<n; i++) gel(b,i) = gen_0;
960 211460353 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
961 : }
962 108987352 : return b;
963 : }
964 : /* return (x * X^n). */
965 : GEN
966 1570805 : RgX_shift(GEN a, long n)
967 : {
968 1570805 : long i, l = lg(a);
969 : GEN b;
970 1570805 : if (l == 2 || !n) return RgX_copy(a);
971 1569853 : l += n;
972 1569853 : if (n < 0)
973 : {
974 1442 : if (l <= 2) return pol_0(varn(a));
975 1372 : b = cgetg(l, t_POL); b[1] = a[1];
976 1372 : a -= n;
977 9534 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
978 : } else {
979 1568411 : b = cgetg(l, t_POL); b[1] = a[1];
980 1568411 : a -= n; n += 2;
981 3948983 : for (i=2; i<n; i++) gel(b,i) = gen_0;
982 4333042 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
983 : }
984 1569783 : return b;
985 : }
986 :
987 : GEN
988 317037 : RgX_rotate_shallow(GEN P, long k, long p)
989 : {
990 317037 : long i, l = lgpol(P);
991 : GEN r;
992 317037 : if (signe(P)==0)
993 1365 : return pol_0(varn(P));
994 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
995 2100644 : for(i=0; i<p; i++)
996 : {
997 1784972 : long s = 2+(i+k)%p;
998 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
999 : }
1000 315672 : return RgX_renormalize(r);
1001 : }
1002 :
1003 : GEN
1004 2989488 : RgX_mulXn(GEN x, long d)
1005 : {
1006 : pari_sp av;
1007 : GEN z;
1008 : long v;
1009 2989488 : if (d >= 0) return RgX_shift(x, d);
1010 1471950 : d = -d;
1011 1471950 : v = RgX_val(x);
1012 1471950 : if (v >= d) return RgX_shift(x, -d);
1013 1471936 : av = avma;
1014 1471936 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1015 1471936 : return gerepileupto(av, z);
1016 : }
1017 :
1018 : long
1019 588 : RgXV_maxdegree(GEN x)
1020 : {
1021 588 : long d = -1, i, l = lg(x);
1022 4494 : for (i = 1; i < l; i++)
1023 3906 : d = maxss(d, degpol(gel(x,i)));
1024 588 : return d;
1025 : }
1026 :
1027 : long
1028 3649348 : RgX_val(GEN x)
1029 : {
1030 3649348 : long i, lx = lg(x);
1031 3649348 : if (lx == 2) return LONG_MAX;
1032 4580502 : for (i = 2; i < lx; i++)
1033 4580446 : if (!isexactzero(gel(x,i))) break;
1034 3649187 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1035 3649131 : return i - 2;
1036 : }
1037 : long
1038 82392167 : RgX_valrem(GEN x, GEN *Z)
1039 : {
1040 82392167 : long v, i, lx = lg(x);
1041 82392167 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1042 128251906 : for (i = 2; i < lx; i++)
1043 128253682 : if (!isexactzero(gel(x,i))) break;
1044 : /* possible with nonrational zeros */
1045 82392163 : if (i == lx)
1046 : {
1047 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1048 14 : return LONG_MAX;
1049 : }
1050 82392149 : v = i - 2;
1051 82392149 : *Z = RgX_shift_shallow(x, -v);
1052 82405031 : return v;
1053 : }
1054 : long
1055 852175 : RgX_valrem_inexact(GEN x, GEN *Z)
1056 : {
1057 : long v;
1058 852175 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1059 869155 : for (v = 0;; v++)
1060 869155 : if (!gequal0(gel(x,2+v))) break;
1061 852168 : if (Z) *Z = RgX_shift_shallow(x, -v);
1062 852168 : return v;
1063 : }
1064 :
1065 : GEN
1066 67291 : RgXQC_red(GEN x, GEN T)
1067 410823 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1068 :
1069 : GEN
1070 1211 : RgXQV_red(GEN x, GEN T)
1071 27832 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1072 :
1073 : GEN
1074 12929 : RgXQM_red(GEN x, GEN T)
1075 80220 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1076 :
1077 : GEN
1078 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1079 : {
1080 322 : return RgXQM_red(RgM_mul(P, Q), T);
1081 : }
1082 :
1083 : GEN
1084 500299 : RgXQX_red(GEN P, GEN T)
1085 : {
1086 500299 : long i, l = lg(P);
1087 500299 : GEN Q = cgetg(l, t_POL);
1088 500299 : Q[1] = P[1];
1089 2600050 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1090 500299 : return normalizepol_lg(Q, l);
1091 : }
1092 :
1093 : GEN
1094 824471 : RgX_deriv(GEN x)
1095 : {
1096 824471 : long i,lx = lg(x)-1;
1097 : GEN y;
1098 :
1099 824471 : if (lx<3) return pol_0(varn(x));
1100 821384 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1101 3637442 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1102 821376 : y[1] = x[1]; return normalizepol_lg(y,i);
1103 : }
1104 :
1105 : GEN
1106 1376542 : RgX_recipspec_shallow(GEN x, long l, long n)
1107 : {
1108 : long i;
1109 1376542 : GEN z = cgetg(n+2,t_POL);
1110 1376547 : z[1] = 0; z += 2;
1111 36074052 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1112 1599862 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1113 1376547 : return normalizepol_lg(z-2,n+2);
1114 : }
1115 :
1116 : GEN
1117 642299 : RgXn_recip_shallow(GEN P, long n)
1118 : {
1119 642299 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1120 642302 : setvarn(Q, varn(P));
1121 642302 : return Q;
1122 : }
1123 :
1124 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1125 : GEN
1126 31353 : RgX_recip(GEN x)
1127 : {
1128 : long lx, i, j;
1129 31353 : GEN y = cgetg_copy(x, &lx);
1130 274414 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1131 31353 : return normalizepol_lg(y,lx);
1132 : }
1133 : /* shallow version */
1134 : GEN
1135 59388 : RgX_recip_shallow(GEN x)
1136 : {
1137 : long lx, i, j;
1138 59388 : GEN y = cgetg_copy(x, &lx);
1139 356118 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1140 59388 : return normalizepol_lg(y,lx);
1141 : }
1142 :
1143 : GEN
1144 5116282 : RgX_recip_i(GEN x)
1145 : {
1146 : long lx, i, j;
1147 5116282 : GEN y = cgetg_copy(x, &lx);
1148 25956147 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1149 5116284 : return y;
1150 : }
1151 : /*******************************************************************/
1152 : /* */
1153 : /* ADDITION / SUBTRACTION */
1154 : /* */
1155 : /*******************************************************************/
1156 : /* same variable */
1157 : GEN
1158 82415008 : RgX_add(GEN x, GEN y)
1159 : {
1160 82415008 : long i, lx = lg(x), ly = lg(y);
1161 : GEN z;
1162 82415008 : if (ly <= lx) {
1163 70036984 : z = cgetg(lx,t_POL); z[1] = x[1];
1164 279858586 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1165 123414531 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1166 70014051 : z = normalizepol_lg(z, lx);
1167 : } else {
1168 12378024 : z = cgetg(ly,t_POL); z[1] = y[1];
1169 49064364 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1170 36403861 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1171 12379252 : z = normalizepol_lg(z, ly);
1172 : }
1173 82406875 : return z;
1174 : }
1175 : GEN
1176 63773182 : RgX_sub(GEN x, GEN y)
1177 : {
1178 63773182 : long i, lx = lg(x), ly = lg(y);
1179 : GEN z;
1180 63773182 : if (ly <= lx) {
1181 30019997 : z = cgetg(lx,t_POL); z[1] = x[1];
1182 157251401 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1183 56390134 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1184 29997655 : z = normalizepol_lg(z, lx);
1185 : } else {
1186 33753185 : z = cgetg(ly,t_POL); z[1] = y[1];
1187 125555202 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1188 74663141 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1189 33720567 : z = normalizepol_lg(z, ly);
1190 : }
1191 63769282 : return z;
1192 : }
1193 : GEN
1194 7783081 : RgX_neg(GEN x)
1195 49950067 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1196 :
1197 : GEN
1198 22066946 : RgX_Rg_add(GEN y, GEN x)
1199 : {
1200 : GEN z;
1201 22066946 : long lz = lg(y), i;
1202 22066946 : if (lz == 2) return scalarpol(x,varn(y));
1203 19018466 : z = cgetg(lz,t_POL); z[1] = y[1];
1204 19018444 : gel(z,2) = gadd(gel(y,2),x);
1205 69288123 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1206 : /* probably useless unless lz = 3, but cannot be skipped if y is
1207 : * an inexact 0 */
1208 19018385 : return normalizepol_lg(z,lz);
1209 : }
1210 : GEN
1211 65024 : RgX_Rg_add_shallow(GEN y, GEN x)
1212 : {
1213 : GEN z;
1214 65024 : long lz = lg(y), i;
1215 65024 : if (lz == 2) return scalarpol(x,varn(y));
1216 65024 : z = cgetg(lz,t_POL); z[1] = y[1];
1217 65024 : gel(z,2) = gadd(gel(y,2),x);
1218 130158 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1219 65025 : return normalizepol_lg(z,lz);
1220 : }
1221 : GEN
1222 182218 : RgX_Rg_sub(GEN y, GEN x)
1223 : {
1224 : GEN z;
1225 182218 : long lz = lg(y), i;
1226 182218 : if (lz == 2)
1227 : { /* scalarpol(gneg(x),varn(y)) optimized */
1228 2653 : long v = varn(y);
1229 2653 : if (isrationalzero(x)) return pol_0(v);
1230 49 : z = cgetg(3,t_POL);
1231 49 : z[1] = gequal0(x)? evalvarn(v)
1232 49 : : evalvarn(v) | evalsigne(1);
1233 49 : gel(z,2) = gneg(x); return z;
1234 : }
1235 179565 : z = cgetg(lz,t_POL); z[1] = y[1];
1236 179565 : gel(z,2) = gsub(gel(y,2),x);
1237 475042 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1238 179564 : return normalizepol_lg(z,lz);
1239 : }
1240 : GEN
1241 4866602 : Rg_RgX_sub(GEN x, GEN y)
1242 : {
1243 : GEN z;
1244 4866602 : long lz = lg(y), i;
1245 4866602 : if (lz == 2) return scalarpol(x,varn(y));
1246 4865433 : z = cgetg(lz,t_POL); z[1] = y[1];
1247 4865416 : gel(z,2) = gsub(x, gel(y,2));
1248 7297226 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1249 4865391 : return normalizepol_lg(z,lz);
1250 : }
1251 : /*******************************************************************/
1252 : /* */
1253 : /* KARATSUBA MULTIPLICATION */
1254 : /* */
1255 : /*******************************************************************/
1256 : #if 0
1257 : /* to debug Karatsuba-like routines */
1258 : GEN
1259 : zx_debug_spec(GEN x, long nx)
1260 : {
1261 : GEN z = cgetg(nx+2,t_POL);
1262 : long i;
1263 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1264 : z[1] = evalsigne(1); return z;
1265 : }
1266 :
1267 : GEN
1268 : RgX_debug_spec(GEN x, long nx)
1269 : {
1270 : GEN z = cgetg(nx+2,t_POL);
1271 : long i;
1272 : for (i=0; i<nx; i++) z[i+2] = x[i];
1273 : z[1] = evalsigne(1); return z;
1274 : }
1275 : #endif
1276 :
1277 : /* generic multiplication */
1278 : GEN
1279 8905448 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1280 : {
1281 : GEN z, t;
1282 : long i;
1283 8905448 : if (nx == ny) {
1284 1535059 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1285 4950812 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1286 1535040 : return normalizepol_lg(z, nx+2);
1287 : }
1288 7370389 : if (ny < nx) {
1289 7183612 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1290 26439790 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1291 17169722 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1292 7183139 : return normalizepol_lg(z, nx+2);
1293 : } else {
1294 186777 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1295 3678393 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1296 451878 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1297 186794 : return normalizepol_lg(z, ny+2);
1298 : }
1299 : }
1300 : GEN
1301 223477 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1302 : {
1303 : GEN z, t;
1304 : long i;
1305 223477 : if (nx == ny) {
1306 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1307 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1308 12824 : return normalizepol_lg(z, nx+2);
1309 : }
1310 210653 : if (ny < nx) {
1311 208875 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1312 3729798 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1313 2382583 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1314 208875 : return normalizepol_lg(z, nx+2);
1315 : } else {
1316 1778 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1317 330904 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1318 12222 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1319 1778 : return normalizepol_lg(z, ny+2);
1320 : }
1321 : }
1322 :
1323 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1324 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1325 : * t_VECSMALL, to hold this, which can be left on stack: gerepile
1326 : * will not crash on it. The returned vector itself is not a proper GEN,
1327 : * we access the coefficients as x[i], i = 0..deg(x) */
1328 : static GEN
1329 56593473 : RgXspec_kill0(GEN x, long lx)
1330 : {
1331 56593473 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1332 : long i;
1333 192305637 : for (i=0; i <lx; i++)
1334 : {
1335 135712111 : GEN c = gel(x,i);
1336 135712111 : z[i] = (long)(isrationalzero(c)? NULL: c);
1337 : }
1338 56593526 : return z;
1339 : }
1340 :
1341 : INLINE GEN
1342 94466951 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1343 : {
1344 94466951 : pari_sp av = avma;
1345 94466951 : GEN s = NULL;
1346 : long i;
1347 :
1348 357404328 : for (i=a; i<b; i++)
1349 262946726 : if (gel(y,i) && gel(x,-i))
1350 : {
1351 184797970 : GEN t = gmul(gel(y,i), gel(x,-i));
1352 184793114 : s = s? gadd(s, t): t;
1353 : }
1354 94457602 : return s? gerepileupto(av, s): gen_0;
1355 : }
1356 :
1357 : /* assume nx >= ny > 0, return x * y * t^v */
1358 : static GEN
1359 21664404 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1360 : {
1361 : long i, lz, nz;
1362 : GEN z;
1363 :
1364 21664404 : x = RgXspec_kill0(x,nx);
1365 21664394 : y = RgXspec_kill0(y,ny);
1366 21664372 : lz = nx + ny + 1; nz = lz-2;
1367 21664372 : lz += v;
1368 21664372 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1369 33435111 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1370 57390733 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1371 41502386 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1372 35726492 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1373 21663951 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1374 : }
1375 :
1376 : /* return (x * X^d) + y. Assume d > 0 */
1377 : GEN
1378 10060228 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1379 : {
1380 : GEN x, y, xd, yd, zd;
1381 : long a, lz, nx, ny;
1382 :
1383 10060228 : if (!signe(x0)) return y0;
1384 9426050 : ny = lgpol(y0);
1385 9426025 : nx = lgpol(x0);
1386 9426145 : zd = (GEN)avma;
1387 9426145 : x = x0 + 2; y = y0 + 2; a = ny-d;
1388 9426145 : if (a <= 0)
1389 : {
1390 1478944 : lz = nx+d+2;
1391 1478944 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1392 3396499 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1393 1478944 : x = zd + a;
1394 1496652 : while (zd > x) gel(--zd,0) = gen_0;
1395 : }
1396 : else
1397 : {
1398 7947201 : xd = new_chunk(d); yd = y+d;
1399 7947200 : x = RgX_addspec_shallow(x,yd, nx,a);
1400 7947153 : lz = (a>nx)? ny+2: lg(x)+d;
1401 39605049 : x += 2; while (xd > x) *--zd = *--xd;
1402 : }
1403 22968267 : while (yd > y) *--zd = *--yd;
1404 9426097 : *--zd = x0[1];
1405 9426097 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1406 : }
1407 : GEN
1408 516050 : RgX_addmulXn(GEN x0, GEN y0, long d)
1409 : {
1410 : GEN x, y, xd, yd, zd;
1411 : long a, lz, nx, ny;
1412 :
1413 516050 : if (!signe(x0)) return RgX_copy(y0);
1414 515266 : nx = lgpol(x0);
1415 515266 : ny = lgpol(y0);
1416 515266 : zd = (GEN)avma;
1417 515266 : x = x0 + 2; y = y0 + 2; a = ny-d;
1418 515266 : if (a <= 0)
1419 : {
1420 291789 : lz = nx+d+2;
1421 291789 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1422 4294362 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1423 291789 : x = zd + a;
1424 757955 : while (zd > x) gel(--zd,0) = gen_0;
1425 : }
1426 : else
1427 : {
1428 223477 : xd = new_chunk(d); yd = y+d;
1429 223477 : x = RgX_addspec(x,yd, nx,a);
1430 223477 : lz = (a>nx)? ny+2: lg(x)+d;
1431 8430632 : x += 2; while (xd > x) *--zd = *--xd;
1432 : }
1433 2623131 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1434 515266 : *--zd = x0[1];
1435 515266 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1436 : }
1437 :
1438 : /* return x * y mod t^n */
1439 : static GEN
1440 6627959 : RgXn_mul_basecase(GEN x, GEN y, long n)
1441 : {
1442 6627959 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1443 : GEN z;
1444 6627959 : if (lx < 0) return pol_0(varn(x));
1445 6627959 : if (ly < 0) return pol_0(varn(x));
1446 6627959 : z = cgetg(lz, t_POL) + 2;
1447 6627959 : x+=2; if (lx > n) lx = n;
1448 6627959 : y+=2; if (ly > n) ly = n;
1449 6627959 : z[-1] = x[-1];
1450 6627959 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1451 6627959 : x = RgXspec_kill0(x, lx);
1452 6627959 : y = RgXspec_kill0(y, ly);
1453 : /* x:y:z [i] = term of degree i */
1454 26223885 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1455 11826184 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1456 6674031 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1457 6627959 : return normalizepol_lg(z - 2, lz);
1458 : }
1459 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1460 : static GEN
1461 9413940 : RgXn_mul2(GEN f, GEN g, long n)
1462 : {
1463 9413940 : pari_sp av = avma;
1464 : GEN fe,fo, ge,go, l,h,m;
1465 : long n0, n1;
1466 9413940 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1467 6663449 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1468 35490 : n0 = n>>1; n1 = n-n0;
1469 35490 : RgX_even_odd(f, &fe, &fo);
1470 35490 : RgX_even_odd(g, &ge, &go);
1471 35490 : l = RgXn_mul2(fe,ge,n1);
1472 35490 : h = RgXn_mul2(fo,go,n0);
1473 35490 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1474 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1475 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1476 35490 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1477 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1478 35490 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1479 35490 : m = RgX_inflate(m,2);
1480 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1481 35490 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1482 35490 : h = RgX_inflate(h,2);
1483 35490 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1484 35490 : return gerepileupto(av, h);
1485 : }
1486 : /* (f*g) \/ x^n */
1487 : static GEN
1488 1589258 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1489 : {
1490 1589258 : long d = degpol(f)+degpol(g) + 1 - n;
1491 : GEN h;
1492 1589258 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1493 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1494 : RgX_recip_i(g), d));
1495 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1496 : }
1497 :
1498 : /* (f*g) \/ x^n */
1499 : static GEN
1500 0 : RgX_sqrhigh_i2(GEN f, long n)
1501 : {
1502 0 : long d = 2*degpol(f)+ 1 - n;
1503 : GEN h;
1504 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1505 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1506 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1507 : }
1508 :
1509 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1510 : * b+2 were sent instead. na, nb = number of terms of a, b.
1511 : * Only c, c0, c1, c2 are genuine GEN.
1512 : */
1513 : GEN
1514 22393863 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1515 : {
1516 : GEN a0, c, c0;
1517 22393863 : long n0, n0a, i, v = 0;
1518 : pari_sp av;
1519 :
1520 28786889 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1521 28305915 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1522 22393838 : if (na < nb) swapspec(a,b, na,nb);
1523 22393838 : if (!nb) return pol_0(0);
1524 :
1525 22144882 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1526 480510 : RgX_shift_inplace_init(v);
1527 480521 : i = (na>>1); n0 = na-i; na = i;
1528 480521 : av = avma; a0 = a+n0; n0a = n0;
1529 1336037 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1530 :
1531 480521 : if (nb > n0)
1532 : {
1533 : GEN b0,c1,c2;
1534 : long n0b;
1535 :
1536 479122 : nb -= n0; b0 = b+n0; n0b = n0;
1537 1447101 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1538 479122 : c = RgX_mulspec(a,b,n0a,n0b);
1539 479122 : c0 = RgX_mulspec(a0,b0, na,nb);
1540 :
1541 479122 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1542 479122 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1543 :
1544 479122 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1545 479122 : c2 = RgX_sub(c1, RgX_add(c0,c));
1546 479122 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1547 : }
1548 : else
1549 : {
1550 1399 : c = RgX_mulspec(a,b,n0a,nb);
1551 1399 : c0 = RgX_mulspec(a0,b,na,nb);
1552 : }
1553 480521 : c0 = RgX_addmulXn(c0,c,n0);
1554 480521 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1555 : }
1556 :
1557 : INLINE GEN
1558 50046 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1559 : {
1560 50046 : pari_sp av = avma;
1561 50046 : GEN s = NULL;
1562 50046 : long j, l = (i+1)>>1;
1563 138971 : for (j=a; j<l; j++)
1564 : {
1565 88925 : GEN xj = gel(x,j), xx = gel(x,i-j);
1566 88925 : if (xj && xx)
1567 : {
1568 81281 : GEN t = gmul(xj, xx);
1569 81281 : s = s? gadd(s, t): t;
1570 : }
1571 : }
1572 50046 : if (s) s = gshift(s,1);
1573 50046 : if ((i&1) == 0)
1574 : {
1575 29491 : GEN t = gel(x, i>>1);
1576 29491 : if (t) {
1577 27097 : t = gsqr(t);
1578 27097 : s = s? gadd(s, t): t;
1579 : }
1580 : }
1581 50046 : return s? gerepileupto(av,s): gen_0;
1582 : }
1583 : static GEN
1584 8936 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1585 : {
1586 : long i, lz, nz;
1587 : GEN z;
1588 :
1589 8936 : if (!nx) return pol_0(0);
1590 8936 : x = RgXspec_kill0(x,nx);
1591 8936 : lz = (nx << 1) + 1, nz = lz-2;
1592 8936 : lz += v;
1593 8936 : z = cgetg(lz,t_POL) + 2;
1594 13304 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1595 38427 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1596 29491 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1597 8936 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1598 : }
1599 : /* return x^2 mod t^n */
1600 : static GEN
1601 0 : RgXn_sqr_basecase(GEN x, long n)
1602 : {
1603 0 : long i, lz = n+2, lx = lgpol(x);
1604 : GEN z;
1605 0 : if (lx < 0) return pol_0(varn(x));
1606 0 : z = cgetg(lz, t_POL);
1607 0 : z[1] = x[1];
1608 0 : x+=2; if (lx > n) lx = n;
1609 0 : x = RgXspec_kill0(x,lx);
1610 0 : z+=2;/* x:z [i] = term of degree i */
1611 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1612 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1613 0 : z -= 2; return normalizepol_lg(z, lz);
1614 : }
1615 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1616 : static GEN
1617 315 : RgXn_sqr2(GEN f, long n)
1618 : {
1619 315 : pari_sp av = avma;
1620 : GEN fe,fo, l,h,m;
1621 : long n0, n1;
1622 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1623 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1624 0 : n0 = n>>1; n1 = n-n0;
1625 0 : RgX_even_odd(f, &fe, &fo);
1626 0 : l = RgXn_sqr(fe,n1);
1627 0 : h = RgXn_sqr(fo,n0);
1628 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1629 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1630 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1631 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1632 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1633 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1634 0 : m = RgX_inflate(m,2);
1635 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1636 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1637 0 : h = RgX_inflate(h,2);
1638 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1639 0 : return gerepileupto(av, h);
1640 : }
1641 : GEN
1642 8975 : RgX_sqrspec(GEN a, long na)
1643 : {
1644 : GEN a0, c, c0, c1;
1645 8975 : long n0, n0a, i, v = 0;
1646 : pari_sp av;
1647 :
1648 11159 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1649 8975 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1650 39 : RgX_shift_inplace_init(v);
1651 39 : i = (na>>1); n0 = na-i; na = i;
1652 39 : av = avma; a0 = a+n0; n0a = n0;
1653 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1654 :
1655 39 : c = RgX_sqrspec(a,n0a);
1656 39 : c0 = RgX_sqrspec(a0,na);
1657 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1658 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1659 39 : c0 = RgX_addmulXn(c0,c,n0);
1660 39 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1661 : }
1662 :
1663 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1664 : GEN
1665 1718723 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1666 : {
1667 1718723 : GEN z = RgX_mul(A, B);
1668 1718713 : if (a < b)
1669 10444 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1670 1708269 : else if (a > b)
1671 1104578 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1672 : else
1673 603691 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1674 1718718 : return z;
1675 : }
1676 :
1677 : GEN
1678 20953705 : RgX_mul_i(GEN x, GEN y)
1679 : {
1680 20953705 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1681 20953382 : setvarn(z, varn(x)); return z;
1682 : }
1683 :
1684 : GEN
1685 8897 : RgX_sqr_i(GEN x)
1686 : {
1687 8897 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1688 8897 : setvarn(z,varn(x)); return z;
1689 : }
1690 :
1691 : /*******************************************************************/
1692 : /* */
1693 : /* DIVISION */
1694 : /* */
1695 : /*******************************************************************/
1696 : GEN
1697 7165640 : RgX_Rg_divexact(GEN x, GEN y) {
1698 7165640 : long i, lx = lg(x);
1699 : GEN z;
1700 7165640 : if (lx == 2) return gcopy(x);
1701 7114683 : switch(typ(y))
1702 : {
1703 6954414 : case t_INT:
1704 6954414 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1705 6579159 : break;
1706 5243 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1707 : }
1708 6734185 : z = cgetg(lx, t_POL); z[1] = x[1];
1709 34393061 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1710 6732592 : return z;
1711 : }
1712 : GEN
1713 30198654 : RgX_Rg_div(GEN x, GEN y) {
1714 30198654 : long i, lx = lg(x);
1715 : GEN z;
1716 30198654 : if (lx == 2) return gcopy(x);
1717 29938764 : switch(typ(y))
1718 : {
1719 20750139 : case t_INT:
1720 20750139 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1721 4136110 : break;
1722 6244 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1723 : }
1724 13318491 : z = cgetg(lx, t_POL); z[1] = x[1];
1725 49715136 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1726 13318072 : return normalizepol_lg(z, lx);
1727 : }
1728 : GEN
1729 39109 : RgX_normalize(GEN x)
1730 : {
1731 39109 : GEN z, d = NULL;
1732 39109 : long i, n = lg(x)-1;
1733 39109 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1734 39109 : if (i == 1) return pol_0(varn(x));
1735 39109 : if (i == n && isint1(d)) return x;
1736 17612 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1737 31612 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1738 17612 : gel(z,n) = Rg_get_1(d); return z;
1739 : }
1740 : GEN
1741 10458 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1742 : GEN
1743 251014 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1744 : {
1745 251014 : long l = lg(a), i;
1746 : GEN a0, z0, z;
1747 :
1748 251014 : if (l <= 3)
1749 : {
1750 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1751 0 : return pol_0(varn(a));
1752 : }
1753 251014 : z = cgetg(l-1, t_POL);
1754 251018 : z[1] = a[1];
1755 251018 : a0 = a + l-1;
1756 251018 : z0 = z + l-2; *z0 = *a0--;
1757 2955468 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1758 : {
1759 2704451 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1760 2704450 : gel(z0,0) = t;
1761 : }
1762 251017 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1763 251017 : return z;
1764 : }
1765 : /* Polynomial division x / y:
1766 : * if pr = ONLY_REM return remainder, otherwise return quotient
1767 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1768 : * if pr != NULL set *pr to remainder, as the last object on stack */
1769 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1770 : static GEN
1771 23507865 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1772 : {
1773 : pari_sp avy, av, av1;
1774 : long dx,dy,dz,i,j,sx,lr;
1775 : GEN z,p1,p2,rem,y_lead,mod,p;
1776 : GEN (*f)(GEN,GEN);
1777 :
1778 23507865 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1779 :
1780 23507865 : dy = degpol(y);
1781 23507857 : y_lead = gel(y,dy+2);
1782 23507857 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1783 : {
1784 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1785 0 : for (dy--; dy>=0; dy--)
1786 : {
1787 0 : y_lead = gel(y,dy+2);
1788 0 : if (!gequal0(y_lead)) break;
1789 : }
1790 : }
1791 23507839 : if (!dy) /* y is constant */
1792 : {
1793 6857 : if (pr == ONLY_REM) return pol_0(varn(x));
1794 6850 : z = RgX_Rg_div(x, y_lead);
1795 6850 : if (pr == ONLY_DIVIDES) return z;
1796 2636 : if (pr) *pr = pol_0(varn(x));
1797 2636 : return z;
1798 : }
1799 23500982 : dx = degpol(x);
1800 23500931 : if (dx < dy)
1801 : {
1802 3845195 : if (pr == ONLY_REM) return RgX_copy(x);
1803 348542 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1804 348521 : z = pol_0(varn(x));
1805 348521 : if (pr) *pr = RgX_copy(x);
1806 348521 : return z;
1807 : }
1808 :
1809 : /* x,y in R[X], y non constant */
1810 19655736 : av = avma;
1811 19655736 : p = NULL;
1812 19655736 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1813 : {
1814 130403 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1815 130403 : if (!z) return gc_NULL(av);
1816 130403 : z = FpX_to_mod(z, p);
1817 130403 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1818 71218 : return gerepileupto(av, z);
1819 59185 : *pr = FpX_to_mod(*pr, p);
1820 59185 : return gc_all(av, 2, &z, pr);
1821 : }
1822 19525448 : switch(typ(y_lead))
1823 : {
1824 0 : case t_REAL:
1825 0 : y_lead = ginv(y_lead);
1826 0 : f = gmul; mod = NULL;
1827 0 : break;
1828 2074 : case t_INTMOD:
1829 2074 : case t_POLMOD: y_lead = ginv(y_lead);
1830 2074 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1831 2074 : break;
1832 19523374 : default: if (gequal1(y_lead)) y_lead = NULL;
1833 19523367 : f = gdiv; mod = NULL;
1834 : }
1835 :
1836 19525441 : if (y_lead == NULL)
1837 17512415 : p2 = gel(x,dx+2);
1838 : else {
1839 : for(;;) {
1840 2013026 : p2 = f(gel(x,dx+2),y_lead);
1841 2013027 : p2 = simplify_shallow(p2);
1842 2013027 : if (!isexactzero(p2) || (--dx < 0)) break;
1843 : }
1844 2013027 : if (dx < dy) /* leading coeff of x was in fact zero */
1845 : {
1846 0 : if (pr == ONLY_DIVIDES) {
1847 0 : set_avma(av);
1848 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1849 : }
1850 0 : if (pr == ONLY_REM)
1851 : {
1852 0 : if (dx < 0)
1853 0 : return gerepilecopy(av, scalarpol(p2, varn(x)));
1854 : else
1855 : {
1856 : GEN t;
1857 0 : set_avma(av);
1858 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1859 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1860 0 : return t;
1861 : }
1862 : }
1863 0 : if (pr) /* cf ONLY_REM above */
1864 : {
1865 0 : if (dx < 0)
1866 : {
1867 0 : p2 = gclone(p2);
1868 0 : set_avma(av);
1869 0 : z = pol_0(varn(x));
1870 0 : x = scalarpol(p2, varn(x));
1871 0 : gunclone(p2);
1872 : }
1873 : else
1874 : {
1875 : GEN t;
1876 0 : set_avma(av);
1877 0 : z = pol_0(varn(x));
1878 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1879 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1880 0 : x = t;
1881 : }
1882 0 : *pr = x;
1883 : }
1884 : else
1885 : {
1886 0 : set_avma(av);
1887 0 : z = pol_0(varn(x));
1888 : }
1889 0 : return z;
1890 : }
1891 : }
1892 : /* dx >= dy */
1893 19525442 : avy = avma;
1894 19525442 : dz = dx-dy;
1895 19525442 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1896 19525416 : x += 2;
1897 19525416 : z += 2;
1898 19525416 : y += 2;
1899 19525416 : gel(z,dz) = gcopy(p2);
1900 :
1901 54074407 : for (i=dx-1; i>=dy; i--)
1902 : {
1903 34549072 : av1=avma; p1=gel(x,i);
1904 1140457847 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1905 34543903 : if (y_lead) p1 = simplify(f(p1,y_lead));
1906 :
1907 34543903 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1908 : else
1909 24553438 : p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1910 34548599 : gel(z,i-dy) = p1;
1911 : }
1912 19525335 : if (!pr) return gerepileupto(av,z-2);
1913 :
1914 12548563 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1915 12913229 : for (sx=0; ; i--)
1916 : {
1917 12913229 : p1 = gel(x,i);
1918 : /* we always enter this loop at least once */
1919 32430445 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1920 12911977 : if (mod && avma==av1) p1 = gmul(p1,mod);
1921 12911977 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1922 3587623 : if (!isexactzero(p1)) break;
1923 3437564 : if (!i) break;
1924 364617 : set_avma(av1);
1925 : }
1926 12547723 : if (pr == ONLY_DIVIDES)
1927 : {
1928 7931 : if (sx) return gc_NULL(av);
1929 7910 : set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1930 : }
1931 12539792 : lr=i+3; rem -= lr;
1932 12539792 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1933 12496187 : else p1 = gerepileupto((pari_sp)rem,p1);
1934 12540577 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1935 12540577 : rem[1] = z[-1];
1936 12540577 : rem += 2;
1937 12540577 : gel(rem,i) = p1;
1938 17268372 : for (i--; i>=0; i--)
1939 : {
1940 4727827 : av1=avma; p1 = gel(x,i);
1941 13589874 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1942 4727283 : if (mod && avma==av1) p1 = gmul(p1,mod);
1943 4727491 : gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1944 : }
1945 12540545 : rem -= 2;
1946 12540545 : if (!sx) (void)normalizepol_lg(rem, lr);
1947 12540750 : if (pr == ONLY_REM) return gerepileupto(av,rem);
1948 6761860 : z -= 2;
1949 : {
1950 6761860 : GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1951 6761860 : gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1952 : }
1953 : }
1954 :
1955 : GEN
1956 14232365 : RgX_divrem(GEN x, GEN y, GEN *pr)
1957 : {
1958 14232365 : if (pr == ONLY_REM) return RgX_rem(x, y);
1959 14232365 : return RgX_divrem_i(x, y, pr);
1960 : }
1961 :
1962 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1963 : GEN
1964 150959 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1965 : {
1966 : long vx, dx, dy, dz, i, j, sx, lr;
1967 : pari_sp av0, av, tetpil;
1968 : GEN z,p1,rem,lead;
1969 :
1970 150959 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1971 150959 : vx = varn(x);
1972 150959 : dx = degpol(x);
1973 150959 : dy = degpol(y);
1974 150959 : if (dx < dy)
1975 : {
1976 37570 : if (pr)
1977 : {
1978 37542 : av0 = avma; x = RgXQX_red(x, T);
1979 37542 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1980 37535 : if (pr == ONLY_REM) return x;
1981 0 : *pr = x;
1982 : }
1983 28 : return pol_0(vx);
1984 : }
1985 113389 : lead = leading_coeff(y);
1986 113389 : if (!dy) /* y is constant */
1987 : {
1988 602 : if (pr && pr != ONLY_DIVIDES)
1989 : {
1990 0 : if (pr == ONLY_REM) return pol_0(vx);
1991 0 : *pr = pol_0(vx);
1992 : }
1993 602 : if (gequal1(lead)) return RgX_copy(x);
1994 0 : av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
1995 0 : return gerepile(av0,tetpil,RgXQX_red(x,T));
1996 : }
1997 112787 : av0 = avma; dz = dx-dy;
1998 112787 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
1999 112787 : set_avma(av0);
2000 112787 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2001 112787 : x += 2; y += 2; z += 2;
2002 :
2003 112787 : p1 = gel(x,dx); av = avma;
2004 112787 : gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2005 575780 : for (i=dx-1; i>=dy; i--)
2006 : {
2007 462993 : av=avma; p1=gel(x,i);
2008 2406142 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2009 462962 : if (lead) p1 = gmul(grem(p1, T), lead);
2010 462960 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
2011 : }
2012 112787 : if (!pr) { guncloneNULL(lead); return z-2; }
2013 :
2014 111513 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2015 187594 : for (sx=0; ; i--)
2016 : {
2017 187594 : p1 = gel(x,i);
2018 649629 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2019 187596 : tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2020 103297 : if (!i) break;
2021 76081 : set_avma(av);
2022 : }
2023 111513 : if (pr == ONLY_DIVIDES)
2024 : {
2025 27202 : guncloneNULL(lead);
2026 27202 : if (sx) return gc_NULL(av0);
2027 24640 : return gc_const((pari_sp)rem, z-2);
2028 : }
2029 84311 : lr=i+3; rem -= lr;
2030 84311 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2031 84311 : rem[1] = z[-1];
2032 84311 : p1 = gerepile((pari_sp)rem,tetpil,p1);
2033 84311 : rem += 2; gel(rem,i) = p1;
2034 184203 : for (i--; i>=0; i--)
2035 : {
2036 99892 : av=avma; p1 = gel(x,i);
2037 334034 : for (j=0; j<=i && j<=dz; j++)
2038 234142 : p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2039 99892 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
2040 : }
2041 84311 : rem -= 2;
2042 84311 : guncloneNULL(lead);
2043 84311 : if (!sx) (void)normalizepol_lg(rem, lr);
2044 84311 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
2045 161 : *pr = rem; return z-2;
2046 : }
2047 :
2048 : /*******************************************************************/
2049 : /* */
2050 : /* PSEUDO-DIVISION */
2051 : /* */
2052 : /*******************************************************************/
2053 : INLINE GEN
2054 4051726 : rem(GEN c, GEN T)
2055 : {
2056 4051726 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2057 4051725 : return c;
2058 : }
2059 :
2060 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2061 : int
2062 17321 : ZXQX_dvd(GEN x, GEN y, GEN T)
2063 : {
2064 : long dx, dy, i, T_ismonic;
2065 17321 : pari_sp av = avma, av2;
2066 : GEN y_lead;
2067 :
2068 17321 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2069 17321 : dy = degpol(y); y_lead = gel(y,dy+2);
2070 17321 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2071 : /* if monic, no point in using pseudo-division */
2072 17321 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2073 14780 : T_ismonic = gequal1(leading_coeff(T));
2074 14780 : dx = degpol(x);
2075 14780 : if (dx < dy) return !signe(x);
2076 14780 : (void)new_chunk(2);
2077 14780 : x = RgX_recip_i(x)+2;
2078 14780 : y = RgX_recip_i(y)+2;
2079 : /* pay attention to sparse divisors */
2080 30061 : for (i = 1; i <= dy; i++)
2081 15281 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2082 14780 : av2 = avma;
2083 : for (;;)
2084 72995 : {
2085 87775 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2086 87774 : x0 = gneg(x0);
2087 87775 : m = gcdii(cx, y0);
2088 87775 : if (!equali1(m))
2089 : {
2090 85273 : x0 = gdiv(x0, m);
2091 85272 : y0 = diviiexact(y0, m);
2092 85271 : if (equali1(y0)) y0 = NULL;
2093 : }
2094 179321 : for (i=1; i<=dy; i++)
2095 : {
2096 91546 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2097 91546 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2098 91548 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2099 91548 : gel(x,i) = c;
2100 : }
2101 688126 : for ( ; i<=dx; i++)
2102 : {
2103 600351 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2104 600351 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2105 600351 : gel(x,i) = c;
2106 : }
2107 102321 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2108 87775 : if (dx < dy) break;
2109 72995 : if (gc_needed(av2,1))
2110 : {
2111 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2112 28 : gerepilecoeffs(av2,x,dx+1);
2113 : }
2114 : }
2115 14780 : return gc_bool(av, dx < 0);
2116 : }
2117 :
2118 : /* T either NULL or a t_POL. */
2119 : GEN
2120 630477 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2121 : {
2122 630477 : long vx = varn(x), dx, dy, dz, i, lx, p;
2123 630477 : pari_sp av = avma, av2;
2124 : GEN y_lead;
2125 :
2126 630477 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2127 630477 : dy = degpol(y); y_lead = gel(y,dy+2);
2128 : /* if monic, no point in using pseudo-division */
2129 630479 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2130 581068 : dx = degpol(x);
2131 581068 : if (dx < dy) return RgX_copy(x);
2132 581061 : (void)new_chunk(2);
2133 581061 : x = RgX_recip_i(x)+2;
2134 581063 : y = RgX_recip_i(y)+2;
2135 : /* pay attention to sparse divisors */
2136 1957725 : for (i = 1; i <= dy; i++)
2137 1376663 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2138 581062 : dz = dx-dy; p = dz+1;
2139 581062 : av2 = avma;
2140 : for (;;)
2141 : {
2142 1249426 : gel(x,0) = gneg(gel(x,0)); p--;
2143 4201105 : for (i=1; i<=dy; i++)
2144 : {
2145 2951690 : GEN c = gmul(y_lead, gel(x,i));
2146 2951616 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2147 2951680 : gel(x,i) = rem(c, T);
2148 : }
2149 2187451 : for ( ; i<=dx; i++)
2150 : {
2151 938031 : GEN c = gmul(y_lead, gel(x,i));
2152 938023 : gel(x,i) = rem(c, T);
2153 : }
2154 1295808 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2155 1249420 : if (dx < dy) break;
2156 668364 : if (gc_needed(av2,1))
2157 : {
2158 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2159 0 : gerepilecoeffs(av2,x,dx+1);
2160 : }
2161 : }
2162 581056 : if (dx < 0) return pol_0(vx);
2163 579873 : lx = dx+3; x -= 2;
2164 579873 : x[0] = evaltyp(t_POL) | _evallg(lx);
2165 579873 : x[1] = evalsigne(1) | evalvarn(vx);
2166 579873 : x = RgX_recip_i(x);
2167 579874 : if (p)
2168 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2169 28664 : GEN t = y_lead;
2170 28664 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2171 0 : t = RgXQ_powu(t, p, T);
2172 : else
2173 28664 : t = gpowgs(t, p);
2174 88718 : for (i=2; i<lx; i++)
2175 : {
2176 60054 : GEN c = gmul(gel(x,i), t);
2177 60054 : gel(x,i) = rem(c,T);
2178 : }
2179 28664 : if (!T) return gerepileupto(av, x);
2180 : }
2181 551210 : return gerepilecopy(av, x);
2182 : }
2183 :
2184 : GEN
2185 630476 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2186 :
2187 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2188 : GEN
2189 12098 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2190 : {
2191 12098 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2192 12098 : pari_sp av = avma, av2;
2193 : GEN z, r, ypow, y_lead;
2194 :
2195 12098 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2196 12098 : dy = degpol(y); y_lead = gel(y,dy+2);
2197 12098 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2198 7603 : dx = degpol(x);
2199 7603 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2200 7603 : if (dx == dy)
2201 : {
2202 98 : GEN x_lead = gel(x,lg(x)-1);
2203 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2204 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2205 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2206 98 : *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
2207 : }
2208 7505 : (void)new_chunk(2);
2209 7505 : x = RgX_recip_i(x)+2;
2210 7505 : y = RgX_recip_i(y)+2;
2211 : /* pay attention to sparse divisors */
2212 39084 : for (i = 1; i <= dy; i++)
2213 31579 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2214 7505 : dz = dx-dy; p = dz+1;
2215 7505 : lz = dz+3;
2216 7505 : z = cgetg(lz, t_POL);
2217 7505 : z[1] = evalsigne(1) | evalvarn(vx);
2218 28808 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2219 7505 : ypow = new_chunk(dz+1);
2220 7505 : gel(ypow,0) = gen_1;
2221 7505 : gel(ypow,1) = y_lead;
2222 13798 : for (i=2; i<=dz; i++)
2223 : {
2224 6293 : GEN c = gmul(gel(ypow,i-1), y_lead);
2225 6293 : gel(ypow,i) = rem(c,T);
2226 : }
2227 7505 : av2 = avma;
2228 7505 : for (iz=2;;)
2229 : {
2230 15731 : p--;
2231 15731 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2232 70167 : for (i=1; i<=dy; i++)
2233 : {
2234 54436 : GEN c = gmul(y_lead, gel(x,i));
2235 54436 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2236 54436 : gel(x,i) = rem(c, T);
2237 : }
2238 41310 : for ( ; i<=dx; i++)
2239 : {
2240 25579 : GEN c = gmul(y_lead, gel(x,i));
2241 25579 : gel(x,i) = rem(c,T);
2242 : }
2243 15731 : x++; dx--;
2244 21303 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2245 15731 : if (dx < dy) break;
2246 8226 : if (gc_needed(av2,1))
2247 : {
2248 0 : GEN X = x-2;
2249 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2250 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2251 0 : gerepileall(av2,2, &X, &z); x = X+2;
2252 : }
2253 : }
2254 14050 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2255 7505 : if (dx < 0)
2256 182 : x = pol_0(vx);
2257 : else
2258 : {
2259 7323 : lx = dx+3; x -= 2;
2260 7323 : x[0] = evaltyp(t_POL) | _evallg(lx);
2261 7323 : x[1] = evalsigne(1) | evalvarn(vx);
2262 7323 : x = RgX_recip_i(x);
2263 : }
2264 7505 : z = RgX_recip_i(z);
2265 7505 : r = x;
2266 7505 : if (p)
2267 : {
2268 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2269 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2270 : }
2271 7505 : *ptr = r; return gc_all(av, 2, &z, ptr);
2272 : }
2273 : GEN
2274 11860 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2275 11860 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2276 :
2277 : GEN
2278 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2279 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2280 : GEN
2281 524704102 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2282 : GEN
2283 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2284 : GEN
2285 26411 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2286 : GEN
2287 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2288 : GEN
2289 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2290 :
2291 : GEN
2292 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2293 :
2294 : GEN
2295 0 : RgXQX_powers(GEN P, long n, GEN T)
2296 : {
2297 0 : GEN v = cgetg(n+2, t_VEC);
2298 : long i;
2299 0 : gel(v, 1) = pol_1(varn(T));
2300 0 : if (n==0) return v;
2301 0 : gel(v, 2) = gcopy(P);
2302 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2303 0 : return v;
2304 : }
2305 :
2306 : static GEN
2307 617696 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2308 : static GEN
2309 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2310 : static GEN
2311 67208 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2312 : static GEN
2313 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2314 : static GEN
2315 275342 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2316 : static GEN
2317 1030578 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2318 : static GEN
2319 1009081 : _one(void *data) { return pol_1(varn((GEN)data)); }
2320 : static GEN
2321 1155 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2322 : static GEN
2323 717966 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2324 :
2325 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2326 : _mul, _sqr, _one, _zero };
2327 :
2328 : GEN
2329 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2330 : {
2331 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2332 : }
2333 :
2334 : GEN
2335 413085 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2336 : {
2337 413085 : int use_sqr = 2*degpol(x) >= degpol(T);
2338 413085 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2339 : }
2340 :
2341 : /* mod X^n */
2342 : struct modXn {
2343 : long v; /* varn(X) */
2344 : long n;
2345 : } ;
2346 : static GEN
2347 11893 : _sqrXn(void *data, GEN x) {
2348 11893 : struct modXn *S = (struct modXn*)data;
2349 11893 : return RgXn_sqr(x, S->n);
2350 : }
2351 : static GEN
2352 4528 : _mulXn(void *data, GEN x, GEN y) {
2353 4528 : struct modXn *S = (struct modXn*)data;
2354 4528 : return RgXn_mul(x,y, S->n);
2355 : }
2356 : static GEN
2357 1939 : _oneXn(void *data) {
2358 1939 : struct modXn *S = (struct modXn*)data;
2359 1939 : return pol_1(S->v);
2360 : }
2361 : static GEN
2362 0 : _zeroXn(void *data) {
2363 0 : struct modXn *S = (struct modXn*)data;
2364 0 : return pol_0(S->v);
2365 : }
2366 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2367 : _oneXn, _zeroXn };
2368 :
2369 : GEN
2370 483 : RgXn_powers(GEN x, long m, long n)
2371 : {
2372 483 : long d = degpol(x);
2373 483 : int use_sqr = (d<<1) >= n;
2374 : struct modXn S;
2375 483 : S.v = varn(x); S.n = n;
2376 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2377 : }
2378 :
2379 : GEN
2380 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2381 : {
2382 : struct modXn S;
2383 : long v;
2384 2286 : if (n == 0) return x;
2385 2286 : v = RgX_valrem(x, &x);
2386 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2387 2265 : S.v = varn(x); S.n = n;
2388 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2389 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2390 2265 : return x;
2391 : }
2392 : GEN
2393 0 : RgXn_powu(GEN x, ulong m, long n)
2394 : {
2395 : pari_sp av;
2396 0 : if (n == 0) return gcopy(x);
2397 0 : av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
2398 : }
2399 :
2400 : GEN
2401 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2402 : {
2403 : struct modXn S;
2404 966 : S.v = varn(gel(x,2)); S.n = n;
2405 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2406 : }
2407 :
2408 : GEN
2409 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2410 : {
2411 0 : int use_sqr = 2*degpol(x) >= n;
2412 : struct modXn S;
2413 0 : S.v = varn(x); S.n = n;
2414 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2415 : }
2416 :
2417 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2418 : GEN
2419 5159 : RgXn_eval(GEN Q, GEN x, long n)
2420 : {
2421 5159 : long d = degpol(x);
2422 : int use_sqr;
2423 : struct modXn S;
2424 5159 : if (d == 1 && isrationalzero(gel(x,2)))
2425 : {
2426 5152 : GEN y = RgX_unscale(Q, gel(x,3));
2427 5152 : setvarn(y, varn(x)); return y;
2428 : }
2429 7 : S.v = varn(x);
2430 7 : S.n = n;
2431 7 : use_sqr = (d<<1) >= n;
2432 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2433 : }
2434 :
2435 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2436 : static GEN
2437 2614491 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2438 : {
2439 2614491 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2440 2614491 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2441 : }
2442 :
2443 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2444 : static GEN
2445 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2446 : {
2447 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2448 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2449 : }
2450 :
2451 : static GEN
2452 2188922 : RgXn_div_gen(GEN g, GEN f, long e)
2453 : {
2454 : pari_sp av;
2455 : ulong mask;
2456 : GEN W, a;
2457 2188922 : long v = varn(f), n = 1;
2458 :
2459 2188922 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2460 2188915 : a = ginv(gel(f,2));
2461 2188912 : if (e == 1 && !g) return scalarpol(a, v);
2462 2184026 : else if (e == 2 && !g)
2463 : {
2464 : GEN b;
2465 836930 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2466 246251 : b = gneg(b);
2467 246249 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2468 246251 : return deg1pol(b, a, v);
2469 : }
2470 1347096 : av = avma;
2471 1347096 : W = scalarpol_shallow(a,v);
2472 1347096 : mask = quadratic_prec_mask(e);
2473 3955623 : while (mask > 1)
2474 : {
2475 : GEN u, fr;
2476 2608527 : long n2 = n;
2477 2608527 : n<<=1; if (mask & 1) n--;
2478 2608527 : mask >>= 1;
2479 2608527 : fr = RgXn_red_shallow(f, n);
2480 2608527 : if (mask>1 || !g)
2481 : {
2482 1324547 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2483 1324547 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2484 : }
2485 : else
2486 : {
2487 1283980 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2488 1283980 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2489 1283980 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2490 : }
2491 2608527 : if (gc_needed(av,2))
2492 : {
2493 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2494 0 : W = gerepileupto(av, W);
2495 : }
2496 : }
2497 1347096 : return W;
2498 : }
2499 :
2500 : static GEN
2501 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2502 : {
2503 : GEN r;
2504 119 : if (lgefint(p) == 3)
2505 : {
2506 119 : ulong pp = uel(p, 2);
2507 119 : if (pp == 2)
2508 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2509 : else
2510 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2511 : }
2512 : else
2513 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2514 119 : return FpX_to_mod(r, p);
2515 : }
2516 :
2517 : static GEN
2518 7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2519 : {
2520 7 : GEN r, T = RgX_to_FpX(pol, p);
2521 7 : if (signe(T) == 0) pari_err_OP("/", x, y);
2522 7 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2523 7 : return FpXQX_to_mod(r, T, p);
2524 : }
2525 :
2526 : static GEN
2527 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2528 : {
2529 : GEN r;
2530 91 : if (lgefint(p) == 3)
2531 : {
2532 91 : ulong pp = uel(p, 2);
2533 91 : if (pp == 2)
2534 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2535 : else
2536 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2537 : }
2538 : else
2539 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2540 91 : return FpX_to_mod(r, p);
2541 : }
2542 :
2543 : static GEN
2544 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2545 : {
2546 0 : GEN r, T = RgX_to_FpX(pol, p);
2547 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2548 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2549 0 : return FpXQX_to_mod(r, T, p);
2550 : }
2551 :
2552 : #define code(t1,t2) ((t1 << 6) | t2)
2553 :
2554 : static GEN
2555 905035 : RgXn_inv_fast(GEN x, long e)
2556 : {
2557 : GEN p, pol;
2558 : long pa;
2559 905035 : long t = RgX_type(x,&p,&pol,&pa);
2560 905034 : switch(t)
2561 : {
2562 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2563 0 : case code(t_POLMOD, t_INTMOD):
2564 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2565 904943 : default: return NULL;
2566 : }
2567 : }
2568 :
2569 : static GEN
2570 1284106 : RgXn_div_fast(GEN x, GEN y, long e)
2571 : {
2572 : GEN p, pol;
2573 : long pa;
2574 1284106 : long t = RgX_type2(x,y,&p,&pol,&pa);
2575 1284106 : switch(t)
2576 : {
2577 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2578 7 : case code(t_POLMOD, t_INTMOD):
2579 7 : return RgXn_div_FpXQX(x, y, e, pol, p);
2580 1283980 : default: return NULL;
2581 : }
2582 : }
2583 : #undef code
2584 :
2585 : GEN
2586 1284106 : RgXn_div_i(GEN g, GEN f, long e)
2587 : {
2588 1284106 : GEN h = RgXn_div_fast(g, f, e);
2589 1284106 : if (h) return h;
2590 1283980 : return RgXn_div_gen(g, f, e);
2591 : }
2592 :
2593 : GEN
2594 584443 : RgXn_div(GEN g, GEN f, long e)
2595 : {
2596 584443 : pari_sp av = avma;
2597 584443 : return gerepileupto(av, RgXn_div_i(g, f, e));
2598 : }
2599 :
2600 : GEN
2601 905036 : RgXn_inv_i(GEN f, long e)
2602 : {
2603 905036 : GEN h = RgXn_inv_fast(f, e);
2604 905034 : if (h) return h;
2605 904943 : return RgXn_div_gen(NULL, f, e);
2606 : }
2607 :
2608 : GEN
2609 807457 : RgXn_inv(GEN f, long e)
2610 : {
2611 807457 : pari_sp av = avma;
2612 807457 : return gerepileupto(av, RgXn_inv_i(f, e));
2613 : }
2614 :
2615 : /* intformal(x^n*S) / x^(n+1) */
2616 : static GEN
2617 56392 : RgX_integXn(GEN x, long n)
2618 114387 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2619 :
2620 : GEN
2621 52927 : RgXn_expint(GEN h, long e)
2622 : {
2623 52927 : pari_sp av = avma, av2;
2624 52927 : long v = varn(h), n;
2625 52927 : GEN f = pol_1(v), g;
2626 : ulong mask;
2627 :
2628 52927 : if (!signe(h)) return f;
2629 50435 : g = pol_1(v);
2630 50434 : n = 1; mask = quadratic_prec_mask(e);
2631 50434 : av2 = avma;
2632 56391 : for (;mask>1;)
2633 : {
2634 : GEN u, w;
2635 56391 : long n2 = n;
2636 56391 : n<<=1; if (mask & 1) n--;
2637 56391 : mask >>= 1;
2638 56391 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2639 56392 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2640 56392 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2641 56392 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2642 56390 : if (mask<=1) break;
2643 5957 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2644 5957 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2645 5957 : if (gc_needed(av2,2))
2646 : {
2647 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2648 0 : gerepileall(av2, 2, &f, &g);
2649 : }
2650 : }
2651 50433 : return gerepileupto(av, f);
2652 : }
2653 :
2654 : GEN
2655 0 : RgXn_exp(GEN h, long e)
2656 : {
2657 0 : long d = degpol(h);
2658 0 : if (d < 0) return pol_1(varn(h));
2659 0 : if (!d || !gequal0(gel(h,2)))
2660 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2661 0 : return RgXn_expint(RgX_deriv(h), e);
2662 : }
2663 :
2664 : GEN
2665 154 : RgXn_reverse(GEN f, long e)
2666 : {
2667 154 : pari_sp av = avma, av2;
2668 : ulong mask;
2669 : GEN fi, a, df, W, an;
2670 154 : long v = varn(f), n=1;
2671 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2672 0 : pari_err_INV("serreverse",f);
2673 154 : fi = ginv(gel(f,3));
2674 154 : a = deg1pol_shallow(fi,gen_0,v);
2675 154 : if (e <= 2) return gerepilecopy(av, a);
2676 133 : W = scalarpol(fi,v);
2677 133 : df = RgX_deriv(f);
2678 133 : mask = quadratic_prec_mask(e);
2679 133 : av2 = avma;
2680 616 : for (;mask>1;)
2681 : {
2682 : GEN u, fa, fr;
2683 483 : long n2 = n, rt;
2684 483 : n<<=1; if (mask & 1) n--;
2685 483 : mask >>= 1;
2686 483 : fr = RgXn_red_shallow(f, n);
2687 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2688 483 : an = RgXn_powers(a, rt, n);
2689 483 : if (n>1)
2690 : {
2691 483 : long n4 = (n2+1)>>1;
2692 483 : GEN dfr = RgXn_red_shallow(df, n2);
2693 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2694 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2695 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2696 : }
2697 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2698 483 : fa = RgX_shift(fa, -n2);
2699 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2700 483 : if (gc_needed(av2,2))
2701 : {
2702 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2703 0 : gerepileall(av2, 2, &a, &W);
2704 : }
2705 : }
2706 133 : return gerepileupto(av, a);
2707 : }
2708 :
2709 : GEN
2710 7 : RgXn_sqrt(GEN h, long e)
2711 : {
2712 7 : pari_sp av = avma, av2;
2713 7 : long v = varn(h), n = 1;
2714 7 : GEN f = scalarpol(gen_1, v), df = f;
2715 7 : ulong mask = quadratic_prec_mask(e);
2716 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2717 0 : pari_err_SQRTN("RgXn_sqrt",h);
2718 7 : av2 = avma;
2719 : while(1)
2720 7 : {
2721 14 : long n2 = n, m;
2722 : GEN g;
2723 14 : n<<=1; if (mask & 1) n--;
2724 14 : mask >>= 1;
2725 14 : m = n-n2;
2726 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2727 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2728 14 : if (mask==1) return gerepileupto(av, f);
2729 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2730 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2731 7 : if (gc_needed(av2,2))
2732 : {
2733 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2734 0 : gerepileall(av2, 2, &f, &df);
2735 : }
2736 : }
2737 : }
2738 :
2739 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2740 : GEN
2741 116554 : RgXQ_powu(GEN x, ulong n, GEN T)
2742 : {
2743 116554 : pari_sp av = avma;
2744 :
2745 116554 : if (!n) return pol_1(varn(x));
2746 68385 : if (n == 1) return RgX_copy(x);
2747 19831 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2748 19831 : return gerepilecopy(av, x);
2749 : }
2750 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2751 : GEN
2752 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2753 : {
2754 : pari_sp av;
2755 102160 : long s = signe(n);
2756 :
2757 102160 : if (!s) return pol_1(varn(x));
2758 102160 : if (is_pm1(n) == 1)
2759 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2760 13853 : av = avma;
2761 13853 : if (s < 0) x = RgXQ_inv(x, T);
2762 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2763 13853 : return gerepilecopy(av, x);
2764 : }
2765 : static GEN
2766 200577 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2767 : static GEN
2768 111189 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2769 :
2770 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2771 : GEN
2772 11529 : ZXQ_powers(GEN x, long l, GEN T)
2773 : {
2774 11529 : int use_sqr = 2*degpol(x) >= degpol(T);
2775 11529 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2776 : }
2777 :
2778 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2779 : GEN
2780 168509 : ZXQ_powu(GEN x, ulong n, GEN T)
2781 : {
2782 168509 : pari_sp av = avma;
2783 :
2784 168509 : if (!n) return pol_1(varn(x));
2785 168509 : if (n == 1) return ZX_copy(x);
2786 113874 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2787 113876 : return gerepilecopy(av, x);
2788 : }
2789 :
2790 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2791 : GEN
2792 8561 : RgXQ_powers(GEN x, long l, GEN T)
2793 : {
2794 8561 : int use_sqr = 2*degpol(x) >= degpol(T);
2795 8561 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2796 : }
2797 :
2798 : GEN
2799 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2800 : {
2801 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2802 : }
2803 :
2804 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2805 : GEN
2806 8974 : QXQ_powers(GEN a, long n, GEN T)
2807 : {
2808 : GEN den, v;
2809 8974 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2810 8946 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2811 : /* den*a integral; v[i+1] = (den*a)^i in K */
2812 8946 : if (den)
2813 : { /* restore denominators */
2814 5526 : GEN d = den;
2815 : long i;
2816 5526 : gel(v,2) = a;
2817 27125 : for (i=3; i<=n+1; i++) {
2818 21599 : d = mulii(d,den);
2819 21599 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2820 : }
2821 : }
2822 8946 : return v;
2823 : }
2824 :
2825 : static GEN
2826 3010 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2827 : {
2828 3010 : long l, i, m = 0;
2829 : GEN dz, z;
2830 3010 : GEN V = cgetg_copy(v, &l);
2831 10206 : for (i = imin; i < l; i++)
2832 : {
2833 7196 : GEN c = gel(v, i);
2834 7196 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2835 : }
2836 3010 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2837 3241 : for (i = 1; i < imin; i++) V[i] = v[i];
2838 10206 : for (i = imin; i < l; i++)
2839 : {
2840 7196 : GEN c = gel(v,i);
2841 7196 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2842 7196 : gel(V,i) = c;
2843 : }
2844 3010 : return V;
2845 : }
2846 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2847 : GEN
2848 2779 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2849 2779 : { return do_QXQ_eval(v, 1, a, T); }
2850 :
2851 : GEN
2852 231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2853 231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2854 :
2855 : GEN
2856 2947 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2857 : {
2858 2947 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2859 : }
2860 :
2861 : GEN
2862 5144 : RgXQ_norm(GEN x, GEN T)
2863 : {
2864 : pari_sp av;
2865 5144 : long dx = degpol(x);
2866 : GEN L, y;
2867 5144 : if (degpol(T)==0) return gpowgs(x,0);
2868 5137 : av = avma; y = resultant(T, x);
2869 5137 : L = leading_coeff(T);
2870 5137 : if (gequal1(L) || !signe(x)) return y;
2871 0 : return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2872 : }
2873 :
2874 : GEN
2875 476 : RgXQ_trace(GEN x, GEN T)
2876 : {
2877 476 : pari_sp av = avma;
2878 : GEN dT, z;
2879 : long n;
2880 476 : if (degpol(T)==0) return gmulgs(x,0);
2881 469 : dT = RgX_deriv(T); n = degpol(dT);
2882 469 : z = RgXQ_mul(x, dT, T);
2883 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2884 420 : return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2885 : }
2886 :
2887 : GEN
2888 2787616 : RgX_blocks(GEN P, long n, long m)
2889 : {
2890 2787616 : GEN z = cgetg(m+1,t_VEC);
2891 2787616 : long i,j, k=2, l = lg(P);
2892 8556426 : for(i=1; i<=m; i++)
2893 : {
2894 5768810 : GEN zi = cgetg(n+2,t_POL);
2895 5768810 : zi[1] = P[1];
2896 5768810 : gel(z,i) = zi;
2897 16990254 : for(j=2; j<n+2; j++)
2898 11221444 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2899 5768810 : zi = RgX_renormalize_lg(zi, n+2);
2900 : }
2901 2787616 : return z;
2902 : }
2903 :
2904 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2905 : void
2906 49718153 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2907 : {
2908 49718153 : long n = degpol(p), v = varn(p), n0, n1, i;
2909 : GEN p0, p1;
2910 :
2911 49718250 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2912 :
2913 49624973 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2914 49624973 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2915 49627671 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2916 146667839 : for (i=0; i<n1; i++)
2917 : {
2918 97037698 : p0[2+i] = p[2+(i<<1)];
2919 97037698 : p1[2+i] = p[3+(i<<1)];
2920 : }
2921 49630141 : if (n1 != n0)
2922 17988052 : p0[2+i] = p[2+(i<<1)];
2923 49630141 : *pe = normalizepol(p0);
2924 49633336 : *po = normalizepol(p1);
2925 : }
2926 :
2927 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2928 : GEN
2929 43631 : RgX_splitting(GEN p, long k)
2930 : {
2931 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
2932 : GEN r;
2933 :
2934 43631 : m = n/k;
2935 43631 : r = cgetg(k+1,t_VEC);
2936 234059 : for(i=1; i<=k; i++)
2937 : {
2938 190428 : gel(r,i) = cgetg(m+3, t_POL);
2939 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2940 : }
2941 571991 : for (j=1, i=0, l=2; i<=n; i++)
2942 : {
2943 528360 : gmael(r,j,l) = gel(p,2+i);
2944 528360 : if (j==k) { j=1; l++; } else j++;
2945 : }
2946 234059 : for(i=1; i<=k; i++)
2947 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2948 43631 : return r;
2949 : }
2950 :
2951 : /*******************************************************************/
2952 : /* */
2953 : /* Kronecker form */
2954 : /* */
2955 : /*******************************************************************/
2956 :
2957 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2958 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2959 : * normalized coefficients */
2960 : GEN
2961 187323 : Kronecker_to_mod(GEN z, GEN T)
2962 : {
2963 187323 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2964 187323 : GEN x, t = cgetg(N,t_POL);
2965 187323 : t[1] = T[1];
2966 187323 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2967 187323 : x[1] = z[1];
2968 187323 : T = RgX_copy(T);
2969 1406590 : for (i=2; i<lx+2; i++, z+= N-2)
2970 : {
2971 5793796 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2972 1219267 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2973 : }
2974 187323 : N = (l-2) % (N-2) + 2;
2975 475688 : for (j=2; j<N; j++) t[j] = z[j];
2976 187323 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2977 187323 : return normalizepol_lg(x, i+1);
2978 : }
2979 :
2980 : /*******************************************************************/
2981 : /* */
2982 : /* Domain detection */
2983 : /* */
2984 : /*******************************************************************/
2985 :
2986 : static GEN
2987 666063 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2988 : {
2989 666063 : pari_sp av = avma;
2990 : GEN r;
2991 666063 : if (lgefint(p) == 3)
2992 : {
2993 614707 : ulong pp = uel(p, 2);
2994 614707 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2995 : RgX_to_Flx(y, pp), pp));
2996 : }
2997 : else
2998 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2999 666063 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3000 561230 : return gerepileupto(av, FpX_to_mod(r, p));
3001 : }
3002 :
3003 : static GEN
3004 616 : zero_FpXQX_mod(GEN pol, GEN p, long v)
3005 : {
3006 616 : GEN r = cgetg(3,t_POL);
3007 616 : r[1] = evalvarn(v);
3008 616 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
3009 616 : return r;
3010 : }
3011 :
3012 : static GEN
3013 2079 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3014 : {
3015 2079 : pari_sp av = avma;
3016 : long dT;
3017 : GEN kx, ky, r;
3018 2079 : GEN T = RgX_to_FpX(pol, p);
3019 2079 : if (signe(T)==0) pari_err_OP("*", x, y);
3020 2072 : dT = degpol(T);
3021 2072 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3022 2072 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3023 2072 : r = FpX_mul(kx, ky, p);
3024 2072 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3025 1463 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3026 : }
3027 :
3028 : static GEN
3029 450787 : RgX_liftred(GEN x, GEN T)
3030 450787 : { return RgXQX_red(liftpol_shallow(x), T); }
3031 :
3032 : static GEN
3033 172246 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3034 : {
3035 172246 : pari_sp av = avma;
3036 172246 : long dT = degpol(T);
3037 172246 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3038 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3039 172246 : return gerepileupto(av, Kronecker_to_mod(r, T));
3040 : }
3041 :
3042 : static GEN
3043 1421 : RgX_sqr_FpX(GEN x, GEN p)
3044 : {
3045 1421 : pari_sp av = avma;
3046 : GEN r;
3047 1421 : if (lgefint(p) == 3)
3048 : {
3049 1217 : ulong pp = uel(p, 2);
3050 1217 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3051 : }
3052 : else
3053 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3054 1421 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3055 1281 : return gerepileupto(av, FpX_to_mod(r, p));
3056 : }
3057 :
3058 : static GEN
3059 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3060 : {
3061 196 : pari_sp av = avma;
3062 : long dT;
3063 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3064 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3065 189 : dT = degpol(T);
3066 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3067 189 : r = FpX_sqr(kx, p);
3068 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3069 189 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3070 : }
3071 :
3072 : static GEN
3073 13425 : RgX_sqr_QXQX(GEN x, GEN T)
3074 : {
3075 13425 : pari_sp av = avma;
3076 13425 : long dT = degpol(T);
3077 13425 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3078 13425 : return gerepileupto(av, Kronecker_to_mod(r, T));
3079 : }
3080 :
3081 : static GEN
3082 68320 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3083 : {
3084 68320 : pari_sp av = avma;
3085 : GEN r;
3086 68320 : if (lgefint(p) == 3)
3087 : {
3088 51684 : ulong pp = uel(p, 2);
3089 51684 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3090 : RgX_to_Flx(y, pp), pp));
3091 : }
3092 : else
3093 16636 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3094 68320 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3095 31024 : return gerepileupto(av, FpX_to_mod(r, p));
3096 : }
3097 :
3098 : static GEN
3099 46435 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3100 : {
3101 46435 : pari_sp av = avma;
3102 : GEN r;
3103 46435 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3104 46435 : return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
3105 : }
3106 : static GEN
3107 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3108 : {
3109 70 : pari_sp av = avma;
3110 : GEN r;
3111 70 : GEN T = RgX_to_FpX(pol, p);
3112 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3113 63 : if (lgefint(p) == 3)
3114 : {
3115 55 : ulong pp = uel(p, 2);
3116 55 : GEN Tp = ZX_to_Flx(T, pp);
3117 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3118 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3119 : }
3120 : else
3121 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3122 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3123 56 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3124 : }
3125 :
3126 : #define code(t1,t2) ((t1 << 6) | t2)
3127 : static GEN
3128 77519331 : RgX_mul_fast(GEN x, GEN y)
3129 : {
3130 : GEN p, pol;
3131 : long pa;
3132 77519331 : long t = RgX_type2(x,y, &p,&pol,&pa);
3133 77519596 : switch(t)
3134 : {
3135 41410420 : case t_INT: return ZX_mul(x,y);
3136 3342874 : case t_FRAC: return QX_mul(x,y);
3137 105176 : case t_FFELT: return FFX_mul(x, y, pol);
3138 666063 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3139 172246 : case code(t_POLMOD, t_INT):
3140 : case code(t_POLMOD, t_FRAC):
3141 172246 : return RgX_mul_QXQX(x, y, pol);
3142 2039 : case code(t_POLMOD, t_INTMOD):
3143 2039 : return RgX_mul_FpXQX(x, y, pol, p);
3144 31820778 : default: return NULL;
3145 : }
3146 : }
3147 : static GEN
3148 1344308 : RgX_sqr_fast(GEN x)
3149 : {
3150 : GEN p, pol;
3151 : long pa;
3152 1344308 : long t = RgX_type(x,&p,&pol,&pa);
3153 1344310 : switch(t)
3154 : {
3155 1239839 : case t_INT: return ZX_sqr(x);
3156 78034 : case t_FRAC: return QX_sqr(x);
3157 2499 : case t_FFELT: return FFX_sqr(x, pol);
3158 1421 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3159 13425 : case code(t_POLMOD, t_INT):
3160 : case code(t_POLMOD, t_FRAC):
3161 13425 : return RgX_sqr_QXQX(x, pol);
3162 195 : case code(t_POLMOD, t_INTMOD):
3163 195 : return RgX_sqr_FpXQX(x, pol, p);
3164 8897 : default: return NULL;
3165 : }
3166 : }
3167 :
3168 : static GEN
3169 14591447 : RgX_rem_fast(GEN x, GEN y)
3170 : {
3171 : GEN p, pol;
3172 : long pa;
3173 14591447 : long t = RgX_type2(x,y, &p,&pol,&pa);
3174 14591568 : switch(t)
3175 : {
3176 3402948 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3177 1814641 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3178 84 : case t_FFELT: return FFX_rem(x, y, pol);
3179 68320 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3180 46435 : case code(t_POLMOD, t_INT):
3181 : case code(t_POLMOD, t_FRAC):
3182 46435 : return RgX_rem_QXQX(x, y, pol);
3183 67 : case code(t_POLMOD, t_INTMOD):
3184 67 : return RgX_rem_FpXQX(x, y, pol, p);
3185 9259073 : default: return NULL;
3186 : }
3187 : }
3188 :
3189 : #undef code
3190 :
3191 : GEN
3192 63787156 : RgX_mul(GEN x, GEN y)
3193 : {
3194 63787156 : GEN z = RgX_mul_fast(x,y);
3195 63787286 : if (!z) z = RgX_mul_i(x,y);
3196 63786970 : return z;
3197 : }
3198 :
3199 : GEN
3200 1332156 : RgX_sqr(GEN x)
3201 : {
3202 1332156 : GEN z = RgX_sqr_fast(x);
3203 1332149 : if (!z) z = RgX_sqr_i(x);
3204 1332149 : return z;
3205 : }
3206 :
3207 : GEN
3208 14591448 : RgX_rem(GEN x, GEN y)
3209 : {
3210 14591448 : GEN z = RgX_rem_fast(x, y);
3211 14591546 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3212 14591509 : return z;
3213 : }
3214 :
3215 : static GEN
3216 0 : _RgX_mul(void* E, GEN x, GEN y)
3217 0 : { (void) E; return RgX_mul(x, y); }
3218 :
3219 : GEN
3220 0 : RgXV_prod(GEN V)
3221 0 : { return gen_product(V, NULL, &_RgX_mul); }
3222 :
3223 : GEN
3224 11061285 : RgXn_mul(GEN f, GEN g, long n)
3225 : {
3226 11061285 : pari_sp av = avma;
3227 11061285 : GEN h = RgX_mul_fast(f,g);
3228 11061286 : if (!h) return RgXn_mul2(f,g,n);
3229 1783470 : if (degpol(h) < n) return h;
3230 382941 : return gerepilecopy(av, RgXn_red_shallow(h, n));
3231 : }
3232 :
3233 : GEN
3234 12152 : RgXn_sqr(GEN f, long n)
3235 : {
3236 12152 : pari_sp av = avma;
3237 12152 : GEN g = RgX_sqr_fast(f);
3238 12152 : if (!g) return RgXn_sqr2(f,n);
3239 11837 : if (degpol(g) < n) return g;
3240 10640 : return gerepilecopy(av, RgXn_red_shallow(g, n));
3241 : }
3242 :
3243 : /* (f*g) \/ x^n */
3244 : GEN
3245 2670895 : RgX_mulhigh_i(GEN f, GEN g, long n)
3246 : {
3247 2670895 : GEN h = RgX_mul_fast(f,g);
3248 2670897 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3249 : }
3250 :
3251 : /* (f*g) \/ x^n */
3252 : GEN
3253 0 : RgX_sqrhigh_i(GEN f, long n)
3254 : {
3255 0 : GEN h = RgX_sqr_fast(f);
3256 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3257 : }
|
