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 19613826 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19613826 : long pold=1, rold=n*(d-1);
35 91839226 : for(p=2; p<=d; p++)
36 : {
37 72225400 : r = m*(p-1) + n*((d-1)/p);
38 72225400 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19613826 : return pold;
41 : }
42 :
43 : static GEN
44 14041062 : 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 14041062 : pari_sp av = avma;
48 : long i;
49 14041062 : GEN z = cmul(E,P,a,ff->one(E));
50 14025074 : if (!z) z = gen_0;
51 74628623 : for (i=1; i<=n; i++)
52 : {
53 60605748 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 60603512 : if (t) {
55 44638343 : z = ff->add(E, z, t);
56 44617184 : if (gc_needed(av,2)) z = gc_upto(av, z);
57 : }
58 : }
59 14022875 : 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 12123836 : 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 12123836 : pari_sp av = avma;
72 12123836 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12123836 : if (d < 0) return ff->zero(E);
76 11538243 : if (d < l) return gc_upto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1005423 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1005423 : 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 1005423 : d -= l;
84 1005423 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2502082 : while (d >= l-1)
86 : {
87 1495762 : d -= l-1;
88 1495762 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1495803 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1495753 : if (gc_needed(av,2))
91 91 : z = gc_upto(av, z);
92 : }
93 1006320 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1006321 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1006339 : return gc_upto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 881584 : 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 881584 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 881584 : if (d < 0) return ff->zero(E);
106 880310 : rtd = (long) sqrt((double)d);
107 880310 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 880319 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 880323 : return gc_upto(av, z);
110 : }
111 :
112 : static GEN
113 2151988 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19681222 : _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 1975812 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 635808 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2194583 : _gen_one(void *E) { (void)E; return gen_1; }
124 : static GEN
125 301 : _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 142885 : RgX_RgV_eval(GEN Q, GEN x)
136 : {
137 142885 : 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 395408 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 395408 : pari_sp av = avma, btop;
163 395408 : long i, d = degpol(P), o;
164 : GEN s;
165 395406 : if (signe(P)==0) return pol_0(varn(P));
166 395406 : s = gel(P, d+2);
167 395406 : if (d == 0) return gcopy(s);
168 392088 : o = RgX_deflate_order(P);
169 392096 : if (o > 1) A = gpowgs(A, o);
170 392099 : btop = avma;
171 1503297 : for (i = d-o; i >= 0; i-=o)
172 : {
173 1111201 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 1111198 : if (gc_needed(btop,1))
175 : {
176 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
177 13 : s = gc_upto(btop, s);
178 : }
179 : }
180 392096 : return gc_upto(av, s);
181 : }
182 :
183 : GEN
184 1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
185 : {
186 1652 : pari_sp av = avma;
187 1652 : long i, d = degpol(P), v = varn(A);
188 : GEN s;
189 1652 : if (signe(P)==0) return pol_0(v);
190 1652 : if (d == 0) return scalarpol(gel(P, d+2), v);
191 1232 : s = scalarpol_shallow(gel(P, d+2), v);
192 4963 : for (i = d-1; i >= 0; i--)
193 : {
194 3731 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
195 3731 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
196 3731 : if (gc_needed(av,1))
197 : {
198 0 : if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
199 0 : s = gc_upto(av, s);
200 : }
201 : }
202 1232 : return gc_upto(av, s);
203 : }
204 :
205 : const struct bb_algebra *
206 295021 : get_Rg_algebra(void)
207 : {
208 295021 : return &Rg_algebra;
209 : }
210 :
211 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
212 :
213 : static GEN
214 11809 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
215 : {
216 : (void) E;
217 11809 : return RgX_divrem(x, y, r);
218 : }
219 :
220 : GEN
221 3234 : RgX_digits(GEN x, GEN T)
222 : {
223 3234 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
224 3234 : if (signe(x)==0) return(cgetg(1, t_VEC));
225 3234 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
226 : }
227 :
228 : /*******************************************************************/
229 : /* */
230 : /* RgX */
231 : /* */
232 : /*******************************************************************/
233 :
234 : long
235 24795880 : RgX_equal(GEN x, GEN y)
236 : {
237 24795880 : long i = lg(x);
238 :
239 24795880 : if (i != lg(y)) return 0;
240 106699611 : for (i--; i > 1; i--)
241 82231755 : if (!gequal(gel(x,i),gel(y,i))) return 0;
242 24467856 : return 1;
243 : }
244 :
245 : /* Returns 1 in the base ring over which x is defined */
246 : /* HACK: this also works for t_SER */
247 : GEN
248 156745060 : Rg_get_1(GEN x)
249 : {
250 : GEN p, T;
251 156745060 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
252 156745062 : if (RgX_type_is_composite(tx))
253 11194785 : RgX_type_decode(tx, &i /*junk*/, &tx);
254 156745062 : switch(tx)
255 : {
256 421932 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
257 1015 : case t_PADIC: return cvtop(gen_1, p, lx);
258 6076 : case t_FFELT: return FF_1(T);
259 156316039 : default: return gen_1;
260 : }
261 : }
262 : /* Returns 0 in the base ring over which x is defined */
263 : /* HACK: this also works for t_SER */
264 : GEN
265 1862498 : Rg_get_0(GEN x)
266 : {
267 : GEN p, T;
268 1862498 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
269 1862498 : if (RgX_type_is_composite(tx))
270 25053 : RgX_type_decode(tx, &i /*junk*/, &tx);
271 1862498 : switch(tx)
272 : {
273 504 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
274 42 : case t_PADIC: return zeropadic(p, lx);
275 210 : case t_FFELT: return FF_zero(T);
276 1861742 : default: return gen_0;
277 : }
278 : }
279 :
280 : GEN
281 7455 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
282 : {
283 7455 : long i, n = degpol(P);
284 : GEN z, dz, dP;
285 7455 : if (n < 0) return gen_0;
286 7455 : P = Q_remove_denom(P, &dP);
287 7455 : z = gel(P,2); if (n == 0) return icopy(z);
288 4214 : if (dV) z = mulii(dV, z); /* V[1] = dV */
289 4214 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
290 7931 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
291 4214 : dz = mul_denom(dP, dV);
292 4214 : return dz? RgX_Rg_div(z, dz): z;
293 : }
294 :
295 : /* Return P(h * x), not memory clean */
296 : GEN
297 52066 : RgX_unscale(GEN P, GEN h)
298 : {
299 52066 : long i, l = lg(P);
300 52066 : GEN hi = gen_1, Q = cgetg(l, t_POL);
301 52067 : Q[1] = P[1];
302 52067 : if (l == 2) return Q;
303 40825 : gel(Q,2) = gcopy(gel(P,2));
304 95266 : for (i=3; i<l; i++)
305 : {
306 54441 : hi = gmul(hi,h);
307 54439 : gel(Q,i) = gmul(gel(P,i), hi);
308 : }
309 40825 : return Q;
310 : }
311 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 and 2^k */
312 : GEN
313 1172573 : ZX_z_unscale(GEN P, long h)
314 : {
315 1172573 : long i, l = lg(P);
316 1172573 : GEN Q = cgetg(l, t_POL);
317 1172573 : Q[1] = P[1];
318 1172573 : if (l == 2) return Q;
319 1129628 : gel(Q,2) = gel(P,2);
320 1129628 : if (l == 3) return Q;
321 1109034 : if (h == -1)
322 242049 : for (i = 3; i < l; i++)
323 : {
324 199263 : gel(Q,i) = negi(gel(P,i));
325 199264 : if (++i == l) break;
326 148232 : gel(Q,i) = gel(P,i);
327 : }
328 1015217 : else if (h > 0 && !(h & (h-1))) return ZX_unscale2n(P, vals(h));
329 : else
330 : {
331 : GEN hi;
332 493179 : gel(Q,3) = mulis(gel(P,3), h);
333 493179 : hi = sqrs(h);
334 2458823 : for (i = 4; i < l; i++)
335 : {
336 1965645 : gel(Q,i) = mulii(gel(P,i), hi);
337 1965644 : if (i != l-1) hi = mulis(hi,h);
338 : }
339 : }
340 586996 : return Q;
341 : }
342 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean */
343 : GEN
344 674724 : ZX_unscale(GEN P, GEN h)
345 : {
346 : long i, l;
347 : GEN Q, hi;
348 674724 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
349 888 : l = lg(P); Q = cgetg(l, t_POL);
350 888 : Q[1] = P[1];
351 888 : if (l == 2) return Q;
352 888 : gel(Q,2) = gel(P,2);
353 888 : if (l == 3) return Q;
354 888 : hi = h;
355 888 : gel(Q,3) = mulii(gel(P,3), hi);
356 2783 : for (i = 4; i < l; i++)
357 : {
358 1895 : hi = mulii(hi,h);
359 1895 : gel(Q,i) = mulii(gel(P,i), hi);
360 : }
361 888 : return Q;
362 : }
363 : /* P a ZX. Return P(x << n), not memory clean */
364 : GEN
365 1975089 : ZX_unscale2n(GEN P, long n)
366 : {
367 1975089 : long i, ni = n, l = lg(P);
368 1975089 : GEN Q = cgetg(l, t_POL);
369 1975088 : Q[1] = P[1];
370 1975088 : if (l == 2) return Q;
371 1929291 : gel(Q,2) = gel(P,2);
372 1929291 : if (l == 3) return Q;
373 1923537 : gel(Q,3) = shifti(gel(P,3), ni);
374 6725562 : for (i=4; i<l; i++)
375 : {
376 4802065 : ni += n;
377 4802065 : gel(Q,i) = shifti(gel(P,i), ni);
378 : }
379 1923497 : return Q;
380 : }
381 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
382 : GEN
383 12997 : ZX_unscale_div(GEN P, GEN h)
384 : {
385 12997 : long i, l = lg(P);
386 12997 : GEN hi, Q = cgetg(l, t_POL);
387 12997 : Q[1] = P[1];
388 12997 : if (l == 2) return Q;
389 12997 : gel(Q,2) = diviiexact(gel(P,2), h);
390 12997 : if (l == 3) return Q;
391 12997 : gel(Q,3) = gel(P,3);
392 12997 : if (l == 4) return Q;
393 12997 : hi = h;
394 12997 : gel(Q,4) = mulii(gel(P,4), hi);
395 64121 : for (i=5; i<l; i++)
396 : {
397 51124 : hi = mulii(hi,h);
398 51124 : gel(Q,i) = mulii(gel(P,i), hi);
399 : }
400 12997 : return Q;
401 : }
402 : /* P(h*X) / h^k, assuming the result is a ZX */
403 : GEN
404 1393 : ZX_unscale_divpow(GEN P, GEN h, long k)
405 : {
406 1393 : long i, j, l = lg(P);
407 1393 : GEN H, Q = cgetg(l, t_POL);
408 1393 : Q[1] = P[1]; if (l == 2) return Q;
409 1393 : H = gpowers(h, maxss(k, l - 3 - k));
410 5572 : for (i = 2, j = k+1; j > 1 && i < l; i++)
411 4179 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
412 1393 : if (i == l) return Q;
413 1393 : gel(Q, i) = gel(P, i); i++;
414 5082 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
415 1393 : return Q;
416 : }
417 :
418 : GEN
419 6488 : RgXV_unscale(GEN x, GEN h)
420 : {
421 6488 : if (isint1(h)) return gcopy(x);
422 18429 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
423 : }
424 :
425 : /* Return h^degpol(P) P(x / h), not memory clean */
426 : GEN
427 4345961 : RgX_rescale(GEN P, GEN h)
428 : {
429 4345961 : long i, l = lg(P);
430 4345961 : GEN Q = cgetg(l,t_POL), hi = h;
431 4345991 : gel(Q,l-1) = gel(P,l-1);
432 11437093 : for (i=l-2; i>=2; i--)
433 : {
434 11434529 : gel(Q,i) = gmul(gel(P,i), hi);
435 11434425 : if (i == 2) break;
436 7090871 : hi = gmul(hi,h);
437 : }
438 4346118 : Q[1] = P[1]; return Q;
439 : }
440 :
441 : GEN
442 2401 : RgXV_rescale(GEN x, GEN h)
443 : {
444 2401 : if (isint1(h)) return RgX_copy(x);
445 16086 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
446 : }
447 :
448 : /* A(X^d) --> A(X) */
449 : GEN
450 1186162 : RgX_deflate(GEN x0, long d)
451 : {
452 : GEN z, y, x;
453 1186162 : long i,id, dy, dx = degpol(x0);
454 1186180 : if (d == 1 || dx <= 0) return x0;
455 462546 : dy = dx/d;
456 462546 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
457 462541 : z = y + 2;
458 462541 : x = x0+ 2;
459 1772695 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
460 462541 : return y;
461 : }
462 :
463 : GEN
464 1260 : RgX_homogenize_deg(GEN P, long d, long v)
465 : {
466 : long i, l;
467 1260 : GEN Q = cgetg_copy(P, &l);
468 1260 : Q[1] = P[1];
469 4018 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
470 1260 : return Q;
471 : }
472 :
473 : GEN
474 18102 : RgX_homogenize(GEN P, long v)
475 : {
476 : long i, l, d;
477 18102 : GEN Q = cgetg_copy(P, &l);
478 18102 : Q[1] = P[1]; d = l-3;
479 159348 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
480 18102 : return Q;
481 : }
482 :
483 : /* F a t_RFRAC */
484 : long
485 140 : rfrac_deflate_order(GEN F)
486 : {
487 140 : GEN N = gel(F,1), D = gel(F,2);
488 140 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
489 140 : if (m == 1) return 1;
490 49 : if (typ(N) == t_POL && varn(N) == varn(D))
491 28 : m = cgcd(m, RgX_deflate_order(N));
492 49 : return m;
493 : }
494 : /* F a t_RFRAC */
495 : GEN
496 140 : rfrac_deflate_max(GEN F, long *m)
497 : {
498 140 : *m = rfrac_deflate_order(F);
499 140 : return rfrac_deflate(F, *m);
500 : }
501 : /* F a t_RFRAC */
502 : GEN
503 140 : rfrac_deflate(GEN F, long m)
504 : {
505 140 : GEN N = gel(F,1), D = gel(F,2);
506 140 : if (m == 1) return F;
507 49 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
508 49 : D = RgX_deflate(D, m); return mkrfrac(N, D);
509 : }
510 :
511 : /* return x0(X^d) */
512 : GEN
513 884763 : RgX_inflate(GEN x0, long d)
514 : {
515 884763 : long i, id, dy, dx = degpol(x0);
516 884763 : GEN x = x0 + 2, z, y;
517 884763 : if (dx <= 0) return leafcopy(x0);
518 807473 : dy = dx*d;
519 807473 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
520 807473 : z = y + 2;
521 28101801 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
522 11499711 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
523 807473 : return y;
524 : }
525 :
526 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
527 : static GEN
528 5709627 : RgX_Rg_translate_basecase(GEN P, GEN c)
529 : {
530 5709627 : pari_sp av = avma;
531 : GEN Q;
532 : long i, k, n;
533 :
534 5709627 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
535 5707696 : Q = leafcopy(P); n = degpol(P);
536 5707694 : if (isint1(c))
537 : {
538 6867 : for (i=1; i<=n; i++)
539 : {
540 20265 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gel(Q,2+k+1));
541 5390 : if (gc_needed(av,2))
542 : {
543 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(1), i = %ld/%ld", i,n);
544 0 : Q = gc_GEN(av, Q);
545 : }
546 : }
547 : }
548 5706214 : else if (isintm1(c))
549 : {
550 15246 : for (i=1; i<=n; i++)
551 : {
552 49630 : for (k=n-i; k<n; k++) gel(Q,2+k) = gsub(gel(Q,2+k), gel(Q,2+k+1));
553 12299 : if (gc_needed(av,2))
554 : {
555 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(-1), i = %ld/%ld", i,n);
556 0 : Q = gc_GEN(av, Q);
557 : }
558 : }
559 : }
560 : else
561 : {
562 19453925 : for (i=1; i<=n; i++)
563 : {
564 45505375 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1)));
565 13750657 : if (gc_needed(av,2))
566 : {
567 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate, i = %ld/%ld", i,n);
568 0 : Q = gc_GEN(av, Q);
569 : }
570 : }
571 : }
572 5707477 : return gc_GEN(av, Q);
573 : }
574 :
575 : static GEN
576 82656 : zero_FpX_mod(GEN p, long v)
577 : {
578 82656 : GEN r = cgetg(3,t_POL);
579 82656 : r[1] = evalvarn(v);
580 82656 : gel(r,2) = mkintmod(gen_0, icopy(p));
581 82656 : return r;
582 : }
583 :
584 : static GEN
585 91 : zero_FpXQX_mod(GEN pol, GEN p, long v)
586 : {
587 91 : GEN r = cgetg(3,t_POL);
588 91 : r[1] = evalvarn(v);
589 91 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
590 91 : return r;
591 : }
592 :
593 : static GEN
594 0 : RgX_Rg_translate_FpX(GEN P, GEN c, GEN p)
595 : {
596 0 : pari_sp av = avma;
597 : GEN r;
598 : #if 0
599 : /* 'divide by 0' error if p is not prime and c not invertible */
600 : if (lgefint(p) == 3)
601 : {
602 : ulong pp = uel(p, 2);
603 : r = Flx_to_ZX_inplace(Flx_Fl_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
604 : }
605 : else
606 : #endif
607 0 : r = FpX_Fp_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
608 0 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
609 0 : return gc_upto(av, FpX_to_mod(r, p));
610 : }
611 :
612 : static GEN
613 7 : RgX_Rg_translate_FpXQX(GEN x, GEN c, GEN pol, GEN p)
614 : {
615 7 : pari_sp av = avma;
616 7 : GEN r, T = RgX_to_FpX(pol, p);
617 7 : if (signe(T)==0) pari_err_OP("subst", x, c);
618 7 : r = FpXQX_FpXQ_translate(RgX_to_FpXQX(x, T, p), Rg_to_FpXQ(c, T, p), T, p);
619 7 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
620 7 : return gc_upto(av, FpXQX_to_mod(r, T, p));
621 : }
622 :
623 : static GEN
624 6070292 : RgX_Rg_translate_fast(GEN P, GEN c)
625 : {
626 : GEN p, pol;
627 : long pa;
628 6070292 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
629 6070322 : switch(t)
630 : {
631 360875 : case t_INT: return ZX_Z_translate(P, c);
632 0 : case t_INTMOD: return RgX_Rg_translate_FpX(P, c, p);
633 6 : case RgX_type_code(t_POLMOD, t_INTMOD):
634 6 : return RgX_Rg_translate_FpXQX(P, c, pol, p);
635 5709441 : default: return NULL;
636 : }
637 : }
638 :
639 : static GEN
640 5709819 : RgX_Rg_translate_i(GEN P, GEN c)
641 : {
642 5709819 : pari_sp av = avma;
643 : long n;
644 5709819 : n = degpol(P);
645 5709817 : if (n < 40)
646 5709628 : return RgX_Rg_translate_basecase(P, c);
647 : else
648 : {
649 189 : long d = n >> 1;
650 189 : GEN Q = RgX_Rg_translate_i(RgX_shift_shallow(P, -d), c);
651 189 : GEN R = RgX_Rg_translate_i(RgXn_red_shallow(P, d), c);
652 189 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
653 189 : return gc_upto(av, RgX_add(RgX_mul(Q, S), R));
654 : }
655 : }
656 :
657 : GEN
658 6070326 : RgX_Rg_translate(GEN P, GEN c)
659 : {
660 6070326 : GEN R = RgX_Rg_translate_fast(P, c);
661 6070325 : return R ? R: RgX_Rg_translate_i(P,c);
662 : }
663 : /* P(ax + b) */
664 : GEN
665 30226 : RgX_affine(GEN P, GEN a, GEN b)
666 : {
667 30226 : if (!gequal0(b)) P = RgX_Rg_translate(P, b);
668 30226 : return RgX_unscale(P, a);
669 : }
670 :
671 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
672 : GEN
673 33584 : RgXQX_RgXQ_translate(GEN P, GEN c, GEN T)
674 : {
675 33584 : pari_sp av = avma;
676 : GEN Q;
677 : long i, k, n;
678 :
679 33584 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
680 33241 : Q = leafcopy(P); n = degpol(P);
681 105778 : for (i=1; i<=n; i++)
682 : {
683 303686 : for (k=n-i; k<n; k++)
684 : {
685 231149 : pari_sp av2 = avma;
686 231149 : gel(Q,2+k) = gc_upto(av2,
687 231149 : RgX_rem(gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1))), T));
688 : }
689 72537 : if (gc_needed(av,2))
690 : {
691 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_RgXQ_translate, i = %ld/%ld", i,n);
692 0 : Q = gc_GEN(av, Q);
693 : }
694 : }
695 33241 : return gc_GEN(av, Q);
696 : }
697 :
698 : /********************************************************************/
699 : /** **/
700 : /** CONVERSIONS **/
701 : /** (not memory clean) **/
702 : /** **/
703 : /********************************************************************/
704 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
705 : * but everything else is */
706 : static GEN
707 168329 : QXQ_to_mod(GEN x, GEN T)
708 : {
709 : long d;
710 168329 : switch(typ(x))
711 : {
712 68391 : case t_INT: return icopy(x);
713 2079 : case t_FRAC: return gcopy(x);
714 97859 : case t_POL:
715 97859 : d = degpol(x);
716 97859 : if (d < 0) return gen_0;
717 97831 : if (d == 0) return gcopy(gel(x,2));
718 97761 : return mkpolmod(RgX_copy(x), T);
719 0 : default: pari_err_TYPE("QXQ_to_mod",x);
720 : return NULL;/* LCOV_EXCL_LINE */
721 : }
722 : }
723 : /* pure shallow version */
724 : GEN
725 840695 : QXQ_to_mod_shallow(GEN x, GEN T)
726 : {
727 : long d;
728 840695 : switch(typ(x))
729 : {
730 546003 : case t_INT:
731 546003 : case t_FRAC: return x;
732 294692 : case t_POL:
733 294692 : d = degpol(x);
734 294692 : if (d < 0) return gen_0;
735 249702 : if (d == 0) return gel(x,2);
736 233361 : return mkpolmod(x, T);
737 0 : default: pari_err_TYPE("QXQ_to_mod",x);
738 : return NULL;/* LCOV_EXCL_LINE */
739 : }
740 : }
741 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
742 : * Not memory clean because T not copied, but everything else is */
743 : static GEN
744 37212 : QXQX_to_mod(GEN z, GEN T)
745 : {
746 37212 : long i,l = lg(z);
747 37212 : GEN x = cgetg(l,t_POL);
748 186046 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
749 37212 : x[1] = z[1]; return normalizepol_lg(x,l);
750 : }
751 : /* pure shallow version */
752 : GEN
753 204096 : QXQX_to_mod_shallow(GEN z, GEN T)
754 : {
755 204096 : long i,l = lg(z);
756 204096 : GEN x = cgetg(l,t_POL);
757 961337 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
758 204096 : x[1] = z[1]; return normalizepol_lg(x,l);
759 : }
760 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
761 : GEN
762 12572 : QXQXV_to_mod(GEN V, GEN T)
763 : {
764 12572 : long i, l = lg(V);
765 12572 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
766 49784 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
767 12572 : return z;
768 : }
769 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
770 : GEN
771 22339 : QXQV_to_mod(GEN V, GEN T)
772 : {
773 22339 : long i, l = lg(V);
774 22339 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
775 41834 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
776 22339 : return z;
777 : }
778 :
779 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
780 : GEN
781 14854 : QXQC_to_mod_shallow(GEN V, GEN T)
782 : {
783 14854 : long i, l = lg(V);
784 14854 : GEN z = cgetg(l, t_COL);
785 98308 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
786 14854 : return z;
787 : }
788 :
789 : GEN
790 6720 : QXQM_to_mod_shallow(GEN V, GEN T)
791 : {
792 6720 : long i, l = lg(V);
793 6720 : GEN z = cgetg(l, t_MAT);
794 21574 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
795 6720 : return z;
796 : }
797 :
798 : GEN
799 7740076 : RgX_renormalize_lg(GEN x, long lx)
800 : {
801 : long i;
802 10730931 : for (i = lx-1; i>1; i--)
803 10329969 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
804 7740077 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
805 7740078 : setlg(x, i+1); setsigne(x, i != 1); return x;
806 : }
807 :
808 : GEN
809 2280969 : RgV_to_RgX(GEN x, long v)
810 : {
811 2280969 : long i, k = lg(x);
812 : GEN p;
813 :
814 6352629 : while (--k && gequal0(gel(x,k)));
815 2280969 : if (!k) return pol_0(v);
816 2248100 : i = k+2; p = cgetg(i,t_POL);
817 2248101 : p[1] = evalsigne(1) | evalvarn(v);
818 24372287 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
819 2248101 : return p;
820 : }
821 : GEN
822 188283 : RgV_to_RgX_reverse(GEN x, long v)
823 : {
824 188283 : long j, k, l = lg(x);
825 : GEN p;
826 :
827 189851 : for (k = 1; k < l; k++)
828 189851 : if (!gequal0(gel(x,k))) break;
829 188283 : if (k == l) return pol_0(v);
830 188283 : k -= 1;
831 188283 : l -= k;
832 188283 : x += k;
833 188283 : p = cgetg(l+1,t_POL);
834 188283 : p[1] = evalsigne(1) | evalvarn(v);
835 985038 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
836 188283 : return p;
837 : }
838 :
839 : /* return the (N-dimensional) vector of coeffs of p */
840 : GEN
841 15053115 : RgX_to_RgC(GEN x, long N)
842 : {
843 : long i, l;
844 : GEN z;
845 15053115 : l = lg(x)-1; x++;
846 15053115 : if (l > N+1) l = N+1; /* truncate higher degree terms */
847 15053115 : z = cgetg(N+1,t_COL);
848 110958632 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
849 31113468 : for ( ; i<=N; i++) gel(z,i) = gen_0;
850 15053272 : return z;
851 : }
852 : GEN
853 1360424 : Rg_to_RgC(GEN x, long N)
854 : {
855 1360424 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
856 : }
857 :
858 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
859 : GEN
860 302964 : RgM_to_RgXV(GEN x, long v)
861 1289771 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
862 : GEN
863 7202 : RgM_to_RgXV_reverse(GEN x, long v)
864 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
865 :
866 : /* matrix whose entries are given by the coeffs of the polynomials in
867 : * vector v (considered as degree n-1 polynomials) */
868 : GEN
869 336566 : RgV_to_RgM(GEN x, long n)
870 1691788 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
871 :
872 : GEN
873 79010 : RgXV_to_RgM(GEN x, long n)
874 401769 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
875 :
876 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
877 : GEN
878 23714 : RgM_to_RgXX(GEN x, long v,long w)
879 : {
880 23714 : long j, lx = lg(x);
881 23714 : GEN y = cgetg(lx+1, t_POL);
882 23714 : y[1] = evalsigne(1) | evalvarn(v);
883 23714 : y++;
884 132085 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
885 23714 : return normalizepol_lg(--y, lx+1);
886 : }
887 :
888 : /* matrix whose entries are given by the coeffs of the polynomial v in
889 : * two variables (considered as degree n-1 polynomials) */
890 : GEN
891 322 : RgXX_to_RgM(GEN v, long n)
892 : {
893 322 : long j, N = lg(v)-1;
894 322 : GEN y = cgetg(N, t_MAT);
895 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
896 322 : return y;
897 : }
898 :
899 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
900 : GEN
901 30875 : RgXY_swapspec(GEN x, long n, long w, long nx)
902 : {
903 30875 : long j, ly = n+3;
904 30875 : GEN y = cgetg(ly, t_POL);
905 30875 : y[1] = evalsigne(1);
906 385383 : for (j=2; j<ly; j++)
907 : {
908 : long k;
909 354508 : GEN a = cgetg(nx+2,t_POL);
910 354508 : a[1] = evalsigne(1) | evalvarn(w);
911 1971941 : for (k=0; k<nx; k++)
912 : {
913 1617433 : GEN xk = gel(x,k);
914 1617433 : if (typ(xk)==t_POL && varn(xk)==w)
915 1522069 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
916 : else
917 95364 : gel(a,k+2) = j==2 ? xk: gen_0;
918 : }
919 354508 : gel(y,j) = normalizepol_lg(a, nx+2);
920 : }
921 30875 : return normalizepol_lg(y,ly);
922 : }
923 :
924 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
925 : GEN
926 2051 : RgXY_swap(GEN x, long n, long w)
927 : {
928 2051 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
929 2051 : setvarn(z, varn(x)); return z;
930 : }
931 :
932 : long
933 2022351 : RgXY_degreex(GEN b)
934 : {
935 2022351 : long deg = 0, i;
936 2022351 : if (!signe(b)) return -1;
937 8496973 : for (i = 2; i < lg(b); ++i)
938 : {
939 6474622 : GEN bi = gel(b, i);
940 6474622 : if (typ(bi) == t_POL)
941 1023339 : deg = maxss(deg, degpol(bi));
942 : }
943 2022351 : return deg;
944 : }
945 :
946 : GEN
947 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
948 :
949 : /* return (x % X^n). Shallow */
950 : GEN
951 7847363 : RgXn_red_shallow(GEN a, long n)
952 : {
953 7847363 : long i, L = n+2, l = lg(a);
954 : GEN b;
955 7847363 : if (L >= l) return a; /* deg(x) < n */
956 5652658 : b = cgetg(L, t_POL); b[1] = a[1];
957 37217260 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
958 5652659 : return normalizepol_lg(b,L);
959 : }
960 :
961 : GEN
962 483 : RgXnV_red_shallow(GEN x, long n)
963 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
964 :
965 : /* return (x * X^n). Shallow */
966 : GEN
967 175068056 : RgX_shift_shallow(GEN a, long n)
968 : {
969 175068056 : long i, l = lg(a);
970 : GEN b;
971 175068056 : if (l == 2 || !n) return a;
972 108836181 : l += n;
973 108836181 : if (n < 0)
974 : {
975 53804281 : if (l <= 2) return pol_0(varn(a));
976 52379378 : b = cgetg(l, t_POL); b[1] = a[1];
977 52379983 : a -= n;
978 163166057 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
979 : } else {
980 55031900 : b = cgetg(l, t_POL); b[1] = a[1];
981 55035431 : a -= n; n += 2;
982 119597705 : for (i=2; i<n; i++) gel(b,i) = gen_0;
983 212336762 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
984 : }
985 107415414 : return b;
986 : }
987 : /* return (x * X^n). */
988 : GEN
989 1577707 : RgX_shift(GEN a, long n)
990 : {
991 1577707 : long i, l = lg(a);
992 : GEN b;
993 1577707 : if (l == 2 || !n) return RgX_copy(a);
994 1576755 : l += n;
995 1576755 : if (n < 0)
996 : {
997 1442 : if (l <= 2) return pol_0(varn(a));
998 1372 : b = cgetg(l, t_POL); b[1] = a[1];
999 1372 : a -= n;
1000 9534 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1001 : } else {
1002 1575313 : b = cgetg(l, t_POL); b[1] = a[1];
1003 1575313 : a -= n; n += 2;
1004 3976311 : for (i=2; i<n; i++) gel(b,i) = gen_0;
1005 4343570 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1006 : }
1007 1576685 : return b;
1008 : }
1009 :
1010 : GEN
1011 314657 : RgX_rotate_shallow(GEN P, long k, long p)
1012 : {
1013 314657 : long i, l = lgpol(P);
1014 : GEN r;
1015 314657 : if (signe(P)==0)
1016 805 : return pol_0(varn(P));
1017 313852 : r = cgetg(p+2,t_POL); r[1] = P[1];
1018 2089724 : for(i=0; i<p; i++)
1019 : {
1020 1775872 : long s = 2+(i+k)%p;
1021 1775872 : gel(r,s) = i<l? gel(P,2+i): gen_0;
1022 : }
1023 313852 : return RgX_renormalize(r);
1024 : }
1025 :
1026 : GEN
1027 2997489 : RgX_mulXn(GEN x, long d)
1028 : {
1029 : pari_sp av;
1030 : GEN z;
1031 : long v;
1032 2997489 : if (d >= 0) return RgX_shift(x, d);
1033 1472762 : d = -d;
1034 1472762 : v = RgX_val(x);
1035 1472762 : if (v >= d) return RgX_shift(x, -d);
1036 1472748 : av = avma;
1037 1472748 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1038 1472748 : return gc_upto(av, z);
1039 : }
1040 :
1041 : long
1042 588 : RgXV_maxdegree(GEN x)
1043 : {
1044 588 : long d = -1, i, l = lg(x);
1045 4494 : for (i = 1; i < l; i++)
1046 3906 : d = maxss(d, degpol(gel(x,i)));
1047 588 : return d;
1048 : }
1049 :
1050 : long
1051 3667633 : RgX_val(GEN x)
1052 : {
1053 3667633 : long i, lx = lg(x);
1054 3667633 : if (lx == 2) return LONG_MAX;
1055 4631267 : for (i = 2; i < lx; i++)
1056 4631225 : if (!isexactzero(gel(x,i))) break;
1057 3667479 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1058 3667437 : return i - 2;
1059 : }
1060 : long
1061 80712911 : RgX_valrem(GEN x, GEN *Z)
1062 : {
1063 80712911 : long v, i, lx = lg(x);
1064 80712911 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1065 125068257 : for (i = 2; i < lx; i++)
1066 125068738 : if (!isexactzero(gel(x,i))) break;
1067 : /* possible with nonrational zeros */
1068 80713318 : if (i == lx)
1069 : {
1070 21 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1071 21 : return LONG_MAX;
1072 : }
1073 80713297 : v = i - 2;
1074 80713297 : *Z = RgX_shift_shallow(x, -v);
1075 80722065 : return v;
1076 : }
1077 : long
1078 852175 : RgX_valrem_inexact(GEN x, GEN *Z)
1079 : {
1080 : long v;
1081 852175 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1082 869162 : for (v = 0;; v++)
1083 869162 : if (!gequal0(gel(x,2+v))) break;
1084 852168 : if (Z) *Z = RgX_shift_shallow(x, -v);
1085 852168 : return v;
1086 : }
1087 :
1088 : GEN
1089 68187 : RgXQC_red(GEN x, GEN T)
1090 414505 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1091 :
1092 : GEN
1093 1449 : RgXQV_red(GEN x, GEN T)
1094 35007 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1095 :
1096 : GEN
1097 13195 : RgXQM_red(GEN x, GEN T)
1098 81382 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1099 :
1100 : GEN
1101 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1102 : {
1103 322 : return RgXQM_red(RgM_mul(P, Q), T);
1104 : }
1105 :
1106 : GEN
1107 497504 : RgXQX_red(GEN P, GEN T)
1108 : {
1109 497504 : long i, l = lg(P);
1110 497504 : GEN Q = cgetg(l, t_POL);
1111 497504 : Q[1] = P[1];
1112 2590699 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1113 497504 : return normalizepol_lg(Q, l);
1114 : }
1115 :
1116 : GEN
1117 826852 : RgX_deriv(GEN x)
1118 : {
1119 826852 : long i,lx = lg(x)-1;
1120 : GEN y;
1121 :
1122 826852 : if (lx<3) return pol_0(varn(x));
1123 823870 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1124 3647297 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1125 823864 : y[1] = x[1]; return normalizepol_lg(y,i);
1126 : }
1127 :
1128 : GEN
1129 2583276 : RgX_recipspec_shallow(GEN x, long l, long n)
1130 : {
1131 : long i;
1132 2583276 : GEN z = cgetg(n+2,t_POL);
1133 2583277 : z[1] = 0; z += 2;
1134 144288138 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1135 2828800 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1136 2583277 : return normalizepol_lg(z-2,n+2);
1137 : }
1138 :
1139 : GEN
1140 642591 : RgXn_recip_shallow(GEN P, long n)
1141 : {
1142 642591 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1143 642595 : setvarn(Q, varn(P));
1144 642595 : return Q;
1145 : }
1146 :
1147 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1148 : GEN
1149 30541 : RgX_recip(GEN x)
1150 : {
1151 : long lx, i, j;
1152 30541 : GEN y = cgetg_copy(x, &lx);
1153 267519 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1154 30541 : return normalizepol_lg(y,lx);
1155 : }
1156 : /* shallow version */
1157 : GEN
1158 59388 : RgX_recip_shallow(GEN x)
1159 : {
1160 : long lx, i, j;
1161 59388 : GEN y = cgetg_copy(x, &lx);
1162 356118 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1163 59388 : return normalizepol_lg(y,lx);
1164 : }
1165 :
1166 : GEN
1167 3751325 : RgX_recip_i(GEN x)
1168 : {
1169 : long lx, i, j;
1170 3751325 : GEN y = cgetg_copy(x, &lx);
1171 20232684 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1172 3751324 : return y;
1173 : }
1174 : /*******************************************************************/
1175 : /* */
1176 : /* ADDITION / SUBTRACTION */
1177 : /* */
1178 : /*******************************************************************/
1179 : /* cf RgX_coeff_simplify */
1180 : INLINE GEN
1181 440810798 : RgX_coeff_add(GEN a, GEN b)
1182 : {
1183 440810798 : pari_sp av = avma;
1184 440810798 : GEN s = gadd(a, b);
1185 440777313 : if (typ(s) != t_POL) return s;
1186 50951076 : switch(lg(s))
1187 : {
1188 1505140 : case 2: return gc_const(av, gen_0);
1189 131895 : case 3: return gc_GEN(av, gel(s,2));
1190 49314041 : default: return s;
1191 : }
1192 : }
1193 : INLINE GEN
1194 246193496 : RgX_coeff_sub(GEN a, GEN b)
1195 : {
1196 246193496 : pari_sp av = avma;
1197 246193496 : GEN s = gsub(a, b);
1198 246093385 : if (typ(s) != t_POL) return s;
1199 9748550 : switch(lg(s))
1200 : {
1201 2275077 : case 2: return gc_const(av, gen_0);
1202 249840 : case 3: return gc_GEN(av, gel(s,2));
1203 7223633 : default: return s;
1204 : }
1205 : }
1206 : /* same variable */
1207 : GEN
1208 117232706 : RgX_add(GEN x, GEN y)
1209 : {
1210 117232706 : long i, lx = lg(x), ly = lg(y);
1211 : GEN z;
1212 117232706 : if (ly <= lx) {
1213 104359641 : z = cgetg(lx,t_POL); z[1] = x[1];
1214 480372089 : for (i=2; i < ly; i++) gel(z,i) = RgX_coeff_add(gel(x,i),gel(y,i));
1215 160865054 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1216 104340513 : z = normalizepol_lg(z, lx);
1217 : } else {
1218 12873065 : z = cgetg(ly,t_POL); z[1] = y[1];
1219 52029703 : for (i=2; i < lx; i++) gel(z,i) = RgX_coeff_add(gel(x,i),gel(y,i));
1220 37488422 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1221 12873967 : z = normalizepol_lg(z, ly);
1222 : }
1223 117226440 : return z;
1224 : }
1225 : GEN
1226 68732616 : RgX_sub(GEN x, GEN y)
1227 : {
1228 68732616 : long i, lx = lg(x), ly = lg(y);
1229 : GEN z;
1230 68732616 : if (ly <= lx) {
1231 35022732 : z = cgetg(lx,t_POL); z[1] = x[1];
1232 184267005 : for (i=2; i < ly; i++) gel(z,i) = RgX_coeff_sub(gel(x,i),gel(y,i));
1233 61024313 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1234 35003504 : z = normalizepol_lg(z, lx);
1235 : } else {
1236 33709884 : z = cgetg(ly,t_POL); z[1] = y[1];
1237 125622038 : for (i=2; i < lx; i++) gel(z,i) = RgX_coeff_sub(gel(x,i),gel(y,i));
1238 74305717 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1239 33682938 : z = normalizepol_lg(z, ly);
1240 : }
1241 68728739 : return z;
1242 : }
1243 : GEN
1244 6490469 : RgX_neg(GEN x)
1245 46158264 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1246 :
1247 : GEN
1248 26288027 : RgX_Rg_add(GEN y, GEN x)
1249 : {
1250 : GEN z;
1251 26288027 : long lz = lg(y), i;
1252 26288027 : if (lz == 2) return scalarpol(x,varn(y));
1253 25572692 : z = cgetg(lz,t_POL); z[1] = y[1];
1254 25572664 : gel(z,2) = RgX_coeff_add(gel(y,2),x);
1255 86268978 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1256 : /* probably useless unless lz = 3, but cannot be skipped if y is
1257 : * an inexact 0 */
1258 25572569 : return normalizepol_lg(z,lz);
1259 : }
1260 : GEN
1261 65078 : RgX_Rg_add_shallow(GEN y, GEN x)
1262 : {
1263 : GEN z;
1264 65078 : long lz = lg(y), i;
1265 65078 : if (lz == 2) return scalarpol(x,varn(y));
1266 65078 : z = cgetg(lz,t_POL); z[1] = y[1];
1267 65078 : gel(z,2) = RgX_coeff_add(gel(y,2),x);
1268 130292 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1269 65078 : return normalizepol_lg(z,lz);
1270 : }
1271 : GEN
1272 223071 : RgX_Rg_sub(GEN y, GEN x)
1273 : {
1274 : GEN z;
1275 223071 : long lz = lg(y), i;
1276 223071 : if (lz == 2)
1277 : { /* scalarpol(gneg(x),varn(y)) optimized */
1278 133 : long v = varn(y);
1279 133 : if (isrationalzero(x)) return pol_0(v);
1280 7 : z = cgetg(3,t_POL);
1281 7 : z[1] = gequal0(x)? evalvarn(v)
1282 7 : : evalvarn(v) | evalsigne(1);
1283 7 : gel(z,2) = gneg(x); return z;
1284 : }
1285 222938 : z = cgetg(lz,t_POL); z[1] = y[1];
1286 222937 : gel(z,2) = RgX_coeff_sub(gel(y,2),x);
1287 589422 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1288 222937 : return normalizepol_lg(z,lz);
1289 : }
1290 : GEN
1291 4880158 : Rg_RgX_sub(GEN x, GEN y)
1292 : {
1293 : GEN z;
1294 4880158 : long lz = lg(y), i;
1295 4880158 : if (lz == 2) return scalarpol(x,varn(y));
1296 4880018 : z = cgetg(lz,t_POL); z[1] = y[1];
1297 4879999 : gel(z,2) = RgX_coeff_sub(x, gel(y,2));
1298 7330498 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1299 4879988 : return normalizepol_lg(z,lz);
1300 : }
1301 : /*******************************************************************/
1302 : /* */
1303 : /* KARATSUBA MULTIPLICATION */
1304 : /* */
1305 : /*******************************************************************/
1306 : #if 0
1307 : /* to debug Karatsuba-like routines */
1308 : GEN
1309 : zx_debug_spec(GEN x, long nx)
1310 : {
1311 : GEN z = cgetg(nx+2,t_POL);
1312 : long i;
1313 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1314 : z[1] = evalsigne(1); return z;
1315 : }
1316 :
1317 : GEN
1318 : RgX_debug_spec(GEN x, long nx)
1319 : {
1320 : GEN z = cgetg(nx+2,t_POL);
1321 : long i;
1322 : for (i=0; i<nx; i++) z[i+2] = x[i];
1323 : z[1] = evalsigne(1); return z;
1324 : }
1325 : #endif
1326 :
1327 : /* generic multiplication */
1328 : GEN
1329 8856939 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1330 : {
1331 : GEN z, t;
1332 : long i;
1333 8856939 : if (nx == ny) {
1334 1480615 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1335 4799995 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1336 1480608 : return normalizepol_lg(z, nx+2);
1337 : }
1338 7376324 : if (ny < nx) {
1339 7190973 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1340 26463940 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1341 17181942 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1342 7190607 : return normalizepol_lg(z, nx+2);
1343 : } else {
1344 185351 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1345 3672720 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1346 447626 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1347 185364 : return normalizepol_lg(z, ny+2);
1348 : }
1349 : }
1350 : GEN
1351 222462 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1352 : {
1353 : GEN z, t;
1354 : long i;
1355 222462 : if (nx == ny) {
1356 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1357 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1358 12824 : return normalizepol_lg(z, nx+2);
1359 : }
1360 209638 : if (ny < nx) {
1361 207853 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1362 3725563 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1363 2373002 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1364 207853 : return normalizepol_lg(z, nx+2);
1365 : } else {
1366 1785 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1367 331478 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1368 12236 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1369 1785 : return normalizepol_lg(z, ny+2);
1370 : }
1371 : }
1372 :
1373 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1374 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1375 : * t_VECSMALL, to hold this, which can be left on stack: GC functions
1376 : * will not crash on it. The returned vector itself is not a proper GEN,
1377 : * we access the coefficients as x[i], i = 0..deg(x) */
1378 : static GEN
1379 79643936 : RgXspec_kill0(GEN x, long lx)
1380 : {
1381 79643936 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit GC-wise */
1382 : long i;
1383 242190673 : for (i=0; i <lx; i++)
1384 : {
1385 162546717 : GEN c = gel(x,i);
1386 162546717 : z[i] = (long)(isrationalzero(c)? NULL: c);
1387 : }
1388 79643956 : return z;
1389 : }
1390 :
1391 : /* For consistency we want operations on t_POL of the same variable v
1392 : * to return a t_POL in v; but this becomes a nuisance for coefficients,
1393 : * e.g., multivariate t_POL, which should have simplest possible types */
1394 : INLINE GEN
1395 109634420 : RgX_coeff_simplify(pari_sp av, GEN s)
1396 : {
1397 109634420 : if (!s) return gen_0;
1398 102584371 : if (typ(s) != t_POL) return gc_upto(av, s);
1399 24892597 : switch(lg(s))
1400 : {
1401 12222 : case 2: return gen_0;
1402 2151997 : case 3: return gc_GEN(av, gel(s,2));
1403 22728378 : default: return gc_upto(av, s);
1404 : }
1405 : }
1406 : INLINE GEN
1407 109543364 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1408 : {
1409 109543364 : pari_sp av = avma;
1410 109543364 : GEN s = NULL;
1411 : long i;
1412 :
1413 391567328 : for (i=a; i<b; i++)
1414 282031018 : if (gel(y,i) && gel(x,-i))
1415 : {
1416 203535753 : GEN t = gmul(gel(y,i), gel(x,-i));
1417 203532007 : s = s? gadd(s, t): t;
1418 : }
1419 109536310 : return RgX_coeff_simplify(av, s);
1420 : }
1421 :
1422 : /* assume nx >= ny > 0, return x * y * t^v */
1423 : static GEN
1424 33146769 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1425 : {
1426 : long i, lz, nz;
1427 : GEN z;
1428 :
1429 33146769 : x = RgXspec_kill0(x,nx);
1430 33146763 : y = RgXspec_kill0(y,ny);
1431 33146744 : lz = nx + ny + 1; nz = lz-2;
1432 33146744 : lz += v;
1433 33146744 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1434 71040814 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1435 81057495 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1436 54918975 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1437 47910831 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1438 33146438 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1439 : }
1440 :
1441 : /* return (x * X^d) + y. Assume d > 0 */
1442 : GEN
1443 9532927 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1444 : {
1445 : GEN x, y, xd, yd, zd;
1446 : long a, lz, nx, ny;
1447 :
1448 9532927 : if (!signe(x0)) return y0;
1449 9378919 : ny = lgpol(y0);
1450 9378914 : nx = lgpol(x0);
1451 9378982 : zd = (GEN)avma;
1452 9378982 : x = x0 + 2; y = y0 + 2; a = ny-d;
1453 9378982 : if (a <= 0)
1454 : {
1455 1478063 : lz = nx+d+2;
1456 1478063 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1457 3395305 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1458 1478063 : x = zd + a;
1459 1495802 : while (zd > x) gel(--zd,0) = gen_0;
1460 : }
1461 : else
1462 : {
1463 7900919 : xd = new_chunk(d); yd = y+d;
1464 7900916 : x = RgX_addspec_shallow(x,yd, nx,a);
1465 7900849 : lz = (a>nx)? ny+2: lg(x)+d;
1466 39487239 : x += 2; while (xd > x) *--zd = *--xd;
1467 : }
1468 22869330 : while (yd > y) *--zd = *--yd;
1469 9378912 : *--zd = x0[1];
1470 9378912 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1471 : }
1472 : GEN
1473 514932 : RgX_addmulXn(GEN x0, GEN y0, long d)
1474 : {
1475 : GEN x, y, xd, yd, zd;
1476 : long a, lz, nx, ny;
1477 :
1478 514932 : if (!signe(x0)) return RgX_copy(y0);
1479 514148 : nx = lgpol(x0);
1480 514148 : ny = lgpol(y0);
1481 514148 : zd = (GEN)avma;
1482 514148 : x = x0 + 2; y = y0 + 2; a = ny-d;
1483 514148 : if (a <= 0)
1484 : {
1485 291686 : lz = nx+d+2;
1486 291686 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1487 4293251 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1488 291686 : x = zd + a;
1489 757796 : while (zd > x) gel(--zd,0) = gen_0;
1490 : }
1491 : else
1492 : {
1493 222462 : xd = new_chunk(d); yd = y+d;
1494 222462 : x = RgX_addspec(x,yd, nx,a);
1495 222462 : lz = (a>nx)? ny+2: lg(x)+d;
1496 8418419 : x += 2; while (xd > x) *--zd = *--xd;
1497 : }
1498 2616086 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1499 514148 : *--zd = x0[1];
1500 514148 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1501 : }
1502 :
1503 : /* return x * y mod t^n */
1504 : static GEN
1505 6658643 : RgXn_mul_basecase(GEN x, GEN y, long n)
1506 : {
1507 6658643 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1508 : GEN z;
1509 6658643 : if (lx < 0) return pol_0(varn(x));
1510 6658643 : if (ly < 0) return pol_0(varn(x));
1511 6658643 : z = cgetg(lz, t_POL) + 2;
1512 6658643 : x+=2; if (lx > n) lx = n;
1513 6658643 : y+=2; if (ly > n) ly = n;
1514 6658643 : z[-1] = x[-1];
1515 6658643 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1516 6658643 : x = RgXspec_kill0(x, lx);
1517 6658643 : y = RgXspec_kill0(y, ly);
1518 : /* x:y:z [i] = term of degree i */
1519 26495993 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1520 11871517 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1521 6704792 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1522 6658643 : return normalizepol_lg(z - 2, lz);
1523 : }
1524 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1525 : static GEN
1526 9379769 : RgXn_mul2(GEN f, GEN g, long n)
1527 : {
1528 9379769 : pari_sp av = avma;
1529 : GEN fe,fo, ge,go, l,h,m;
1530 : long n0, n1;
1531 9379769 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1532 6694140 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1533 35497 : n0 = n>>1; n1 = n-n0;
1534 35497 : RgX_even_odd(f, &fe, &fo);
1535 35497 : RgX_even_odd(g, &ge, &go);
1536 35497 : l = RgXn_mul2(fe,ge,n1);
1537 35497 : h = RgXn_mul2(fo,go,n0);
1538 35497 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1539 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1540 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1541 35497 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1542 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1543 35497 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1544 35497 : m = RgX_inflate(m,2);
1545 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1546 35497 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1547 35497 : h = RgX_inflate(h,2);
1548 35497 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1549 35497 : return gc_upto(av, h);
1550 : }
1551 : /* (f*g) \/ x^n */
1552 : static GEN
1553 1595393 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1554 : {
1555 1595393 : long d = degpol(f)+degpol(g) + 1 - n;
1556 : GEN h;
1557 1595393 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1558 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1559 : RgX_recip_i(g), d));
1560 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1561 : }
1562 :
1563 : /* (f*g) \/ x^n */
1564 : static GEN
1565 0 : RgX_sqrhigh_i2(GEN f, long n)
1566 : {
1567 0 : long d = 2*degpol(f)+ 1 - n;
1568 : GEN h;
1569 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1570 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1571 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1572 : }
1573 :
1574 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1575 : * b+2 were sent instead. na, nb = number of terms of a, b.
1576 : * Only c, c0, c1, c2 are genuine GEN.
1577 : */
1578 : GEN
1579 33785891 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1580 : {
1581 : GEN a0, c, c0;
1582 33785891 : long n0, n0a, i, v = 0;
1583 : pari_sp av;
1584 :
1585 51860609 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1586 54044600 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1587 33785870 : if (na < nb) swapspec(a,b, na,nb);
1588 33785870 : if (!nb) return pol_0(0);
1589 :
1590 33626140 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1591 479388 : RgX_shift_inplace_init(v);
1592 479396 : i = (na>>1); n0 = na-i; na = i;
1593 479396 : av = avma; a0 = a+n0; n0a = n0;
1594 1334255 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1595 :
1596 479396 : if (nb > n0)
1597 : {
1598 : GEN b0,c1,c2;
1599 : long n0b;
1600 :
1601 478013 : nb -= n0; b0 = b+n0; n0b = n0;
1602 1445143 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1603 478013 : c = RgX_mulspec(a,b,n0a,n0b);
1604 478013 : c0 = RgX_mulspec(a0,b0, na,nb);
1605 :
1606 478013 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1607 478013 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1608 :
1609 478013 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1610 478013 : c2 = RgX_sub(c1, RgX_add(c0,c));
1611 478013 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1612 : }
1613 : else
1614 : {
1615 1383 : c = RgX_mulspec(a,b,n0a,nb);
1616 1383 : c0 = RgX_mulspec(a0,b,na,nb);
1617 : }
1618 479396 : c0 = RgX_addmulXn(c0,c,n0);
1619 479396 : return RgX_shift_inplace(gc_upto(av,c0), v);
1620 : }
1621 :
1622 : INLINE GEN
1623 94825 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1624 : {
1625 94825 : pari_sp av = avma;
1626 94825 : GEN s = NULL;
1627 94825 : long j, l = (i+1)>>1;
1628 197435 : for (j=a; j<l; j++)
1629 : {
1630 102610 : GEN xj = gel(x,j), xx = gel(x,i-j);
1631 102610 : if (xj && xx)
1632 : {
1633 93475 : GEN t = gmul(xj, xx);
1634 93475 : s = s? gadd(s, t): t;
1635 : }
1636 : }
1637 94825 : if (s) s = gshift(s,1);
1638 94825 : if ((i&1) == 0)
1639 : {
1640 64050 : GEN t = gel(x, i>>1);
1641 64050 : if (t) {
1642 60991 : t = gsqr(t);
1643 60991 : s = s? gadd(s, t): t;
1644 : }
1645 : }
1646 94825 : return RgX_coeff_simplify(av, s);
1647 : }
1648 : static GEN
1649 33275 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1650 : {
1651 : long i, lz, nz;
1652 : GEN z;
1653 :
1654 33275 : if (!nx) return pol_0(0);
1655 33275 : x = RgXspec_kill0(x,nx);
1656 33275 : lz = (nx << 1) + 1, nz = lz-2;
1657 33275 : lz += v;
1658 33275 : z = cgetg(lz,t_POL) + 2;
1659 88911 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1660 97325 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1661 64050 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1662 33275 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1663 : }
1664 : /* return x^2 mod t^n */
1665 : static GEN
1666 0 : RgXn_sqr_basecase(GEN x, long n)
1667 : {
1668 0 : long i, lz = n+2, lx = lgpol(x);
1669 : GEN z;
1670 0 : if (lx < 0) return pol_0(varn(x));
1671 0 : z = cgetg(lz, t_POL);
1672 0 : z[1] = x[1];
1673 0 : x+=2; if (lx > n) lx = n;
1674 0 : x = RgXspec_kill0(x,lx);
1675 0 : z+=2;/* x:z [i] = term of degree i */
1676 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1677 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1678 0 : z -= 2; return normalizepol_lg(z, lz);
1679 : }
1680 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1681 : static GEN
1682 315 : RgXn_sqr2(GEN f, long n)
1683 : {
1684 315 : pari_sp av = avma;
1685 : GEN fe,fo, l,h,m;
1686 : long n0, n1;
1687 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1688 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1689 0 : n0 = n>>1; n1 = n-n0;
1690 0 : RgX_even_odd(f, &fe, &fo);
1691 0 : l = RgXn_sqr(fe,n1);
1692 0 : h = RgXn_sqr(fo,n0);
1693 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1694 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1695 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1696 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1697 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1698 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1699 0 : m = RgX_inflate(m,2);
1700 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1701 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1702 0 : h = RgX_inflate(h,2);
1703 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1704 0 : return gc_upto(av, h);
1705 : }
1706 : GEN
1707 33314 : RgX_sqrspec(GEN a, long na)
1708 : {
1709 : GEN a0, c, c0, c1;
1710 33314 : long n0, n0a, i, v = 0;
1711 : pari_sp av;
1712 :
1713 61132 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1714 33314 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1715 39 : RgX_shift_inplace_init(v);
1716 39 : i = (na>>1); n0 = na-i; na = i;
1717 39 : av = avma; a0 = a+n0; n0a = n0;
1718 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1719 :
1720 39 : c = RgX_sqrspec(a,n0a);
1721 39 : c0 = RgX_sqrspec(a0,na);
1722 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1723 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1724 39 : c0 = RgX_addmulXn(c0,c,n0);
1725 39 : return RgX_shift_inplace(gc_upto(av,c0), v);
1726 : }
1727 :
1728 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1729 : GEN
1730 1720985 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1731 : {
1732 1720985 : GEN z = RgX_mul(A, B);
1733 1720977 : if (a < b)
1734 10334 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1735 1710643 : else if (a > b)
1736 1106964 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1737 : else
1738 603679 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1739 1720980 : return z;
1740 : }
1741 :
1742 : GEN
1743 32349109 : RgX_mul_i(GEN x, GEN y)
1744 : {
1745 32349109 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1746 32348799 : setvarn(z, varn(x)); return z;
1747 : }
1748 :
1749 : GEN
1750 33236 : RgX_sqr_i(GEN x)
1751 : {
1752 33236 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1753 33236 : setvarn(z,varn(x)); return z;
1754 : }
1755 :
1756 : /*******************************************************************/
1757 : /* */
1758 : /* DIVISION */
1759 : /* */
1760 : /*******************************************************************/
1761 : GEN
1762 5229745 : RgX_Rg_divexact(GEN x, GEN y) {
1763 5229745 : long i, lx = lg(x);
1764 : GEN z;
1765 5229745 : if (lx == 2) return gcopy(x);
1766 5199513 : switch(typ(y))
1767 : {
1768 5088933 : case t_INT:
1769 5088933 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1770 5023490 : break;
1771 5257 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1772 : }
1773 5128813 : z = cgetg(lx, t_POL); z[1] = x[1];
1774 29046073 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1775 5127538 : return z;
1776 : }
1777 : /* cf RgX_coeff_simplify */
1778 : INLINE GEN
1779 36712801 : RgX_coeff_div(GEN a, GEN b)
1780 : {
1781 36712801 : pari_sp av = avma;
1782 36712801 : GEN s = gdiv(a, b);
1783 36712750 : if (typ(s) != t_POL) return s;
1784 19885072 : switch(lg(s))
1785 : {
1786 7265418 : case 2: return gc_const(av, gen_0);
1787 9559363 : case 3: return gc_GEN(av, gel(s,2));
1788 3060291 : default: return s;
1789 : }
1790 : }
1791 : GEN
1792 29549746 : RgX_Rg_div(GEN x, GEN y) {
1793 29549746 : if (lg(x) == 2) return gcopy(x);
1794 29493917 : switch(typ(y))
1795 : {
1796 20326108 : case t_INT:
1797 20326108 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1798 4219274 : break;
1799 5740 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1800 : }
1801 50094091 : pari_APPLY_pol(RgX_coeff_div(gel(x,i), y));
1802 : }
1803 : GEN
1804 39711 : RgX_normalize(GEN x)
1805 : {
1806 39711 : GEN z, d = NULL;
1807 39711 : long i, n = lg(x)-1;
1808 39711 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1809 39711 : if (i == 1) return pol_0(varn(x));
1810 39711 : if (i == n && isint1(d)) return x;
1811 18200 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1812 32571 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1813 18200 : gel(z,n) = Rg_get_1(d); return z;
1814 : }
1815 : GEN
1816 10948 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1817 : GEN
1818 265228 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1819 : {
1820 265228 : long l = lg(a), i;
1821 : GEN a0, z0, z;
1822 :
1823 265228 : if (l <= 3)
1824 : {
1825 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1826 0 : return pol_0(varn(a));
1827 : }
1828 265228 : z = cgetg(l-1, t_POL);
1829 265227 : z[1] = a[1];
1830 265227 : a0 = a + l-1;
1831 265227 : z0 = z + l-2; *z0 = *a0--;
1832 3025205 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1833 : {
1834 2759979 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1835 2759978 : gel(z0,0) = t;
1836 : }
1837 265226 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1838 265226 : return z;
1839 : }
1840 : /* Polynomial division x / y:
1841 : * if pr = ONLY_REM return remainder, otherwise return quotient
1842 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1843 : * if pr != NULL set *pr to remainder, as the last object on stack */
1844 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1845 : static GEN
1846 23430270 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1847 : {
1848 : pari_sp avy, av, av1;
1849 : long dx,dy,dz,i,j,sx,lr;
1850 : GEN z,p1,p2,rem,y_lead,mod,p;
1851 : GEN (*f)(GEN,GEN);
1852 :
1853 23430270 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1854 :
1855 23430270 : dy = degpol(y);
1856 23430248 : y_lead = gel(y,dy+2);
1857 23430248 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1858 : {
1859 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1860 0 : for (dy--; dy>=0; dy--)
1861 : {
1862 0 : y_lead = gel(y,dy+2);
1863 0 : if (!gequal0(y_lead)) break;
1864 : }
1865 : }
1866 23430247 : if (!dy) /* y is constant */
1867 : {
1868 6816 : if (pr == ONLY_REM) return pol_0(varn(x));
1869 6809 : z = RgX_Rg_div(x, y_lead);
1870 6809 : if (pr == ONLY_DIVIDES) return z;
1871 2595 : if (pr) *pr = pol_0(varn(x));
1872 2595 : return z;
1873 : }
1874 23423431 : dx = degpol(x);
1875 23423395 : if (dx < dy)
1876 : {
1877 3845776 : if (pr == ONLY_REM) return RgX_copy(x);
1878 346015 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1879 345994 : z = pol_0(varn(x));
1880 345994 : if (pr) *pr = RgX_copy(x);
1881 345994 : return z;
1882 : }
1883 :
1884 : /* x,y in R[X], y non constant */
1885 19577619 : av = avma;
1886 19577619 : p = NULL;
1887 19577619 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1888 : {
1889 127267 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1890 127267 : if (!z) return gc_NULL(av);
1891 127267 : z = FpX_to_mod(z, p);
1892 127267 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1893 68579 : return gc_upto(av, z);
1894 58688 : *pr = FpX_to_mod(*pr, p);
1895 58688 : return gc_all(av, 2, &z, pr);
1896 : }
1897 19450463 : switch(typ(y_lead))
1898 : {
1899 0 : case t_REAL:
1900 0 : y_lead = ginv(y_lead);
1901 0 : f = gmul; mod = NULL;
1902 0 : break;
1903 2788 : case t_INTMOD:
1904 2788 : case t_POLMOD: y_lead = ginv(y_lead);
1905 2788 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1906 2788 : break;
1907 19447675 : default: if (gequal1(y_lead)) y_lead = NULL;
1908 19447669 : f = gdiv; mod = NULL;
1909 : }
1910 :
1911 19450457 : if (y_lead == NULL)
1912 17489803 : p2 = gel(x,dx+2);
1913 : else {
1914 : for(;;) {
1915 1960654 : p2 = f(gel(x,dx+2),y_lead);
1916 1960654 : p2 = simplify_shallow(p2);
1917 1960654 : if (!isexactzero(p2) || (--dx < 0)) break;
1918 : }
1919 1960654 : if (dx < dy) /* leading coeff of x was in fact zero */
1920 : {
1921 0 : if (pr == ONLY_DIVIDES) {
1922 0 : set_avma(av);
1923 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1924 : }
1925 0 : if (pr == ONLY_REM)
1926 : {
1927 0 : if (dx < 0)
1928 0 : return gc_GEN(av, scalarpol(p2, varn(x)));
1929 : else
1930 : {
1931 : GEN t;
1932 0 : set_avma(av);
1933 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1934 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1935 0 : return t;
1936 : }
1937 : }
1938 0 : if (pr) /* cf ONLY_REM above */
1939 : {
1940 0 : if (dx < 0)
1941 : {
1942 0 : p2 = gclone(p2);
1943 0 : set_avma(av);
1944 0 : z = pol_0(varn(x));
1945 0 : x = scalarpol(p2, varn(x));
1946 0 : gunclone(p2);
1947 : }
1948 : else
1949 : {
1950 : GEN t;
1951 0 : set_avma(av);
1952 0 : z = pol_0(varn(x));
1953 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1954 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1955 0 : x = t;
1956 : }
1957 0 : *pr = x;
1958 : }
1959 : else
1960 : {
1961 0 : set_avma(av);
1962 0 : z = pol_0(varn(x));
1963 : }
1964 0 : return z;
1965 : }
1966 : }
1967 : /* dx >= dy */
1968 19450457 : avy = avma;
1969 19450457 : dz = dx-dy;
1970 19450457 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1971 19450452 : x += 2;
1972 19450452 : z += 2;
1973 19450452 : y += 2;
1974 19450452 : gel(z,dz) = gcopy(p2);
1975 :
1976 53708065 : for (i=dx-1; i>=dy; i--)
1977 : {
1978 34257815 : av1=avma; p1=gel(x,i);
1979 1139425652 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1980 34254087 : if (y_lead) p1 = simplify(f(p1,y_lead));
1981 :
1982 34254087 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1983 : else
1984 24266293 : p1 = avma==av1? gcopy(p1): gc_upto(av1,p1);
1985 34257220 : gel(z,i-dy) = p1;
1986 : }
1987 19450250 : if (!pr) return gc_upto(av,z-2);
1988 :
1989 12526180 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1990 12881246 : for (sx=0; ; i--)
1991 : {
1992 12881246 : p1 = gel(x,i);
1993 : /* we always enter this loop at least once */
1994 32143809 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1995 12880507 : if (mod && avma==av1) p1 = gmul(p1,mod);
1996 12880507 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1997 3574812 : if (!isexactzero(p1)) break;
1998 3424615 : if (!i) break;
1999 355035 : set_avma(av1);
2000 : }
2001 12525748 : if (pr == ONLY_DIVIDES)
2002 : {
2003 8078 : if (sx) return gc_NULL(av);
2004 8057 : set_avma((pari_sp)rem); return gc_upto(av,z-2);
2005 : }
2006 12517670 : lr=i+3; rem -= lr;
2007 12517670 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
2008 12408374 : else p1 = gc_upto((pari_sp)rem,p1);
2009 12518100 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2010 12518100 : rem[1] = z[-1];
2011 12518100 : rem += 2;
2012 12518100 : gel(rem,i) = p1;
2013 17228702 : for (i--; i>=0; i--)
2014 : {
2015 4710612 : av1=avma; p1 = gel(x,i);
2016 13575683 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2017 4710249 : if (mod && avma==av1) p1 = gmul(p1,mod);
2018 4710425 : gel(rem,i) = avma==av1? gcopy(p1):gc_upto(av1,p1);
2019 : }
2020 12518090 : rem -= 2;
2021 12518090 : if (!sx) (void)normalizepol_lg(rem, lr);
2022 12518264 : if (pr == ONLY_REM) return gc_upto(av,rem);
2023 6763584 : z -= 2; *pr = rem; return gc_all_unsafe(av,avy,2,&z,pr);
2024 : }
2025 :
2026 : GEN
2027 14175868 : RgX_divrem(GEN x, GEN y, GEN *pr)
2028 : {
2029 14175868 : if (pr == ONLY_REM) return RgX_rem(x, y);
2030 14175868 : return RgX_divrem_i(x, y, pr);
2031 : }
2032 :
2033 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
2034 : GEN
2035 156752 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
2036 : {
2037 156752 : long vx = varn(x), dx = degpol(x), dy = degpol(y), dz, i, j, sx, lr;
2038 : pari_sp av0, av;
2039 : GEN z, p1, rem, lead;
2040 :
2041 156752 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
2042 156752 : if (dx < dy)
2043 : {
2044 41644 : if (pr)
2045 : {
2046 41616 : av0 = avma; x = RgXQX_red(x, T);
2047 41616 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
2048 41609 : if (pr == ONLY_REM) return x;
2049 0 : *pr = x;
2050 : }
2051 28 : return pol_0(vx);
2052 : }
2053 115108 : lead = leading_coeff(y);
2054 115108 : if (!dy) /* y is constant */
2055 : {
2056 602 : if (pr && pr != ONLY_DIVIDES)
2057 : {
2058 0 : if (pr == ONLY_REM) return pol_0(vx);
2059 0 : *pr = pol_0(vx);
2060 : }
2061 602 : if (gequal1(lead)) return RgX_copy(x);
2062 0 : av0 = avma; x = gmul(x, ginvmod(lead,T));
2063 0 : return gc_upto(av0, RgXQX_red(x,T));
2064 : }
2065 114506 : av0 = avma; dz = dx-dy;
2066 114506 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
2067 114506 : set_avma(av0);
2068 114506 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2069 114506 : x += 2; y += 2; z += 2;
2070 :
2071 114506 : p1 = gel(x,dx); av = avma;
2072 114506 : gel(z,dz) = lead? gc_upto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2073 581322 : for (i=dx-1; i>=dy; i--)
2074 : {
2075 466816 : av = avma; p1 = gel(x,i);
2076 2416854 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2077 466796 : if (lead) p1 = gmul(grem(p1, T), lead);
2078 466795 : gel(z,i-dy) = gc_upto(av, grem(p1, T));
2079 : }
2080 114506 : if (!pr) { guncloneNULL(lead); return z-2; }
2081 :
2082 113232 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2083 191267 : for (sx=0; ; i--)
2084 : {
2085 191267 : p1 = gel(x,i);
2086 661291 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2087 191269 : p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2088 105538 : if (!i) break;
2089 78035 : set_avma(av);
2090 : }
2091 113232 : if (pr == ONLY_DIVIDES)
2092 : {
2093 27038 : guncloneNULL(lead);
2094 27038 : if (sx) return gc_NULL(av0);
2095 24479 : return gc_const((pari_sp)rem, z-2);
2096 : }
2097 86194 : lr=i+3; rem -= lr; av = (pari_sp)rem;
2098 86194 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2099 86194 : rem[1] = z[-1];
2100 86194 : rem += 2; gel(rem,i) = gc_upto(av, p1);
2101 187610 : for (i--; i>=0; i--)
2102 : {
2103 101416 : av = avma; p1 = gel(x,i);
2104 337537 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2105 101416 : gel(rem,i) = gc_upto(av, grem(p1, T));
2106 : }
2107 86194 : rem -= 2;
2108 86194 : guncloneNULL(lead);
2109 86194 : if (!sx) (void)normalizepol_lg(rem, lr);
2110 86194 : if (pr == ONLY_REM) return gc_upto(av0,rem);
2111 168 : *pr = rem; return z-2;
2112 : }
2113 :
2114 : /*******************************************************************/
2115 : /* */
2116 : /* PSEUDO-DIVISION */
2117 : /* */
2118 : /*******************************************************************/
2119 : INLINE GEN
2120 1050120 : rem(GEN c, GEN T)
2121 : {
2122 1050120 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2123 1050120 : return c;
2124 : }
2125 :
2126 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2127 : int
2128 17720 : ZXQX_dvd(GEN x, GEN y, GEN T)
2129 : {
2130 : long dx, dy, i, T_ismonic;
2131 17720 : pari_sp av = avma, av2;
2132 : GEN y_lead;
2133 :
2134 17720 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2135 17720 : dy = degpol(y); y_lead = gel(y,dy+2);
2136 17720 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2137 : /* if monic, no point in using pseudo-division */
2138 17720 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2139 14738 : T_ismonic = gequal1(leading_coeff(T));
2140 14738 : dx = degpol(x);
2141 14738 : if (dx < dy) return !signe(x);
2142 14738 : (void)new_chunk(2);
2143 14738 : x = RgX_recip_i(x)+2;
2144 14738 : y = RgX_recip_i(y)+2;
2145 : /* pay attention to sparse divisors */
2146 29977 : for (i = 1; i <= dy; i++)
2147 15239 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2148 14738 : av2 = avma;
2149 : for (;;)
2150 72953 : {
2151 87691 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2152 87691 : x0 = gneg(x0);
2153 87691 : m = gcdii(cx, y0);
2154 87691 : if (!equali1(m))
2155 : {
2156 85189 : x0 = gdiv(x0, m);
2157 85189 : y0 = diviiexact(y0, m);
2158 85189 : if (equali1(y0)) y0 = NULL;
2159 : }
2160 179155 : for (i=1; i<=dy; i++)
2161 : {
2162 91464 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2163 91464 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2164 91464 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2165 91464 : gel(x,i) = c;
2166 : }
2167 687972 : for ( ; i<=dx; i++)
2168 : {
2169 600281 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2170 600281 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2171 600281 : gel(x,i) = c;
2172 : }
2173 102167 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2174 87691 : if (dx < dy) break;
2175 72953 : if (gc_needed(av2,1))
2176 : {
2177 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2178 28 : gc_slice(av2,x,dx+1);
2179 : }
2180 : }
2181 14738 : return gc_bool(av, dx < 0);
2182 : }
2183 :
2184 : /* T either NULL or a t_POL. */
2185 : GEN
2186 138802 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2187 : {
2188 138802 : long vx = varn(x), dx, dy, dz, i, lx, p;
2189 138802 : pari_sp av = avma, av2;
2190 : GEN y_lead;
2191 :
2192 138802 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2193 138802 : dy = degpol(y); y_lead = gel(y,dy+2);
2194 : /* if monic, no point in using pseudo-division */
2195 138802 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2196 125088 : dx = degpol(x);
2197 125088 : if (dx < dy) return RgX_copy(x);
2198 125081 : (void)new_chunk(2);
2199 125081 : x = RgX_recip_i(x)+2;
2200 125081 : y = RgX_recip_i(y)+2;
2201 : /* pay attention to sparse divisors */
2202 535552 : for (i = 1; i <= dy; i++)
2203 410471 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2204 125081 : dz = dx-dy; p = dz+1;
2205 125081 : av2 = avma;
2206 : for (;;)
2207 : {
2208 247146 : gel(x,0) = gneg(gel(x,0)); p--;
2209 1061829 : for (i=1; i<=dy; i++)
2210 : {
2211 814686 : GEN c = gmul(y_lead, gel(x,i));
2212 814687 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2213 814682 : gel(x,i) = rem(c, T);
2214 : }
2215 373890 : for ( ; i<=dx; i++)
2216 : {
2217 126742 : GEN c = gmul(y_lead, gel(x,i));
2218 126741 : gel(x,i) = rem(c, T);
2219 : }
2220 253347 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2221 247146 : if (dx < dy) break;
2222 122065 : if (gc_needed(av2,1))
2223 : {
2224 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2225 0 : gc_slice(av2,x,dx+1);
2226 : }
2227 : }
2228 125081 : if (dx < 0) return pol_0(vx);
2229 124920 : lx = dx+3; x -= 2;
2230 124920 : x[0] = evaltyp(t_POL) | _evallg(lx);
2231 124920 : x[1] = evalsigne(1) | evalvarn(vx);
2232 124920 : x = RgX_recip_i(x);
2233 124920 : if (p)
2234 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2235 3884 : GEN t = y_lead;
2236 3884 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2237 0 : t = RgXQ_powu(t, p, T);
2238 : else
2239 3884 : t = gpowgs(t, p);
2240 11803 : for (i=2; i<lx; i++)
2241 : {
2242 7919 : GEN c = gmul(gel(x,i), t);
2243 7919 : gel(x,i) = rem(c,T);
2244 : }
2245 3884 : if (!T) return gc_upto(av, x);
2246 : }
2247 121036 : return gc_GEN(av, x);
2248 : }
2249 :
2250 : GEN
2251 138802 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2252 :
2253 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2254 : GEN
2255 11923 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2256 : {
2257 11923 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2258 11923 : pari_sp av = avma, av2;
2259 : GEN z, r, ypow, y_lead;
2260 :
2261 11923 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2262 11923 : dy = degpol(y); y_lead = gel(y,dy+2);
2263 11923 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2264 7421 : dx = degpol(x);
2265 7421 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2266 7421 : if (dx == dy)
2267 : {
2268 98 : GEN x_lead = gel(x,lg(x)-1);
2269 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2270 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2271 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2272 98 : *ptr = gc_upto(av, r); return scalarpol(x_lead, vx);
2273 : }
2274 7323 : (void)new_chunk(2);
2275 7323 : x = RgX_recip_i(x)+2;
2276 7323 : y = RgX_recip_i(y)+2;
2277 : /* pay attention to sparse divisors */
2278 38650 : for (i = 1; i <= dy; i++)
2279 31327 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2280 7323 : dz = dx-dy; p = dz+1;
2281 7323 : lz = dz+3;
2282 7323 : z = cgetg(lz, t_POL);
2283 7323 : z[1] = evalsigne(1) | evalvarn(vx);
2284 28220 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2285 7323 : ypow = new_chunk(dz+1);
2286 7323 : gel(ypow,0) = gen_1;
2287 7323 : gel(ypow,1) = y_lead;
2288 13574 : for (i=2; i<=dz; i++)
2289 : {
2290 6251 : GEN c = gmul(gel(ypow,i-1), y_lead);
2291 6251 : gel(ypow,i) = rem(c,T);
2292 : }
2293 7323 : av2 = avma;
2294 7323 : for (iz=2;;)
2295 : {
2296 15325 : p--;
2297 15325 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2298 69215 : for (i=1; i<=dy; i++)
2299 : {
2300 53890 : GEN c = gmul(y_lead, gel(x,i));
2301 53890 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2302 53890 : gel(x,i) = rem(c, T);
2303 : }
2304 40638 : for ( ; i<=dx; i++)
2305 : {
2306 25313 : GEN c = gmul(y_lead, gel(x,i));
2307 25313 : gel(x,i) = rem(c,T);
2308 : }
2309 15325 : x++; dx--;
2310 20897 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2311 15325 : if (dx < dy) break;
2312 8002 : if (gc_needed(av2,1))
2313 : {
2314 0 : GEN X = x-2;
2315 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2316 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2317 0 : (void)gc_all(av2,2, &X, &z); x = X+2;
2318 : }
2319 : }
2320 13868 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2321 7323 : if (dx < 0)
2322 182 : x = pol_0(vx);
2323 : else
2324 : {
2325 7141 : lx = dx+3; x -= 2;
2326 7141 : x[0] = evaltyp(t_POL) | _evallg(lx);
2327 7141 : x[1] = evalsigne(1) | evalvarn(vx);
2328 7141 : x = RgX_recip_i(x);
2329 : }
2330 7323 : z = RgX_recip_i(z);
2331 7323 : r = x;
2332 7323 : if (p)
2333 : {
2334 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2335 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2336 : }
2337 7323 : *ptr = r; return gc_all(av, 2, &z, ptr);
2338 : }
2339 : GEN
2340 11678 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2341 11678 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2342 :
2343 : GEN
2344 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2345 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2346 :
2347 : /* cf RgX_coeff_simplify */
2348 : INLINE GEN
2349 727135933 : RgX_coeff_mul(GEN a, GEN b)
2350 : {
2351 727135933 : pari_sp av = avma;
2352 727135933 : GEN s = gmul(a, b);
2353 727099957 : if (typ(s) != t_POL) return s;
2354 129325864 : switch(lg(s))
2355 : {
2356 72205141 : case 2: return gc_const(av, gen_0);
2357 729939 : case 3: return gc_GEN(av, gel(s,2));
2358 56390784 : default: return s;
2359 : }
2360 : }
2361 : GEN
2362 963277721 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(RgX_coeff_mul(y, gel(x,i))); }
2363 : GEN
2364 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2365 : GEN
2366 26285 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2367 : GEN
2368 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2369 : GEN
2370 168 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2371 :
2372 : GEN
2373 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2374 :
2375 : GEN
2376 0 : RgXQX_powers(GEN P, long n, GEN T)
2377 : {
2378 0 : GEN v = cgetg(n+2, t_VEC);
2379 : long i;
2380 0 : gel(v, 1) = pol_1(varn(T));
2381 0 : if (n==0) return v;
2382 0 : gel(v, 2) = gcopy(P);
2383 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2384 0 : return v;
2385 : }
2386 :
2387 : static GEN
2388 623238 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2389 : static GEN
2390 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2391 : static GEN
2392 74806 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2393 : static GEN
2394 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2395 : static GEN
2396 280956 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2397 : static GEN
2398 1040375 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2399 : static GEN
2400 1019355 : _one(void *data) { return pol_1(varn((GEN)data)); }
2401 : static GEN
2402 1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2403 : static GEN
2404 724099 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2405 :
2406 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2407 : _mul, _sqr, _one, _zero };
2408 :
2409 : GEN
2410 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2411 : {
2412 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2413 : }
2414 :
2415 : GEN
2416 417402 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2417 : {
2418 417402 : int use_sqr = 2*degpol(x) >= degpol(T);
2419 417401 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2420 : }
2421 :
2422 : /* mod X^n */
2423 : struct modXn {
2424 : long v; /* varn(X) */
2425 : long n;
2426 : } ;
2427 : static GEN
2428 11893 : _sqrXn(void *data, GEN x) {
2429 11893 : struct modXn *S = (struct modXn*)data;
2430 11893 : return RgXn_sqr(x, S->n);
2431 : }
2432 : static GEN
2433 4528 : _mulXn(void *data, GEN x, GEN y) {
2434 4528 : struct modXn *S = (struct modXn*)data;
2435 4528 : return RgXn_mul(x,y, S->n);
2436 : }
2437 : static GEN
2438 1939 : _oneXn(void *data) {
2439 1939 : struct modXn *S = (struct modXn*)data;
2440 1939 : return pol_1(S->v);
2441 : }
2442 : static GEN
2443 0 : _zeroXn(void *data) {
2444 0 : struct modXn *S = (struct modXn*)data;
2445 0 : return pol_0(S->v);
2446 : }
2447 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2448 : _oneXn, _zeroXn };
2449 :
2450 : GEN
2451 483 : RgXn_powers(GEN x, long m, long n)
2452 : {
2453 483 : long d = degpol(x);
2454 483 : int use_sqr = (d<<1) >= n;
2455 : struct modXn S;
2456 483 : S.v = varn(x); S.n = n;
2457 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2458 : }
2459 :
2460 : GEN
2461 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2462 : {
2463 : struct modXn S;
2464 : long v;
2465 2286 : if (n == 0) return x;
2466 2286 : v = RgX_valrem(x, &x);
2467 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2468 2265 : S.v = varn(x); S.n = n;
2469 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2470 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2471 2265 : return x;
2472 : }
2473 : GEN
2474 0 : RgXn_powu(GEN x, ulong m, long n)
2475 : {
2476 : pari_sp av;
2477 0 : if (n == 0) return gcopy(x);
2478 0 : av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
2479 : }
2480 :
2481 : GEN
2482 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2483 : {
2484 : struct modXn S;
2485 966 : S.v = varn(gel(x,2)); S.n = n;
2486 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2487 : }
2488 :
2489 : GEN
2490 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2491 : {
2492 0 : int use_sqr = 2*degpol(x) >= n;
2493 : struct modXn S;
2494 0 : S.v = varn(x); S.n = n;
2495 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2496 : }
2497 :
2498 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2499 : GEN
2500 5593 : RgXn_eval(GEN Q, GEN x, long n)
2501 : {
2502 5593 : long d = degpol(x);
2503 : int use_sqr;
2504 : struct modXn S;
2505 5593 : if (d == 1 && isrationalzero(gel(x,2)))
2506 : {
2507 5586 : GEN y = RgX_unscale(Q, gel(x,3));
2508 5586 : setvarn(y, varn(x)); return y;
2509 : }
2510 7 : S.v = varn(x);
2511 7 : S.n = n;
2512 7 : use_sqr = (d<<1) >= n;
2513 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2514 : }
2515 :
2516 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2517 : static GEN
2518 2569747 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2519 : {
2520 2569747 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2521 2569747 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2522 : }
2523 :
2524 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2525 : static GEN
2526 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2527 : {
2528 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2529 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2530 : }
2531 :
2532 : static GEN
2533 2141046 : RgXn_div_gen(GEN g, GEN f, long e)
2534 : {
2535 : pari_sp av;
2536 : ulong mask;
2537 : GEN W, a;
2538 2141046 : long v = varn(f), n = 1;
2539 :
2540 2141046 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2541 2141039 : a = ginv(gel(f,2));
2542 2141034 : if (e == 1 && !g) return scalarpol(a, v);
2543 2135385 : else if (e == 2 && !g)
2544 : {
2545 : GEN b;
2546 836881 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2547 246362 : b = gneg(b);
2548 246363 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2549 246363 : return deg1pol(b, a, v);
2550 : }
2551 1298504 : av = avma;
2552 1298504 : W = scalarpol_shallow(a,v);
2553 1298504 : mask = quadratic_prec_mask(e);
2554 3862294 : while (mask > 1)
2555 : {
2556 : GEN u, fr;
2557 2563790 : long n2 = n;
2558 2563790 : n<<=1; if (mask & 1) n--;
2559 2563790 : mask >>= 1;
2560 2563790 : fr = RgXn_red_shallow(f, n);
2561 2563790 : if (mask>1 || !g)
2562 : {
2563 1329067 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2564 1329067 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2565 : }
2566 : else
2567 : {
2568 1234723 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2569 1234723 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2570 1234723 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2571 : }
2572 2563790 : if (gc_needed(av,2))
2573 : {
2574 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2575 0 : W = gc_upto(av, W);
2576 : }
2577 : }
2578 1298504 : return W;
2579 : }
2580 :
2581 : static GEN
2582 98 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2583 : {
2584 : GEN r;
2585 98 : if (lgefint(p) == 3)
2586 : {
2587 98 : ulong pp = uel(p, 2);
2588 98 : if (pp == 2)
2589 7 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2590 : else
2591 91 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2592 : }
2593 : else
2594 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2595 98 : return FpX_to_mod(r, p);
2596 : }
2597 :
2598 : static GEN
2599 0 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2600 : {
2601 0 : GEN r, T = RgX_to_FpX(pol, p);
2602 0 : if (signe(T) == 0) pari_err_OP("/", x, y);
2603 0 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2604 0 : return FpXQX_to_mod(r, T, p);
2605 : }
2606 :
2607 : static GEN
2608 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2609 : {
2610 : GEN r;
2611 91 : if (lgefint(p) == 3)
2612 : {
2613 91 : ulong pp = uel(p, 2);
2614 91 : if (pp == 2)
2615 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2616 : else
2617 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2618 : }
2619 : else
2620 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2621 91 : return FpX_to_mod(r, p);
2622 : }
2623 :
2624 : static GEN
2625 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2626 : {
2627 0 : GEN r, T = RgX_to_FpX(pol, p);
2628 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2629 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2630 0 : return FpXQX_to_mod(r, T, p);
2631 : }
2632 :
2633 : static GEN
2634 906415 : RgXn_inv_fast(GEN x, long e)
2635 : {
2636 : GEN p, pol;
2637 : long pa;
2638 906415 : long t = RgX_type(x,&p,&pol,&pa);
2639 906415 : switch(t)
2640 : {
2641 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2642 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2643 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2644 906324 : default: return NULL;
2645 : }
2646 : }
2647 :
2648 : static GEN
2649 1234821 : RgXn_div_fast(GEN x, GEN y, long e)
2650 : {
2651 : GEN p, pol;
2652 : long pa;
2653 1234821 : long t = RgX_type2(x,y,&p,&pol,&pa);
2654 1234821 : switch(t)
2655 : {
2656 98 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2657 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2658 0 : return RgXn_div_FpXQX(x, y, e, pol, p);
2659 1234723 : default: return NULL;
2660 : }
2661 : }
2662 :
2663 : GEN
2664 1234821 : RgXn_div_i(GEN g, GEN f, long e)
2665 : {
2666 1234821 : GEN h = RgXn_div_fast(g, f, e);
2667 1234821 : if (h) return h;
2668 1234723 : return RgXn_div_gen(g, f, e);
2669 : }
2670 :
2671 : GEN
2672 532507 : RgXn_div(GEN g, GEN f, long e)
2673 : {
2674 532507 : pari_sp av = avma;
2675 532507 : return gc_upto(av, RgXn_div_i(g, f, e));
2676 : }
2677 :
2678 : GEN
2679 906415 : RgXn_inv_i(GEN f, long e)
2680 : {
2681 906415 : GEN h = RgXn_inv_fast(f, e);
2682 906415 : if (h) return h;
2683 906324 : return RgXn_div_gen(NULL, f, e);
2684 : }
2685 :
2686 : GEN
2687 807482 : RgXn_inv(GEN f, long e)
2688 : {
2689 807482 : pari_sp av = avma;
2690 807482 : return gc_upto(av, RgXn_inv_i(f, e));
2691 : }
2692 :
2693 : /* intformal(x^n*S) / x^(n+1) */
2694 : static GEN
2695 56375 : RgX_integXn(GEN x, long n)
2696 114353 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2697 :
2698 : GEN
2699 52919 : RgXn_expint(GEN h, long e)
2700 : {
2701 52919 : pari_sp av = avma, av2;
2702 52919 : long v = varn(h), n;
2703 52919 : GEN f = pol_1(v), g;
2704 : ulong mask;
2705 :
2706 52919 : if (!signe(h)) return f;
2707 50427 : g = pol_1(v);
2708 50427 : n = 1; mask = quadratic_prec_mask(e);
2709 50427 : av2 = avma;
2710 56377 : for (;mask>1;)
2711 : {
2712 : GEN u, w;
2713 56377 : long n2 = n;
2714 56377 : n<<=1; if (mask & 1) n--;
2715 56377 : mask >>= 1;
2716 56377 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2717 56374 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2718 56375 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2719 56376 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2720 56377 : if (mask<=1) break;
2721 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2722 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2723 5950 : if (gc_needed(av2,2))
2724 : {
2725 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2726 0 : (void)gc_all(av2, 2, &f, &g);
2727 : }
2728 : }
2729 50427 : return gc_upto(av, f);
2730 : }
2731 :
2732 : GEN
2733 0 : RgXn_exp(GEN h, long e)
2734 : {
2735 0 : long d = degpol(h);
2736 0 : if (d < 0) return pol_1(varn(h));
2737 0 : if (!d || !gequal0(gel(h,2)))
2738 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2739 0 : return RgXn_expint(RgX_deriv(h), e);
2740 : }
2741 :
2742 : GEN
2743 154 : RgXn_reverse(GEN f, long e)
2744 : {
2745 154 : pari_sp av = avma, av2;
2746 : ulong mask;
2747 : GEN fi, a, df, W, an;
2748 154 : long v = varn(f), n=1;
2749 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2750 0 : pari_err_INV("serreverse",f);
2751 154 : fi = ginv(gel(f,3));
2752 154 : a = deg1pol_shallow(fi,gen_0,v);
2753 154 : if (e <= 2) return gc_GEN(av, a);
2754 133 : W = scalarpol(fi,v);
2755 133 : df = RgX_deriv(f);
2756 133 : mask = quadratic_prec_mask(e);
2757 133 : av2 = avma;
2758 616 : for (;mask>1;)
2759 : {
2760 : GEN u, fa, fr;
2761 483 : long n2 = n, rt;
2762 483 : n<<=1; if (mask & 1) n--;
2763 483 : mask >>= 1;
2764 483 : fr = RgXn_red_shallow(f, n);
2765 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2766 483 : an = RgXn_powers(a, rt, n);
2767 483 : if (n>1)
2768 : {
2769 483 : long n4 = (n2+1)>>1;
2770 483 : GEN dfr = RgXn_red_shallow(df, n2);
2771 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2772 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2773 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2774 : }
2775 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2776 483 : fa = RgX_shift(fa, -n2);
2777 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2778 483 : if (gc_needed(av2,2))
2779 : {
2780 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2781 0 : (void)gc_all(av2, 2, &a, &W);
2782 : }
2783 : }
2784 133 : return gc_upto(av, a);
2785 : }
2786 :
2787 : GEN
2788 7 : RgXn_sqrt(GEN h, long e)
2789 : {
2790 7 : pari_sp av = avma, av2;
2791 7 : long v = varn(h), n = 1;
2792 7 : GEN f = scalarpol(gen_1, v), df = f;
2793 7 : ulong mask = quadratic_prec_mask(e);
2794 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2795 0 : pari_err_SQRTN("RgXn_sqrt",h);
2796 7 : av2 = avma;
2797 : while(1)
2798 7 : {
2799 14 : long n2 = n, m;
2800 : GEN g;
2801 14 : n<<=1; if (mask & 1) n--;
2802 14 : mask >>= 1;
2803 14 : m = n-n2;
2804 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2805 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2806 14 : if (mask==1) return gc_upto(av, f);
2807 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2808 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2809 7 : if (gc_needed(av2,2))
2810 : {
2811 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2812 0 : (void)gc_all(av2, 2, &f, &df);
2813 : }
2814 : }
2815 : }
2816 :
2817 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2818 : GEN
2819 120400 : RgXQ_powu(GEN x, ulong n, GEN T)
2820 : {
2821 120400 : pari_sp av = avma;
2822 :
2823 120400 : if (!n) return pol_1(varn(x));
2824 72156 : if (n == 1) return RgX_copy(x);
2825 23527 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2826 23527 : return gc_GEN(av, x);
2827 : }
2828 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2829 : GEN
2830 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2831 : {
2832 : pari_sp av;
2833 102160 : long s = signe(n);
2834 :
2835 102160 : if (!s) return pol_1(varn(x));
2836 102160 : if (is_pm1(n) == 1)
2837 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2838 13853 : av = avma;
2839 13853 : if (s < 0) x = RgXQ_inv(x, T);
2840 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2841 13853 : return gc_GEN(av, x);
2842 : }
2843 : static GEN
2844 200590 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2845 : static GEN
2846 111181 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2847 :
2848 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2849 : GEN
2850 12061 : ZXQ_powers(GEN x, long l, GEN T)
2851 : {
2852 12061 : int use_sqr = 2*degpol(x) >= degpol(T);
2853 12061 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2854 : }
2855 :
2856 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2857 : GEN
2858 168511 : ZXQ_powu(GEN x, ulong n, GEN T)
2859 : {
2860 168511 : pari_sp av = avma;
2861 :
2862 168511 : if (!n) return pol_1(varn(x));
2863 168511 : if (n == 1) return ZX_copy(x);
2864 113875 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2865 113873 : return gc_GEN(av, x);
2866 : }
2867 :
2868 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2869 : GEN
2870 8603 : RgXQ_powers(GEN x, long l, GEN T)
2871 : {
2872 8603 : int use_sqr = 2*degpol(x) >= degpol(T);
2873 8603 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2874 : }
2875 :
2876 : GEN
2877 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2878 : {
2879 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2880 : }
2881 :
2882 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2883 : GEN
2884 9506 : QXQ_powers(GEN a, long n, GEN T)
2885 : {
2886 : GEN den, v;
2887 9506 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2888 9478 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2889 : /* den*a integral; v[i+1] = (den*a)^i in K */
2890 9478 : if (den)
2891 : { /* restore denominators */
2892 5942 : GEN d = den;
2893 : long i;
2894 5942 : gel(v,2) = a;
2895 27679 : for (i=3; i<=n+1; i++) {
2896 21737 : d = mulii(d,den);
2897 21737 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2898 : }
2899 : }
2900 9478 : return v;
2901 : }
2902 :
2903 : static GEN
2904 3542 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2905 : {
2906 3542 : long l, i, m = 0;
2907 : GEN dz, z;
2908 3542 : GEN V = cgetg_copy(v, &l);
2909 12054 : for (i = imin; i < l; i++)
2910 : {
2911 8512 : GEN c = gel(v, i);
2912 8512 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2913 : }
2914 3542 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2915 3836 : for (i = 1; i < imin; i++) V[i] = v[i];
2916 12054 : for (i = imin; i < l; i++)
2917 : {
2918 8512 : GEN c = gel(v,i);
2919 8512 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2920 8512 : gel(V,i) = c;
2921 : }
2922 3542 : return V;
2923 : }
2924 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2925 : GEN
2926 3248 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2927 3248 : { return do_QXQ_eval(v, 1, a, T); }
2928 :
2929 : GEN
2930 294 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2931 294 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2932 :
2933 : GEN
2934 2933 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2935 : {
2936 2933 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2937 : }
2938 :
2939 : GEN
2940 2659 : RgXQ_norm(GEN x, GEN T)
2941 : {
2942 : pari_sp av;
2943 2659 : long dx = degpol(x);
2944 : GEN L, y;
2945 2659 : if (degpol(T)==0) return gpowgs(x,0);
2946 2652 : av = avma; y = resultant(T, x);
2947 2652 : L = leading_coeff(T);
2948 2652 : if (gequal1(L) || !signe(x)) return y;
2949 0 : return gc_upto(av, gdiv(y, gpowgs(L, dx)));
2950 : }
2951 :
2952 : GEN
2953 476 : RgXQ_trace(GEN x, GEN T)
2954 : {
2955 476 : pari_sp av = avma;
2956 : GEN dT, z;
2957 : long n;
2958 476 : if (degpol(T)==0) return gmulgs(x,0);
2959 469 : dT = RgX_deriv(T); n = degpol(dT);
2960 469 : z = RgXQ_mul(x, dT, T);
2961 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2962 420 : return gc_upto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2963 : }
2964 :
2965 : GEN
2966 2742515 : RgX_blocks(GEN P, long n, long m)
2967 : {
2968 2742515 : GEN z = cgetg(m+1,t_VEC);
2969 2742515 : long i,j, k=2, l = lg(P);
2970 8419863 : for(i=1; i<=m; i++)
2971 : {
2972 5677348 : GEN zi = cgetg(n+2,t_POL);
2973 5677348 : zi[1] = P[1];
2974 5677348 : gel(z,i) = zi;
2975 16836166 : for(j=2; j<n+2; j++)
2976 11158818 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2977 5677348 : zi = RgX_renormalize_lg(zi, n+2);
2978 : }
2979 2742515 : return z;
2980 : }
2981 :
2982 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2983 : void
2984 49788167 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2985 : {
2986 49788167 : long n = degpol(p), v = varn(p), n0, n1, i;
2987 : GEN p0, p1;
2988 :
2989 49788580 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2990 :
2991 49695226 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2992 49695226 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2993 49697984 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2994 146869082 : for (i=0; i<n1; i++)
2995 : {
2996 97169231 : p0[2+i] = p[2+(i<<1)];
2997 97169231 : p1[2+i] = p[3+(i<<1)];
2998 : }
2999 49699851 : if (n1 != n0)
3000 18039139 : p0[2+i] = p[2+(i<<1)];
3001 49699851 : *pe = normalizepol(p0);
3002 49702258 : *po = normalizepol(p1);
3003 : }
3004 :
3005 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
3006 : GEN
3007 43421 : RgX_splitting(GEN p, long k)
3008 : {
3009 43421 : long n = degpol(p), v = varn(p), m, i, j, l;
3010 : GEN r;
3011 :
3012 43421 : m = n/k;
3013 43421 : r = cgetg(k+1,t_VEC);
3014 232799 : for(i=1; i<=k; i++)
3015 : {
3016 189378 : gel(r,i) = cgetg(m+3, t_POL);
3017 189378 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
3018 : }
3019 568491 : for (j=1, i=0, l=2; i<=n; i++)
3020 : {
3021 525070 : gmael(r,j,l) = gel(p,2+i);
3022 525070 : if (j==k) { j=1; l++; } else j++;
3023 : }
3024 232799 : for(i=1; i<=k; i++)
3025 189378 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
3026 43421 : return r;
3027 : }
3028 :
3029 : /*******************************************************************/
3030 : /* */
3031 : /* Kronecker form */
3032 : /* */
3033 : /*******************************************************************/
3034 :
3035 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
3036 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
3037 : * normalized coefficients */
3038 : GEN
3039 180567 : Kronecker_to_mod(GEN z, GEN T)
3040 : {
3041 180567 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
3042 180567 : GEN x, t = cgetg(N,t_POL);
3043 180567 : t[1] = T[1];
3044 180567 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
3045 180567 : x[1] = z[1];
3046 180567 : T = RgX_copy(T);
3047 1375406 : for (i=2; i<lx+2; i++, z+= N-2)
3048 : {
3049 5704416 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
3050 1194839 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
3051 : }
3052 180567 : N = (l-2) % (N-2) + 2;
3053 463323 : for (j=2; j<N; j++) t[j] = z[j];
3054 180567 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
3055 180567 : return normalizepol_lg(x, i+1);
3056 : }
3057 :
3058 : /*******************************************************************/
3059 : /* */
3060 : /* Domain detection */
3061 : /* */
3062 : /*******************************************************************/
3063 :
3064 : static GEN
3065 581889 : RgX_mul_FpX(GEN x, GEN y, GEN p)
3066 : {
3067 581889 : pari_sp av = avma;
3068 : GEN r;
3069 581889 : if (lgefint(p) == 3)
3070 : {
3071 533648 : ulong pp = uel(p, 2);
3072 533648 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
3073 : RgX_to_Flx(y, pp), pp));
3074 : }
3075 : else
3076 48241 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3077 581889 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3078 536123 : return gc_upto(av, FpX_to_mod(r, p));
3079 : }
3080 :
3081 : static GEN
3082 1477 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3083 : {
3084 1477 : pari_sp av = avma;
3085 : long dT;
3086 : GEN kx, ky, r;
3087 1477 : GEN T = RgX_to_FpX(pol, p);
3088 1477 : if (signe(T)==0) pari_err_OP("*", x, y);
3089 1470 : dT = degpol(T);
3090 1470 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3091 1470 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3092 1470 : r = FpX_mul(kx, ky, p);
3093 1470 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3094 1386 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3095 : }
3096 :
3097 : static GEN
3098 444268 : RgX_liftred(GEN x, GEN T)
3099 444268 : { return RgXQX_red(liftpol_shallow(x), T); }
3100 :
3101 : static GEN
3102 165574 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3103 : {
3104 165574 : pari_sp av = avma;
3105 165574 : long dT = degpol(T);
3106 165574 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3107 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3108 165574 : return gc_upto(av, Kronecker_to_mod(r, T));
3109 : }
3110 :
3111 : static GEN
3112 1323 : RgX_sqr_FpX(GEN x, GEN p)
3113 : {
3114 1323 : pari_sp av = avma;
3115 : GEN r;
3116 1323 : if (lgefint(p) == 3)
3117 : {
3118 1119 : ulong pp = uel(p, 2);
3119 1119 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3120 : }
3121 : else
3122 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3123 1323 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3124 1309 : return gc_upto(av, FpX_to_mod(r, p));
3125 : }
3126 :
3127 : static GEN
3128 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3129 : {
3130 196 : pari_sp av = avma;
3131 : long dT;
3132 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3133 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3134 189 : dT = degpol(T);
3135 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3136 189 : r = FpX_sqr(kx, p);
3137 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3138 189 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3139 : }
3140 :
3141 : static GEN
3142 13418 : RgX_sqr_QXQX(GEN x, GEN T)
3143 : {
3144 13418 : pari_sp av = avma;
3145 13418 : long dT = degpol(T);
3146 13418 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3147 13418 : return gc_upto(av, Kronecker_to_mod(r, T));
3148 : }
3149 :
3150 : static GEN
3151 67781 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3152 : {
3153 67781 : pari_sp av = avma;
3154 : GEN r;
3155 67781 : if (lgefint(p) == 3)
3156 : {
3157 51502 : ulong pp = uel(p, 2);
3158 51502 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3159 : RgX_to_Flx(y, pp), pp));
3160 : }
3161 : else
3162 16279 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3163 67781 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3164 30905 : return gc_upto(av, FpX_to_mod(r, p));
3165 : }
3166 :
3167 : static GEN
3168 49851 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3169 : {
3170 49851 : pari_sp av = avma;
3171 : GEN r;
3172 49851 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3173 49851 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3174 : }
3175 : static GEN
3176 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3177 : {
3178 70 : pari_sp av = avma;
3179 : GEN r;
3180 70 : GEN T = RgX_to_FpX(pol, p);
3181 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3182 63 : if (lgefint(p) == 3)
3183 : {
3184 55 : ulong pp = uel(p, 2);
3185 55 : GEN Tp = ZX_to_Flx(T, pp);
3186 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3187 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3188 : }
3189 : else
3190 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3191 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3192 56 : return gc_upto(av, FpXQX_to_mod(r, T, p));
3193 : }
3194 :
3195 : static GEN
3196 81760924 : RgX_mul_fast(GEN x, GEN y)
3197 : {
3198 : GEN p, pol;
3199 : long pa;
3200 81760924 : long t = RgX_type2(x,y, &p,&pol,&pa);
3201 81761070 : switch(t)
3202 : {
3203 34468244 : case t_INT: return ZX_mul(x,y);
3204 3252064 : case t_FRAC: return QX_mul(x,y);
3205 103707 : case t_FFELT: return FFX_mul(x, y, pol);
3206 581889 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3207 165574 : case RgX_type_code(t_POLMOD, t_INT):
3208 : case RgX_type_code(t_POLMOD, t_FRAC):
3209 165574 : return RgX_mul_QXQX(x, y, pol);
3210 1473 : case RgX_type_code(t_POLMOD, t_INTMOD):
3211 1473 : return RgX_mul_FpXQX(x, y, pol, p);
3212 43188119 : default: return NULL;
3213 : }
3214 : }
3215 : static GEN
3216 1305807 : RgX_sqr_fast(GEN x)
3217 : {
3218 : GEN p, pol;
3219 : long pa;
3220 1305807 : long t = RgX_type(x,&p,&pol,&pa);
3221 1305809 : switch(t)
3222 : {
3223 1172787 : case t_INT: return ZX_sqr(x);
3224 82352 : case t_FRAC: return QX_sqr(x);
3225 2499 : case t_FFELT: return FFX_sqr(x, pol);
3226 1323 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3227 13418 : case RgX_type_code(t_POLMOD, t_INT):
3228 : case RgX_type_code(t_POLMOD, t_FRAC):
3229 13418 : return RgX_sqr_QXQX(x, pol);
3230 194 : case RgX_type_code(t_POLMOD, t_INTMOD):
3231 194 : return RgX_sqr_FpXQX(x, pol, p);
3232 33236 : default: return NULL;
3233 : }
3234 : }
3235 :
3236 : static GEN
3237 14645934 : RgX_rem_fast(GEN x, GEN y)
3238 : {
3239 : GEN p, pol;
3240 : long pa;
3241 14645934 : long t = RgX_type2(x,y, &p,&pol,&pa);
3242 14646079 : switch(t)
3243 : {
3244 3453688 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3245 1836607 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3246 84 : case t_FFELT: return FFX_rem(x, y, pol);
3247 67781 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3248 49851 : case RgX_type_code(t_POLMOD, t_INT):
3249 : case RgX_type_code(t_POLMOD, t_FRAC):
3250 49851 : return RgX_rem_QXQX(x, y, pol);
3251 65 : case RgX_type_code(t_POLMOD, t_INTMOD):
3252 65 : return RgX_rem_FpXQX(x, y, pol, p);
3253 9238003 : default: return NULL;
3254 : }
3255 : }
3256 :
3257 : GEN
3258 68212064 : RgX_mul(GEN x, GEN y)
3259 : {
3260 68212064 : GEN z = RgX_mul_fast(x,y);
3261 68212151 : if (!z) z = RgX_mul_i(x,y);
3262 68211842 : return z;
3263 : }
3264 :
3265 : GEN
3266 1293655 : RgX_sqr(GEN x)
3267 : {
3268 1293655 : GEN z = RgX_sqr_fast(x);
3269 1293653 : if (!z) z = RgX_sqr_i(x);
3270 1293653 : return z;
3271 : }
3272 :
3273 : GEN
3274 14645973 : RgX_rem(GEN x, GEN y)
3275 : {
3276 14645973 : GEN z = RgX_rem_fast(x, y);
3277 14646063 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3278 14645992 : return z;
3279 : }
3280 :
3281 : static GEN
3282 0 : _RgX_mul(void* E, GEN x, GEN y)
3283 0 : { (void) E; return RgX_mul(x, y); }
3284 :
3285 : GEN
3286 0 : RgXV_prod(GEN V)
3287 0 : { return gen_product(V, NULL, &_RgX_mul); }
3288 :
3289 : GEN
3290 10922720 : RgXn_mul(GEN f, GEN g, long n)
3291 : {
3292 10922720 : pari_sp av = avma;
3293 10922720 : GEN h = RgX_mul_fast(f,g);
3294 10922719 : if (!h) return RgXn_mul2(f,g,n);
3295 1679095 : if (degpol(h) < n) return h;
3296 384712 : return gc_GEN(av, RgXn_red_shallow(h, n));
3297 : }
3298 :
3299 : GEN
3300 12152 : RgXn_sqr(GEN f, long n)
3301 : {
3302 12152 : pari_sp av = avma;
3303 12152 : GEN g = RgX_sqr_fast(f);
3304 12152 : if (!g) return RgXn_sqr2(f,n);
3305 11837 : if (degpol(g) < n) return g;
3306 10640 : return gc_GEN(av, RgXn_red_shallow(g, n));
3307 : }
3308 :
3309 : /* (f*g) \/ x^n */
3310 : GEN
3311 2626137 : RgX_mulhigh_i(GEN f, GEN g, long n)
3312 : {
3313 2626137 : GEN h = RgX_mul_fast(f,g);
3314 2626137 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3315 : }
3316 :
3317 : /* (f*g) \/ x^n */
3318 : GEN
3319 0 : RgX_sqrhigh_i(GEN f, long n)
3320 : {
3321 0 : GEN h = RgX_sqr_fast(f);
3322 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3323 : }
|
