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 19605536 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19605536 : long pold=1, rold=n*(d-1);
35 91813353 : for(p=2; p<=d; p++)
36 : {
37 72207817 : r = m*(p-1) + n*((d-1)/p);
38 72207817 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19605536 : return pold;
41 : }
42 :
43 : static GEN
44 13931238 : 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 13931238 : pari_sp av = avma;
48 : long i;
49 13931238 : GEN z = cmul(E,P,a,ff->one(E));
50 13909815 : if (!z) z = gen_0;
51 73116518 : for (i=1; i<=n; i++)
52 : {
53 59207816 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 59231776 : if (t) {
55 44310630 : z = ff->add(E, z, t);
56 44272230 : if (gc_needed(av,2)) z = gc_upto(av, z);
57 : }
58 : }
59 13908702 : 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 12131091 : 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 12131091 : pari_sp av = avma;
72 12131091 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12131091 : if (d < 0) return ff->zero(E);
76 11521491 : if (d < l) return gc_upto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 994277 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 994277 : 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 994277 : d -= l;
84 994277 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2409186 : while (d >= l-1)
86 : {
87 1413925 : d -= l-1;
88 1413925 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1414095 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1413851 : if (gc_needed(av,2))
91 91 : z = gc_upto(av, z);
92 : }
93 995261 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 995291 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 995266 : return gc_upto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 869279 : 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 869279 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 869279 : if (d < 0) return ff->zero(E);
106 868005 : rtd = (long) sqrt((double)d);
107 868005 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 868036 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 868019 : return gc_upto(av, z);
110 : }
111 :
112 : static GEN
113 2034097 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19258783 : _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 1845501 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 593606 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2075652 : _gen_one(void *E) { (void)E; return gen_1; }
124 : static GEN
125 24184 : _gen_zero(void *E) { (void)E; return gen_0; }
126 :
127 : static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
128 : _gen_mul, _gen_sqr,_gen_one,_gen_zero };
129 :
130 : static GEN
131 512818 : _gen_cmul(void *E, GEN P, long a, GEN x)
132 512818 : {(void)E; return gmul(gel(P,a+2), x);}
133 :
134 : GEN
135 166768 : RgX_RgV_eval(GEN Q, GEN x)
136 : {
137 166768 : return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
138 : }
139 :
140 : GEN
141 0 : RgX_Rg_eval_bk(GEN Q, GEN x)
142 : {
143 0 : return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
144 : }
145 :
146 : GEN
147 2947 : RgXV_RgV_eval(GEN Q, GEN x)
148 : {
149 2947 : long i, l = lg(Q), vQ = gvar(Q);
150 2947 : GEN v = cgetg(l, t_VEC);
151 248311 : for (i = 1; i < l; i++)
152 : {
153 245364 : GEN Qi = gel(Q, i);
154 245364 : gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
155 : }
156 2947 : return v;
157 : }
158 :
159 : GEN
160 398923 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 398923 : pari_sp av = avma, btop;
163 398923 : long i, d = degpol(P), o;
164 : GEN s;
165 398922 : if (signe(P)==0) return pol_0(varn(P));
166 398922 : s = gel(P, d+2);
167 398922 : if (d == 0) return gcopy(s);
168 395611 : o = RgX_deflate_order(P);
169 395621 : if (o > 1) A = gpowgs(A, o);
170 395620 : btop = avma;
171 1510161 : for (i = d-o; i >= 0; i-=o)
172 : {
173 1114548 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
174 1114541 : 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 395613 : 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 282510 : get_Rg_algebra(void)
207 : {
208 282510 : return &Rg_algebra;
209 : }
210 :
211 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
212 :
213 : static GEN
214 11333 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
215 : {
216 : (void) E;
217 11333 : return RgX_divrem(x, y, r);
218 : }
219 :
220 : GEN
221 3157 : RgX_digits(GEN x, GEN T)
222 : {
223 3157 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
224 3157 : if (signe(x)==0) return(cgetg(1, t_VEC));
225 3157 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
226 : }
227 :
228 : /*******************************************************************/
229 : /* */
230 : /* RgX */
231 : /* */
232 : /*******************************************************************/
233 :
234 : long
235 24796888 : RgX_equal(GEN x, GEN y)
236 : {
237 24796888 : long i = lg(x);
238 :
239 24796888 : if (i != lg(y)) return 0;
240 106679186 : for (i--; i > 1; i--)
241 82210479 : if (!gequal(gel(x,i),gel(y,i))) return 0;
242 24468707 : 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 156956984 : Rg_get_1(GEN x)
249 : {
250 : GEN p, T;
251 156956984 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
252 156956984 : if (RgX_type_is_composite(tx))
253 11657742 : RgX_type_decode(tx, &i /*junk*/, &tx);
254 156956984 : switch(tx)
255 : {
256 427406 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
257 945 : case t_PADIC: return cvtop(gen_1, p, lx);
258 5950 : case t_FFELT: return FF_1(T);
259 156522683 : 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 6942982 : Rg_get_0(GEN x)
266 : {
267 : GEN p, T;
268 6942982 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
269 6942982 : if (RgX_type_is_composite(tx))
270 61754 : RgX_type_decode(tx, &i /*junk*/, &tx);
271 6942982 : switch(tx)
272 : {
273 497 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
274 42 : case t_PADIC: return zeropadic(p, lx);
275 210 : case t_FFELT: return FF_zero(T);
276 6942233 : default: return gen_0;
277 : }
278 : }
279 :
280 : GEN
281 7595 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
282 : {
283 7595 : long i, n = degpol(P);
284 : GEN z, dz, dP;
285 7595 : if (n < 0) return gen_0;
286 7595 : P = Q_remove_denom(P, &dP);
287 7595 : z = gel(P,2); if (n == 0) return icopy(z);
288 4361 : if (dV) z = mulii(dV, z); /* V[1] = dV */
289 4361 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
290 8169 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
291 4361 : dz = mul_denom(dP, dV);
292 4361 : return dz? RgX_Rg_div(z, dz): z;
293 : }
294 :
295 : /* Return P(h * x), not memory clean */
296 : GEN
297 51885 : RgX_unscale(GEN P, GEN h)
298 : {
299 51885 : long i, l = lg(P);
300 51885 : GEN hi = gen_1, Q = cgetg(l, t_POL);
301 51885 : Q[1] = P[1];
302 51885 : if (l == 2) return Q;
303 40643 : gel(Q,2) = gcopy(gel(P,2));
304 95894 : for (i=3; i<l; i++)
305 : {
306 55251 : hi = gmul(hi,h);
307 55251 : gel(Q,i) = gmul(gel(P,i), hi);
308 : }
309 40643 : return Q;
310 : }
311 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 and 2^k */
312 : GEN
313 1171613 : ZX_z_unscale(GEN P, long h)
314 : {
315 1171613 : long i, l = lg(P);
316 1171613 : GEN Q = cgetg(l, t_POL);
317 1171613 : Q[1] = P[1];
318 1171613 : if (l == 2) return Q;
319 1128745 : gel(Q,2) = gel(P,2);
320 1128745 : if (l == 3) return Q;
321 1108165 : if (h == -1)
322 241717 : for (i = 3; i < l; i++)
323 : {
324 199006 : gel(Q,i) = negi(gel(P,i));
325 199006 : if (++i == l) break;
326 147988 : gel(Q,i) = gel(P,i);
327 : }
328 1014436 : else if (h > 0 && !(h & (h-1))) return ZX_unscale2n(P, vals(h));
329 : else
330 : {
331 : GEN hi;
332 492748 : gel(Q,3) = mulis(gel(P,3), h);
333 492747 : hi = sqrs(h);
334 2456515 : for (i = 4; i < l; i++)
335 : {
336 1963768 : gel(Q,i) = mulii(gel(P,i), hi);
337 1963769 : if (i != l-1) hi = mulis(hi,h);
338 : }
339 : }
340 586476 : return Q;
341 : }
342 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean */
343 : GEN
344 674279 : ZX_unscale(GEN P, GEN h)
345 : {
346 : long i, l;
347 : GEN Q, hi;
348 674279 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
349 881 : l = lg(P); Q = cgetg(l, t_POL);
350 881 : Q[1] = P[1];
351 881 : if (l == 2) return Q;
352 881 : gel(Q,2) = gel(P,2);
353 881 : if (l == 3) return Q;
354 881 : hi = h;
355 881 : gel(Q,3) = mulii(gel(P,3), hi);
356 2741 : for (i = 4; i < l; i++)
357 : {
358 1860 : hi = mulii(hi,h);
359 1860 : gel(Q,i) = mulii(gel(P,i), hi);
360 : }
361 881 : return Q;
362 : }
363 : /* P a ZX. Return P(x << n), not memory clean */
364 : GEN
365 1974124 : ZX_unscale2n(GEN P, long n)
366 : {
367 1974124 : long i, ni = n, l = lg(P);
368 1974124 : GEN Q = cgetg(l, t_POL);
369 1974124 : Q[1] = P[1];
370 1974124 : if (l == 2) return Q;
371 1928446 : gel(Q,2) = gel(P,2);
372 1928446 : if (l == 3) return Q;
373 1922720 : gel(Q,3) = shifti(gel(P,3), ni);
374 6721559 : for (i=4; i<l; i++)
375 : {
376 4798891 : ni += n;
377 4798891 : gel(Q,i) = shifti(gel(P,i), ni);
378 : }
379 1922668 : 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 6558 : RgXV_unscale(GEN x, GEN h)
420 : {
421 6558 : if (isint1(h)) return gcopy(x);
422 18733 : 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 4344478 : RgX_rescale(GEN P, GEN h)
428 : {
429 4344478 : long i, l = lg(P);
430 4344478 : GEN Q = cgetg(l,t_POL), hi = h;
431 4344478 : gel(Q,l-1) = gel(P,l-1);
432 11432390 : for (i=l-2; i>=2; i--)
433 : {
434 11429769 : gel(Q,i) = gmul(gel(P,i), hi);
435 11429665 : if (i == 2) break;
436 7087672 : hi = gmul(hi,h);
437 : }
438 4344614 : 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 1187896 : RgX_deflate(GEN x0, long d)
451 : {
452 : GEN z, y, x;
453 1187896 : long i,id, dy, dx = degpol(x0);
454 1187891 : if (d == 1 || dx <= 0) return x0;
455 463931 : dy = dx/d;
456 463931 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
457 463932 : z = y + 2;
458 463932 : x = x0+ 2;
459 1778338 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
460 463932 : 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 18732 : RgX_homogenize(GEN P, long v)
475 : {
476 : long i, l, d;
477 18732 : GEN Q = cgetg_copy(P, &l);
478 18732 : Q[1] = P[1]; d = l-3;
479 165998 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
480 18732 : 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 885486 : RgX_inflate(GEN x0, long d)
514 : {
515 885486 : long i, id, dy, dx = degpol(x0);
516 885486 : GEN x = x0 + 2, z, y;
517 885486 : if (dx <= 0) return leafcopy(x0);
518 808196 : dy = dx*d;
519 808196 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
520 808196 : z = y + 2;
521 28121538 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
522 11507457 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
523 808196 : return y;
524 : }
525 :
526 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
527 : static GEN
528 5703072 : RgX_Rg_translate_basecase(GEN P, GEN c)
529 : {
530 5703072 : pari_sp av = avma;
531 : GEN Q;
532 : long i, k, n;
533 :
534 5703072 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
535 5701145 : Q = leafcopy(P); n = degpol(P);
536 5701146 : if (isint1(c))
537 : {
538 7595 : 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 5698935 : else if (isintm1(c))
549 : {
550 15974 : 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 19429122 : for (i=1; i<=n; i++)
563 : {
564 45458865 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1)));
565 13733858 : 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 5700896 : return gc_GEN(av, Q);
573 : }
574 :
575 : static GEN
576 102627 : zero_FpX_mod(GEN p, long v)
577 : {
578 102627 : GEN r = cgetg(3,t_POL);
579 102627 : r[1] = evalvarn(v);
580 102627 : gel(r,2) = mkintmod(gen_0, icopy(p));
581 102627 : return r;
582 : }
583 :
584 : static GEN
585 616 : zero_FpXQX_mod(GEN pol, GEN p, long v)
586 : {
587 616 : GEN r = cgetg(3,t_POL);
588 616 : r[1] = evalvarn(v);
589 616 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
590 616 : 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 6064818 : RgX_Rg_translate_fast(GEN P, GEN c)
625 : {
626 : GEN p, pol;
627 : long pa;
628 6064818 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
629 6064811 : switch(t)
630 : {
631 361925 : 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 5702881 : default: return NULL;
636 : }
637 : }
638 :
639 : static GEN
640 5703259 : RgX_Rg_translate_i(GEN P, GEN c)
641 : {
642 5703259 : pari_sp av = avma;
643 : long n;
644 5703259 : n = degpol(P);
645 5703260 : if (n < 40)
646 5703071 : 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 6064819 : RgX_Rg_translate(GEN P, GEN c)
659 : {
660 6064819 : GEN R = RgX_Rg_translate_fast(P, c);
661 6064810 : return R ? R: RgX_Rg_translate_i(P,c);
662 : }
663 : /* P(ax + b) */
664 : GEN
665 30212 : RgX_affine(GEN P, GEN a, GEN b)
666 : {
667 30212 : if (!gequal0(b)) P = RgX_Rg_translate(P, b);
668 30212 : 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 169946 : QXQ_to_mod(GEN x, GEN T)
708 : {
709 : long d;
710 169946 : switch(typ(x))
711 : {
712 64796 : case t_INT: return icopy(x);
713 2079 : case t_FRAC: return gcopy(x);
714 103071 : case t_POL:
715 103071 : d = degpol(x);
716 103071 : if (d < 0) return gen_0;
717 100810 : if (d == 0) return gcopy(gel(x,2));
718 98671 : 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 544623 : case t_INT:
731 544623 : case t_FRAC: return x;
732 296072 : case t_POL:
733 296072 : d = degpol(x);
734 296072 : if (d < 0) return gen_0;
735 249856 : 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 37471 : QXQX_to_mod(GEN z, GEN T)
745 : {
746 37471 : long i,l = lg(z);
747 37471 : GEN x = cgetg(l,t_POL);
748 187922 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
749 37471 : 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 12733 : QXQXV_to_mod(GEN V, GEN T)
763 : {
764 12733 : long i, l = lg(V);
765 12733 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
766 50204 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
767 12733 : return z;
768 : }
769 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
770 : GEN
771 22338 : QXQV_to_mod(GEN V, GEN T)
772 : {
773 22338 : long i, l = lg(V);
774 22338 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
775 41833 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
776 22338 : 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 7886648 : RgX_renormalize_lg(GEN x, long lx)
800 : {
801 : long i;
802 10938497 : for (i = lx-1; i>1; i--)
803 10485235 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
804 7886648 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
805 7886647 : setlg(x, i+1); setsigne(x, i != 1); return x;
806 : }
807 :
808 : GEN
809 2289089 : RgV_to_RgX(GEN x, long v)
810 : {
811 2289089 : long i, k = lg(x);
812 : GEN p;
813 :
814 6366247 : while (--k && gequal0(gel(x,k)));
815 2289088 : if (!k) return pol_0(v);
816 2255948 : i = k+2; p = cgetg(i,t_POL);
817 2255950 : p[1] = evalsigne(1) | evalvarn(v);
818 24359937 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
819 2255950 : 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 15066043 : RgX_to_RgC(GEN x, long N)
842 : {
843 : long i, l;
844 : GEN z;
845 15066043 : l = lg(x)-1; x++;
846 15066043 : if (l > N+1) l = N+1; /* truncate higher degree terms */
847 15066043 : z = cgetg(N+1,t_COL);
848 110956258 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
849 31161061 : for ( ; i<=N; i++) gel(z,i) = gen_0;
850 15066200 : return z;
851 : }
852 : GEN
853 1360835 : Rg_to_RgC(GEN x, long N)
854 : {
855 1360835 : 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 303243 : RgM_to_RgXV(GEN x, long v)
861 1290979 : { 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 336712 : RgV_to_RgM(GEN x, long n)
870 1692348 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
871 :
872 : GEN
873 79032 : RgXV_to_RgM(GEN x, long n)
874 401865 : { 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 23686 : RgM_to_RgXX(GEN x, long v,long w)
879 : {
880 23686 : long j, lx = lg(x);
881 23686 : GEN y = cgetg(lx+1, t_POL);
882 23686 : y[1] = evalsigne(1) | evalvarn(v);
883 23686 : y++;
884 131945 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
885 23686 : 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 32771 : RgXY_swapspec(GEN x, long n, long w, long nx)
902 : {
903 32771 : long j, ly = n+3;
904 32771 : GEN y = cgetg(ly, t_POL);
905 32771 : y[1] = evalsigne(1);
906 401037 : for (j=2; j<ly; j++)
907 : {
908 : long k;
909 368266 : GEN a = cgetg(nx+2,t_POL);
910 368266 : a[1] = evalsigne(1) | evalvarn(w);
911 2028298 : for (k=0; k<nx; k++)
912 : {
913 1660032 : GEN xk = gel(x,k);
914 1660032 : if (typ(xk)==t_POL && varn(xk)==w)
915 1565657 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
916 : else
917 94375 : gel(a,k+2) = j==2 ? xk: gen_0;
918 : }
919 368266 : gel(y,j) = normalizepol_lg(a, nx+2);
920 : }
921 32771 : 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 2022702 : RgXY_degreex(GEN b)
934 : {
935 2022702 : long deg = 0, i;
936 2022702 : if (!signe(b)) return -1;
937 8498253 : for (i = 2; i < lg(b); ++i)
938 : {
939 6475550 : GEN bi = gel(b, i);
940 6475550 : if (typ(bi) == t_POL)
941 1024519 : deg = maxss(deg, degpol(bi));
942 : }
943 2022703 : 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 7944198 : RgXn_red_shallow(GEN a, long n)
952 : {
953 7944198 : long i, L = n+2, l = lg(a);
954 : GEN b;
955 7944198 : if (L >= l) return a; /* deg(x) < n */
956 5751840 : b = cgetg(L, t_POL); b[1] = a[1];
957 37334921 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
958 5751840 : 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 177403083 : RgX_shift_shallow(GEN a, long n)
968 : {
969 177403083 : long i, l = lg(a);
970 : GEN b;
971 177403083 : if (l == 2 || !n) return a;
972 109225704 : l += n;
973 109225704 : if (n < 0)
974 : {
975 54215560 : if (l <= 2) return pol_0(varn(a));
976 52741449 : b = cgetg(l, t_POL); b[1] = a[1];
977 52742224 : a -= n;
978 164024652 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
979 : } else {
980 55010144 : b = cgetg(l, t_POL); b[1] = a[1];
981 55014028 : a -= n; n += 2;
982 119567985 : for (i=2; i<n; i++) gel(b,i) = gen_0;
983 212280007 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
984 : }
985 107756252 : return b;
986 : }
987 : /* return (x * X^n). */
988 : GEN
989 1578036 : RgX_shift(GEN a, long n)
990 : {
991 1578036 : long i, l = lg(a);
992 : GEN b;
993 1578036 : if (l == 2 || !n) return RgX_copy(a);
994 1577084 : l += n;
995 1577084 : 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 1575642 : b = cgetg(l, t_POL); b[1] = a[1];
1003 1575642 : a -= n; n += 2;
1004 3979090 : for (i=2; i<n; i++) gel(b,i) = gen_0;
1005 4346111 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1006 : }
1007 1577014 : return b;
1008 : }
1009 :
1010 : GEN
1011 317037 : RgX_rotate_shallow(GEN P, long k, long p)
1012 : {
1013 317037 : long i, l = lgpol(P);
1014 : GEN r;
1015 317037 : if (signe(P)==0)
1016 1365 : return pol_0(varn(P));
1017 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
1018 2100644 : for(i=0; i<p; i++)
1019 : {
1020 1784972 : long s = 2+(i+k)%p;
1021 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
1022 : }
1023 315672 : return RgX_renormalize(r);
1024 : }
1025 :
1026 : GEN
1027 2997279 : RgX_mulXn(GEN x, long d)
1028 : {
1029 : pari_sp av;
1030 : GEN z;
1031 : long v;
1032 2997279 : if (d >= 0) return RgX_shift(x, d);
1033 1472517 : d = -d;
1034 1472517 : v = RgX_val(x);
1035 1472517 : if (v >= d) return RgX_shift(x, -d);
1036 1472503 : av = avma;
1037 1472503 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1038 1472503 : 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 3667912 : RgX_val(GEN x)
1052 : {
1053 3667912 : long i, lx = lg(x);
1054 3667912 : if (lx == 2) return LONG_MAX;
1055 4619751 : for (i = 2; i < lx; i++)
1056 4619695 : if (!isexactzero(gel(x,i))) break;
1057 3656985 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1058 3656929 : return i - 2;
1059 : }
1060 : long
1061 80759632 : RgX_valrem(GEN x, GEN *Z)
1062 : {
1063 80759632 : long v, i, lx = lg(x);
1064 80759632 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1065 125176791 : for (i = 2; i < lx; i++)
1066 125177750 : if (!isexactzero(gel(x,i))) break;
1067 : /* possible with nonrational zeros */
1068 80759209 : if (i == lx)
1069 : {
1070 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1071 14 : return LONG_MAX;
1072 : }
1073 80759195 : v = i - 2;
1074 80759195 : *Z = RgX_shift_shallow(x, -v);
1075 80770196 : return v;
1076 : }
1077 : long
1078 852357 : RgX_valrem_inexact(GEN x, GEN *Z)
1079 : {
1080 : long v;
1081 852357 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1082 869344 : for (v = 0;; v++)
1083 869344 : if (!gequal0(gel(x,2+v))) break;
1084 852350 : if (Z) *Z = RgX_shift_shallow(x, -v);
1085 852350 : 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 1414 : RgXQV_red(GEN x, GEN T)
1094 34937 : { 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 508139 : RgXQX_red(GEN P, GEN T)
1108 : {
1109 508139 : long i, l = lg(P);
1110 508139 : GEN Q = cgetg(l, t_POL);
1111 508139 : Q[1] = P[1];
1112 2625879 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1113 508139 : return normalizepol_lg(Q, l);
1114 : }
1115 :
1116 : GEN
1117 830930 : RgX_deriv(GEN x)
1118 : {
1119 830930 : long i,lx = lg(x)-1;
1120 : GEN y;
1121 :
1122 830930 : if (lx<3) return pol_0(varn(x));
1123 827815 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1124 3653516 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1125 827801 : y[1] = x[1]; return normalizepol_lg(y,i);
1126 : }
1127 :
1128 : GEN
1129 2584178 : RgX_recipspec_shallow(GEN x, long l, long n)
1130 : {
1131 : long i;
1132 2584178 : GEN z = cgetg(n+2,t_POL);
1133 2584186 : z[1] = 0; z += 2;
1134 144346777 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1135 2829525 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1136 2584186 : return normalizepol_lg(z-2,n+2);
1137 : }
1138 :
1139 : GEN
1140 643038 : RgXn_recip_shallow(GEN P, long n)
1141 : {
1142 643038 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1143 643047 : setvarn(Q, varn(P));
1144 643047 : return Q;
1145 : }
1146 :
1147 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1148 : GEN
1149 31388 : RgX_recip(GEN x)
1150 : {
1151 : long lx, i, j;
1152 31388 : GEN y = cgetg_copy(x, &lx);
1153 274519 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1154 31388 : 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 3751805 : RgX_recip_i(GEN x)
1168 : {
1169 : long lx, i, j;
1170 3751805 : GEN y = cgetg_copy(x, &lx);
1171 20233614 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1172 3751801 : return y;
1173 : }
1174 : /*******************************************************************/
1175 : /* */
1176 : /* ADDITION / SUBTRACTION */
1177 : /* */
1178 : /*******************************************************************/
1179 : /* same variable */
1180 : GEN
1181 143425113 : RgX_add(GEN x, GEN y)
1182 : {
1183 143425113 : long i, lx = lg(x), ly = lg(y);
1184 : GEN z;
1185 143425113 : if (ly <= lx) {
1186 129244916 : z = cgetg(lx,t_POL); z[1] = x[1];
1187 506170973 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1188 188909549 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1189 129220633 : z = normalizepol_lg(z, lx);
1190 : } else {
1191 14180197 : z = cgetg(ly,t_POL); z[1] = y[1];
1192 53723231 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1193 41722727 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1194 14181860 : z = normalizepol_lg(z, ly);
1195 : }
1196 143418056 : return z;
1197 : }
1198 : GEN
1199 72856240 : RgX_sub(GEN x, GEN y)
1200 : {
1201 72856240 : long i, lx = lg(x), ly = lg(y);
1202 : GEN z;
1203 72856240 : if (ly <= lx) {
1204 39005705 : z = cgetg(lx,t_POL); z[1] = x[1];
1205 188430785 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1206 66505674 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1207 38982963 : z = normalizepol_lg(z, lx);
1208 : } else {
1209 33850535 : z = cgetg(ly,t_POL); z[1] = y[1];
1210 125780975 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1211 74979692 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1212 33818068 : z = normalizepol_lg(z, ly);
1213 : }
1214 72851205 : return z;
1215 : }
1216 : GEN
1217 7451306 : RgX_neg(GEN x)
1218 48570530 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1219 :
1220 : GEN
1221 42746530 : RgX_Rg_add(GEN y, GEN x)
1222 : {
1223 : GEN z;
1224 42746530 : long lz = lg(y), i;
1225 42746530 : if (lz == 2) return scalarpol(x,varn(y));
1226 20997459 : z = cgetg(lz,t_POL); z[1] = y[1];
1227 20997432 : gel(z,2) = gadd(gel(y,2),x);
1228 75119097 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1229 : /* probably useless unless lz = 3, but cannot be skipped if y is
1230 : * an inexact 0 */
1231 20997394 : return normalizepol_lg(z,lz);
1232 : }
1233 : GEN
1234 65102 : RgX_Rg_add_shallow(GEN y, GEN x)
1235 : {
1236 : GEN z;
1237 65102 : long lz = lg(y), i;
1238 65102 : if (lz == 2) return scalarpol(x,varn(y));
1239 65102 : z = cgetg(lz,t_POL); z[1] = y[1];
1240 65102 : gel(z,2) = gadd(gel(y,2),x);
1241 130344 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1242 65102 : return normalizepol_lg(z,lz);
1243 : }
1244 : GEN
1245 1226082 : RgX_Rg_sub(GEN y, GEN x)
1246 : {
1247 : GEN z;
1248 1226082 : long lz = lg(y), i;
1249 1226082 : if (lz == 2)
1250 : { /* scalarpol(gneg(x),varn(y)) optimized */
1251 922867 : long v = varn(y);
1252 922867 : if (isrationalzero(x)) return pol_0(v);
1253 2235 : z = cgetg(3,t_POL);
1254 2235 : z[1] = gequal0(x)? evalvarn(v)
1255 2235 : : evalvarn(v) | evalsigne(1);
1256 2235 : gel(z,2) = gneg(x); return z;
1257 : }
1258 303215 : z = cgetg(lz,t_POL); z[1] = y[1];
1259 303215 : gel(z,2) = gsub(gel(y,2),x);
1260 697370 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1261 303216 : return normalizepol_lg(z,lz);
1262 : }
1263 : GEN
1264 4954166 : Rg_RgX_sub(GEN x, GEN y)
1265 : {
1266 : GEN z;
1267 4954166 : long lz = lg(y), i;
1268 4954166 : if (lz == 2) return scalarpol(x,varn(y));
1269 4883772 : z = cgetg(lz,t_POL); z[1] = y[1];
1270 4883757 : gel(z,2) = gsub(x, gel(y,2));
1271 7337288 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1272 4883762 : return normalizepol_lg(z,lz);
1273 : }
1274 : /*******************************************************************/
1275 : /* */
1276 : /* KARATSUBA MULTIPLICATION */
1277 : /* */
1278 : /*******************************************************************/
1279 : #if 0
1280 : /* to debug Karatsuba-like routines */
1281 : GEN
1282 : zx_debug_spec(GEN x, long nx)
1283 : {
1284 : GEN z = cgetg(nx+2,t_POL);
1285 : long i;
1286 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1287 : z[1] = evalsigne(1); return z;
1288 : }
1289 :
1290 : GEN
1291 : RgX_debug_spec(GEN x, long nx)
1292 : {
1293 : GEN z = cgetg(nx+2,t_POL);
1294 : long i;
1295 : for (i=0; i<nx; i++) z[i+2] = x[i];
1296 : z[1] = evalsigne(1); return z;
1297 : }
1298 : #endif
1299 :
1300 : /* generic multiplication */
1301 : GEN
1302 8851433 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1303 : {
1304 : GEN z, t;
1305 : long i;
1306 8851433 : if (nx == ny) {
1307 1478651 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1308 4795888 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1309 1478642 : return normalizepol_lg(z, nx+2);
1310 : }
1311 7372782 : if (ny < nx) {
1312 7187487 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1313 26455612 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1314 17174557 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1315 7186985 : return normalizepol_lg(z, nx+2);
1316 : } else {
1317 185295 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1318 3672372 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1319 447507 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1320 185310 : return normalizepol_lg(z, ny+2);
1321 : }
1322 : }
1323 : GEN
1324 222394 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1325 : {
1326 : GEN z, t;
1327 : long i;
1328 222394 : if (nx == ny) {
1329 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1330 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1331 12824 : return normalizepol_lg(z, nx+2);
1332 : }
1333 209570 : if (ny < nx) {
1334 207785 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1335 3725192 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1336 2372263 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1337 207785 : return normalizepol_lg(z, nx+2);
1338 : } else {
1339 1785 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1340 331478 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1341 12236 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1342 1785 : return normalizepol_lg(z, ny+2);
1343 : }
1344 : }
1345 :
1346 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1347 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1348 : * t_VECSMALL, to hold this, which can be left on stack: GC functions
1349 : * will not crash on it. The returned vector itself is not a proper GEN,
1350 : * we access the coefficients as x[i], i = 0..deg(x) */
1351 : static GEN
1352 81718700 : RgXspec_kill0(GEN x, long lx)
1353 : {
1354 81718700 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit GC-wise */
1355 : long i;
1356 246990245 : for (i=0; i <lx; i++)
1357 : {
1358 165271526 : GEN c = gel(x,i);
1359 165271526 : z[i] = (long)(isrationalzero(c)? NULL: c);
1360 : }
1361 81718719 : return z;
1362 : }
1363 :
1364 : INLINE GEN
1365 111263752 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1366 : {
1367 111263752 : pari_sp av = avma;
1368 111263752 : GEN s = NULL;
1369 : long i;
1370 :
1371 394990174 : for (i=a; i<b; i++)
1372 283735497 : if (gel(y,i) && gel(x,-i))
1373 : {
1374 205204001 : GEN t = gmul(gel(y,i), gel(x,-i));
1375 205199010 : s = s? gadd(s, t): t;
1376 : }
1377 111254677 : return s? gc_upto(av, s): gen_0;
1378 : }
1379 :
1380 : /* assume nx >= ny > 0, return x * y * t^v */
1381 : static GEN
1382 34206985 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1383 : {
1384 : long i, lz, nz;
1385 : GEN z;
1386 :
1387 34206985 : x = RgXspec_kill0(x,nx);
1388 34206965 : y = RgXspec_kill0(y,ny);
1389 34206961 : lz = nx + ny + 1; nz = lz-2;
1390 34206961 : lz += v;
1391 34206961 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1392 72784140 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1393 83211913 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1394 56646131 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1395 49005076 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1396 34206452 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1397 : }
1398 :
1399 : /* return (x * X^d) + y. Assume d > 0 */
1400 : GEN
1401 9527057 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1402 : {
1403 : GEN x, y, xd, yd, zd;
1404 : long a, lz, nx, ny;
1405 :
1406 9527057 : if (!signe(x0)) return y0;
1407 9373021 : ny = lgpol(y0);
1408 9373014 : nx = lgpol(x0);
1409 9373085 : zd = (GEN)avma;
1410 9373085 : x = x0 + 2; y = y0 + 2; a = ny-d;
1411 9373085 : if (a <= 0)
1412 : {
1413 1477551 : lz = nx+d+2;
1414 1477551 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1415 3394245 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1416 1477553 : x = zd + a;
1417 1495250 : while (zd > x) gel(--zd,0) = gen_0;
1418 : }
1419 : else
1420 : {
1421 7895534 : xd = new_chunk(d); yd = y+d;
1422 7895541 : x = RgX_addspec_shallow(x,yd, nx,a);
1423 7895463 : lz = (a>nx)? ny+2: lg(x)+d;
1424 39471890 : x += 2; while (xd > x) *--zd = *--xd;
1425 : }
1426 22857803 : while (yd > y) *--zd = *--yd;
1427 9373016 : *--zd = x0[1];
1428 9373016 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1429 : }
1430 : GEN
1431 514864 : RgX_addmulXn(GEN x0, GEN y0, long d)
1432 : {
1433 : GEN x, y, xd, yd, zd;
1434 : long a, lz, nx, ny;
1435 :
1436 514864 : if (!signe(x0)) return RgX_copy(y0);
1437 514080 : nx = lgpol(x0);
1438 514080 : ny = lgpol(y0);
1439 514080 : zd = (GEN)avma;
1440 514080 : x = x0 + 2; y = y0 + 2; a = ny-d;
1441 514080 : if (a <= 0)
1442 : {
1443 291686 : lz = nx+d+2;
1444 291686 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1445 4293251 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1446 291686 : x = zd + a;
1447 757796 : while (zd > x) gel(--zd,0) = gen_0;
1448 : }
1449 : else
1450 : {
1451 222394 : xd = new_chunk(d); yd = y+d;
1452 222394 : x = RgX_addspec(x,yd, nx,a);
1453 222394 : lz = (a>nx)? ny+2: lg(x)+d;
1454 8417378 : x += 2; while (xd > x) *--zd = *--xd;
1455 : }
1456 2615648 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1457 514080 : *--zd = x0[1];
1458 514080 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1459 : }
1460 :
1461 : /* return x * y mod t^n */
1462 : static GEN
1463 6633454 : RgXn_mul_basecase(GEN x, GEN y, long n)
1464 : {
1465 6633454 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1466 : GEN z;
1467 6633454 : if (lx < 0) return pol_0(varn(x));
1468 6633454 : if (ly < 0) return pol_0(varn(x));
1469 6633454 : z = cgetg(lz, t_POL) + 2;
1470 6633454 : x+=2; if (lx > n) lx = n;
1471 6633454 : y+=2; if (ly > n) ly = n;
1472 6633454 : z[-1] = x[-1];
1473 6633454 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1474 6633454 : x = RgXspec_kill0(x, lx);
1475 6633454 : y = RgXspec_kill0(y, ly);
1476 : /* x:y:z [i] = term of degree i */
1477 26409798 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1478 11832239 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1479 6679526 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1480 6633454 : return normalizepol_lg(z - 2, lz);
1481 : }
1482 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1483 : static GEN
1484 9420016 : RgXn_mul2(GEN f, GEN g, long n)
1485 : {
1486 9420016 : pari_sp av = avma;
1487 : GEN fe,fo, ge,go, l,h,m;
1488 : long n0, n1;
1489 9420016 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1490 6668951 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1491 35497 : n0 = n>>1; n1 = n-n0;
1492 35497 : RgX_even_odd(f, &fe, &fo);
1493 35497 : RgX_even_odd(g, &ge, &go);
1494 35497 : l = RgXn_mul2(fe,ge,n1);
1495 35497 : h = RgXn_mul2(fo,go,n0);
1496 35497 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1497 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1498 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1499 35497 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1500 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1501 35497 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1502 35497 : m = RgX_inflate(m,2);
1503 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1504 35497 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1505 35497 : h = RgX_inflate(h,2);
1506 35497 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1507 35497 : return gc_upto(av, h);
1508 : }
1509 : /* (f*g) \/ x^n */
1510 : static GEN
1511 1589510 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1512 : {
1513 1589510 : long d = degpol(f)+degpol(g) + 1 - n;
1514 : GEN h;
1515 1589510 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1516 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1517 : RgX_recip_i(g), d));
1518 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1519 : }
1520 :
1521 : /* (f*g) \/ x^n */
1522 : static GEN
1523 0 : RgX_sqrhigh_i2(GEN f, long n)
1524 : {
1525 0 : long d = 2*degpol(f)+ 1 - n;
1526 : GEN h;
1527 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1528 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1529 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1530 : }
1531 :
1532 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1533 : * b+2 were sent instead. na, nb = number of terms of a, b.
1534 : * Only c, c0, c1, c2 are genuine GEN.
1535 : */
1536 : GEN
1537 40273308 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1538 : {
1539 : GEN a0, c, c0;
1540 40273308 : long n0, n0a, i, v = 0;
1541 : pari_sp av;
1542 :
1543 65271616 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1544 61964024 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1545 40273288 : if (na < nb) swapspec(a,b, na,nb);
1546 40273288 : if (!nb) return pol_0(0);
1547 :
1548 34686284 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1549 479315 : RgX_shift_inplace_init(v);
1550 479328 : i = (na>>1); n0 = na-i; na = i;
1551 479328 : av = avma; a0 = a+n0; n0a = n0;
1552 1334187 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1553 :
1554 479328 : if (nb > n0)
1555 : {
1556 : GEN b0,c1,c2;
1557 : long n0b;
1558 :
1559 477945 : nb -= n0; b0 = b+n0; n0b = n0;
1560 1445075 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1561 477945 : c = RgX_mulspec(a,b,n0a,n0b);
1562 477945 : c0 = RgX_mulspec(a0,b0, na,nb);
1563 :
1564 477945 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1565 477945 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1566 :
1567 477945 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1568 477945 : c2 = RgX_sub(c1, RgX_add(c0,c));
1569 477945 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1570 : }
1571 : else
1572 : {
1573 1383 : c = RgX_mulspec(a,b,n0a,nb);
1574 1383 : c0 = RgX_mulspec(a0,b,na,nb);
1575 : }
1576 479328 : c0 = RgX_addmulXn(c0,c,n0);
1577 479328 : return RgX_shift_inplace(gc_upto(av,c0), v);
1578 : }
1579 :
1580 : INLINE GEN
1581 100082 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1582 : {
1583 100082 : pari_sp av = avma;
1584 100082 : GEN s = NULL;
1585 100082 : long j, l = (i+1)>>1;
1586 203924 : for (j=a; j<l; j++)
1587 : {
1588 103842 : GEN xj = gel(x,j), xx = gel(x,i-j);
1589 103842 : if (xj && xx)
1590 : {
1591 94707 : GEN t = gmul(xj, xx);
1592 94707 : s = s? gadd(s, t): t;
1593 : }
1594 : }
1595 100082 : if (s) s = gshift(s,1);
1596 100082 : if ((i&1) == 0)
1597 : {
1598 69041 : GEN t = gel(x, i>>1);
1599 69041 : if (t) {
1600 65982 : t = gsqr(t);
1601 65982 : s = s? gadd(s, t): t;
1602 : }
1603 : }
1604 100082 : return s? gc_upto(av,s): gen_0;
1605 : }
1606 : static GEN
1607 38000 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1608 : {
1609 : long i, lz, nz;
1610 : GEN z;
1611 :
1612 38000 : if (!nx) return pol_0(0);
1613 38000 : x = RgXspec_kill0(x,nx);
1614 38000 : lz = (nx << 1) + 1, nz = lz-2;
1615 38000 : lz += v;
1616 38000 : z = cgetg(lz,t_POL) + 2;
1617 94602 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1618 107041 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1619 69041 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1620 38000 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1621 : }
1622 : /* return x^2 mod t^n */
1623 : static GEN
1624 0 : RgXn_sqr_basecase(GEN x, long n)
1625 : {
1626 0 : long i, lz = n+2, lx = lgpol(x);
1627 : GEN z;
1628 0 : if (lx < 0) return pol_0(varn(x));
1629 0 : z = cgetg(lz, t_POL);
1630 0 : z[1] = x[1];
1631 0 : x+=2; if (lx > n) lx = n;
1632 0 : x = RgXspec_kill0(x,lx);
1633 0 : z+=2;/* x:z [i] = term of degree i */
1634 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1635 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1636 0 : z -= 2; return normalizepol_lg(z, lz);
1637 : }
1638 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1639 : static GEN
1640 315 : RgXn_sqr2(GEN f, long n)
1641 : {
1642 315 : pari_sp av = avma;
1643 : GEN fe,fo, l,h,m;
1644 : long n0, n1;
1645 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1646 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1647 0 : n0 = n>>1; n1 = n-n0;
1648 0 : RgX_even_odd(f, &fe, &fo);
1649 0 : l = RgXn_sqr(fe,n1);
1650 0 : h = RgXn_sqr(fo,n0);
1651 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1652 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1653 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1654 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1655 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1656 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1657 0 : m = RgX_inflate(m,2);
1658 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1659 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1660 0 : h = RgX_inflate(h,2);
1661 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1662 0 : return gc_upto(av, h);
1663 : }
1664 : GEN
1665 38039 : RgX_sqrspec(GEN a, long na)
1666 : {
1667 : GEN a0, c, c0, c1;
1668 38039 : long n0, n0a, i, v = 0;
1669 : pari_sp av;
1670 :
1671 66340 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1672 38039 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1673 39 : RgX_shift_inplace_init(v);
1674 39 : i = (na>>1); n0 = na-i; na = i;
1675 39 : av = avma; a0 = a+n0; n0a = n0;
1676 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1677 :
1678 39 : c = RgX_sqrspec(a,n0a);
1679 39 : c0 = RgX_sqrspec(a0,na);
1680 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1681 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1682 39 : c0 = RgX_addmulXn(c0,c,n0);
1683 39 : return RgX_shift_inplace(gc_upto(av,c0), v);
1684 : }
1685 :
1686 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1687 : GEN
1688 1719192 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1689 : {
1690 1719192 : GEN z = RgX_mul(A, B);
1691 1719180 : if (a < b)
1692 10421 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1693 1708759 : else if (a > b)
1694 1104946 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1695 : else
1696 603813 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1697 1719190 : return z;
1698 : }
1699 :
1700 : GEN
1701 38836724 : RgX_mul_i(GEN x, GEN y)
1702 : {
1703 38836724 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1704 38836359 : setvarn(z, varn(x)); return z;
1705 : }
1706 :
1707 : GEN
1708 37961 : RgX_sqr_i(GEN x)
1709 : {
1710 37961 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1711 37961 : setvarn(z,varn(x)); return z;
1712 : }
1713 :
1714 : /*******************************************************************/
1715 : /* */
1716 : /* DIVISION */
1717 : /* */
1718 : /*******************************************************************/
1719 : GEN
1720 5262682 : RgX_Rg_divexact(GEN x, GEN y) {
1721 5262682 : long i, lx = lg(x);
1722 : GEN z;
1723 5262682 : if (lx == 2) return gcopy(x);
1724 5212437 : switch(typ(y))
1725 : {
1726 5101476 : case t_INT:
1727 5101476 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1728 5036063 : break;
1729 5243 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1730 : }
1731 5141781 : z = cgetg(lx, t_POL); z[1] = x[1];
1732 29095726 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1733 5140642 : return z;
1734 : }
1735 : GEN
1736 30727629 : RgX_Rg_div(GEN x, GEN y) {
1737 30727629 : long i, lx = lg(x);
1738 : GEN z;
1739 30727629 : if (lx == 2) return gcopy(x);
1740 30453415 : switch(typ(y))
1741 : {
1742 21220289 : case t_INT:
1743 21220289 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1744 4275893 : break;
1745 6069 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1746 : }
1747 13502950 : z = cgetg(lx, t_POL); z[1] = x[1];
1748 50412464 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1749 13503005 : return normalizepol_lg(z, lx);
1750 : }
1751 : GEN
1752 39711 : RgX_normalize(GEN x)
1753 : {
1754 39711 : GEN z, d = NULL;
1755 39711 : long i, n = lg(x)-1;
1756 39711 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1757 39711 : if (i == 1) return pol_0(varn(x));
1758 39711 : if (i == n && isint1(d)) return x;
1759 18200 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1760 32571 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1761 18200 : gel(z,n) = Rg_get_1(d); return z;
1762 : }
1763 : GEN
1764 10948 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1765 : GEN
1766 265485 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1767 : {
1768 265485 : long l = lg(a), i;
1769 : GEN a0, z0, z;
1770 :
1771 265485 : if (l <= 3)
1772 : {
1773 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1774 0 : return pol_0(varn(a));
1775 : }
1776 265485 : z = cgetg(l-1, t_POL);
1777 265485 : z[1] = a[1];
1778 265485 : a0 = a + l-1;
1779 265485 : z0 = z + l-2; *z0 = *a0--;
1780 3024816 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1781 : {
1782 2759333 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1783 2759331 : gel(z0,0) = t;
1784 : }
1785 265483 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1786 265483 : return z;
1787 : }
1788 : /* Polynomial division x / y:
1789 : * if pr = ONLY_REM return remainder, otherwise return quotient
1790 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1791 : * if pr != NULL set *pr to remainder, as the last object on stack */
1792 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1793 : static GEN
1794 23432217 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1795 : {
1796 : pari_sp avy, av, av1;
1797 : long dx,dy,dz,i,j,sx,lr;
1798 : GEN z,p1,p2,rem,y_lead,mod,p;
1799 : GEN (*f)(GEN,GEN);
1800 :
1801 23432217 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1802 :
1803 23432217 : dy = degpol(y);
1804 23432199 : y_lead = gel(y,dy+2);
1805 23432199 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1806 : {
1807 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1808 0 : for (dy--; dy>=0; dy--)
1809 : {
1810 0 : y_lead = gel(y,dy+2);
1811 0 : if (!gequal0(y_lead)) break;
1812 : }
1813 : }
1814 23432184 : if (!dy) /* y is constant */
1815 : {
1816 6928 : if (pr == ONLY_REM) return pol_0(varn(x));
1817 6928 : z = RgX_Rg_div(x, y_lead);
1818 6928 : if (pr == ONLY_DIVIDES) return z;
1819 2714 : if (pr) *pr = pol_0(varn(x));
1820 2714 : return z;
1821 : }
1822 23425256 : dx = degpol(x);
1823 23425220 : if (dx < dy)
1824 : {
1825 3848456 : if (pr == ONLY_REM) return RgX_copy(x);
1826 348591 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1827 348570 : z = pol_0(varn(x));
1828 348570 : if (pr) *pr = RgX_copy(x);
1829 348570 : return z;
1830 : }
1831 :
1832 : /* x,y in R[X], y non constant */
1833 19576764 : av = avma;
1834 19576764 : p = NULL;
1835 19576764 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1836 : {
1837 130277 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1838 130277 : if (!z) return gc_NULL(av);
1839 130277 : z = FpX_to_mod(z, p);
1840 130277 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1841 71134 : return gc_upto(av, z);
1842 59143 : *pr = FpX_to_mod(*pr, p);
1843 59143 : return gc_all(av, 2, &z, pr);
1844 : }
1845 19446602 : switch(typ(y_lead))
1846 : {
1847 0 : case t_REAL:
1848 0 : y_lead = ginv(y_lead);
1849 0 : f = gmul; mod = NULL;
1850 0 : break;
1851 1878 : case t_INTMOD:
1852 1878 : case t_POLMOD: y_lead = ginv(y_lead);
1853 1878 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1854 1878 : break;
1855 19444724 : default: if (gequal1(y_lead)) y_lead = NULL;
1856 19444715 : f = gdiv; mod = NULL;
1857 : }
1858 :
1859 19446593 : if (y_lead == NULL)
1860 17486401 : p2 = gel(x,dx+2);
1861 : else {
1862 : for(;;) {
1863 1960192 : p2 = f(gel(x,dx+2),y_lead);
1864 1960192 : p2 = simplify_shallow(p2);
1865 1960192 : if (!isexactzero(p2) || (--dx < 0)) break;
1866 : }
1867 1960192 : if (dx < dy) /* leading coeff of x was in fact zero */
1868 : {
1869 0 : if (pr == ONLY_DIVIDES) {
1870 0 : set_avma(av);
1871 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1872 : }
1873 0 : if (pr == ONLY_REM)
1874 : {
1875 0 : if (dx < 0)
1876 0 : return gc_GEN(av, scalarpol(p2, varn(x)));
1877 : else
1878 : {
1879 : GEN t;
1880 0 : set_avma(av);
1881 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1882 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1883 0 : return t;
1884 : }
1885 : }
1886 0 : if (pr) /* cf ONLY_REM above */
1887 : {
1888 0 : if (dx < 0)
1889 : {
1890 0 : p2 = gclone(p2);
1891 0 : set_avma(av);
1892 0 : z = pol_0(varn(x));
1893 0 : x = scalarpol(p2, varn(x));
1894 0 : gunclone(p2);
1895 : }
1896 : else
1897 : {
1898 : GEN t;
1899 0 : set_avma(av);
1900 0 : z = pol_0(varn(x));
1901 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1902 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1903 0 : x = t;
1904 : }
1905 0 : *pr = x;
1906 : }
1907 : else
1908 : {
1909 0 : set_avma(av);
1910 0 : z = pol_0(varn(x));
1911 : }
1912 0 : return z;
1913 : }
1914 : }
1915 : /* dx >= dy */
1916 19446593 : avy = avma;
1917 19446593 : dz = dx-dy;
1918 19446593 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1919 19446573 : x += 2;
1920 19446573 : z += 2;
1921 19446573 : y += 2;
1922 19446573 : gel(z,dz) = gcopy(p2);
1923 :
1924 53704988 : for (i=dx-1; i>=dy; i--)
1925 : {
1926 34258533 : av1=avma; p1=gel(x,i);
1927 1139447644 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1928 34253737 : if (y_lead) p1 = simplify(f(p1,y_lead));
1929 :
1930 34253737 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1931 : else
1932 24266149 : p1 = avma==av1? gcopy(p1): gc_upto(av1,p1);
1933 34257970 : gel(z,i-dy) = p1;
1934 : }
1935 19446455 : if (!pr) return gc_upto(av,z-2);
1936 :
1937 12522125 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1938 12877928 : for (sx=0; ; i--)
1939 : {
1940 12877928 : p1 = gel(x,i);
1941 : /* we always enter this loop at least once */
1942 32134851 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1943 12876784 : if (mod && avma==av1) p1 = gmul(p1,mod);
1944 12876784 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1945 3575760 : if (!isexactzero(p1)) break;
1946 3425602 : if (!i) break;
1947 355742 : set_avma(av1);
1948 : }
1949 12521415 : if (pr == ONLY_DIVIDES)
1950 : {
1951 8078 : if (sx) return gc_NULL(av);
1952 8057 : set_avma((pari_sp)rem); return gc_upto(av,z-2);
1953 : }
1954 12513337 : lr=i+3; rem -= lr;
1955 12513337 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1956 12471048 : else p1 = gc_upto((pari_sp)rem,p1);
1957 12513961 : rem[0] = evaltyp(t_POL) | _evallg(lr);
1958 12513961 : rem[1] = z[-1];
1959 12513961 : rem += 2;
1960 12513961 : gel(rem,i) = p1;
1961 17219670 : for (i--; i>=0; i--)
1962 : {
1963 4705721 : av1=avma; p1 = gel(x,i);
1964 13564027 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1965 4705226 : if (mod && avma==av1) p1 = gmul(p1,mod);
1966 4705444 : gel(rem,i) = avma==av1? gcopy(p1):gc_upto(av1,p1);
1967 : }
1968 12513949 : rem -= 2;
1969 12513949 : if (!sx) (void)normalizepol_lg(rem, lr);
1970 12514181 : if (pr == ONLY_REM) return gc_upto(av,rem);
1971 6762732 : z -= 2; *pr = rem; return gc_all_unsafe(av,avy,2,&z,pr);
1972 : }
1973 :
1974 : GEN
1975 14180963 : RgX_divrem(GEN x, GEN y, GEN *pr)
1976 : {
1977 14180963 : if (pr == ONLY_REM) return RgX_rem(x, y);
1978 14180963 : return RgX_divrem_i(x, y, pr);
1979 : }
1980 :
1981 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1982 : GEN
1983 156857 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1984 : {
1985 156857 : long vx = varn(x), dx = degpol(x), dy = degpol(y), dz, i, j, sx, lr;
1986 : pari_sp av0, av;
1987 : GEN z, p1, rem, lead;
1988 :
1989 156857 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1990 156857 : if (dx < dy)
1991 : {
1992 41644 : if (pr)
1993 : {
1994 41616 : av0 = avma; x = RgXQX_red(x, T);
1995 41616 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1996 41609 : if (pr == ONLY_REM) return x;
1997 0 : *pr = x;
1998 : }
1999 28 : return pol_0(vx);
2000 : }
2001 115213 : lead = leading_coeff(y);
2002 115213 : if (!dy) /* y is constant */
2003 : {
2004 602 : if (pr && pr != ONLY_DIVIDES)
2005 : {
2006 0 : if (pr == ONLY_REM) return pol_0(vx);
2007 0 : *pr = pol_0(vx);
2008 : }
2009 602 : if (gequal1(lead)) return RgX_copy(x);
2010 0 : av0 = avma; x = gmul(x, ginvmod(lead,T));
2011 0 : return gc_upto(av0, RgXQX_red(x,T));
2012 : }
2013 114611 : av0 = avma; dz = dx-dy;
2014 114611 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
2015 114611 : set_avma(av0);
2016 114611 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2017 114611 : x += 2; y += 2; z += 2;
2018 :
2019 114611 : p1 = gel(x,dx); av = avma;
2020 114611 : gel(z,dz) = lead? gc_upto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2021 581827 : for (i=dx-1; i>=dy; i--)
2022 : {
2023 467216 : av = avma; p1 = gel(x,i);
2024 2418679 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2025 467185 : if (lead) p1 = gmul(grem(p1, T), lead);
2026 467191 : gel(z,i-dy) = gc_upto(av, grem(p1, T));
2027 : }
2028 114611 : if (!pr) { guncloneNULL(lead); return z-2; }
2029 :
2030 113337 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2031 191645 : for (sx=0; ; i--)
2032 : {
2033 191645 : p1 = gel(x,i);
2034 663081 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2035 191641 : p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2036 105909 : if (!i) break;
2037 78308 : set_avma(av);
2038 : }
2039 113337 : if (pr == ONLY_DIVIDES)
2040 : {
2041 27143 : guncloneNULL(lead);
2042 27143 : if (sx) return gc_NULL(av0);
2043 24577 : return gc_const((pari_sp)rem, z-2);
2044 : }
2045 86194 : lr=i+3; rem -= lr; av = (pari_sp)rem;
2046 86194 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2047 86194 : rem[1] = z[-1];
2048 86194 : rem += 2; gel(rem,i) = gc_upto(av, p1);
2049 187610 : for (i--; i>=0; i--)
2050 : {
2051 101416 : av = avma; p1 = gel(x,i);
2052 337537 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2053 101416 : gel(rem,i) = gc_upto(av, grem(p1, T));
2054 : }
2055 86194 : rem -= 2;
2056 86194 : guncloneNULL(lead);
2057 86194 : if (!sx) (void)normalizepol_lg(rem, lr);
2058 86194 : if (pr == ONLY_REM) return gc_upto(av0,rem);
2059 168 : *pr = rem; return z-2;
2060 : }
2061 :
2062 : /*******************************************************************/
2063 : /* */
2064 : /* PSEUDO-DIVISION */
2065 : /* */
2066 : /*******************************************************************/
2067 : INLINE GEN
2068 1050984 : rem(GEN c, GEN T)
2069 : {
2070 1050984 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2071 1050984 : return c;
2072 : }
2073 :
2074 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2075 : int
2076 17734 : ZXQX_dvd(GEN x, GEN y, GEN T)
2077 : {
2078 : long dx, dy, i, T_ismonic;
2079 17734 : pari_sp av = avma, av2;
2080 : GEN y_lead;
2081 :
2082 17734 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2083 17734 : dy = degpol(y); y_lead = gel(y,dy+2);
2084 17734 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2085 : /* if monic, no point in using pseudo-division */
2086 17734 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2087 14752 : T_ismonic = gequal1(leading_coeff(T));
2088 14752 : dx = degpol(x);
2089 14752 : if (dx < dy) return !signe(x);
2090 14752 : (void)new_chunk(2);
2091 14752 : x = RgX_recip_i(x)+2;
2092 14752 : y = RgX_recip_i(y)+2;
2093 : /* pay attention to sparse divisors */
2094 30005 : for (i = 1; i <= dy; i++)
2095 15253 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2096 14752 : av2 = avma;
2097 : for (;;)
2098 72967 : {
2099 87719 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2100 87719 : x0 = gneg(x0);
2101 87719 : m = gcdii(cx, y0);
2102 87719 : if (!equali1(m))
2103 : {
2104 85217 : x0 = gdiv(x0, m);
2105 85217 : y0 = diviiexact(y0, m);
2106 85217 : if (equali1(y0)) y0 = NULL;
2107 : }
2108 179211 : for (i=1; i<=dy; i++)
2109 : {
2110 91492 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2111 91492 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2112 91492 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2113 91492 : gel(x,i) = c;
2114 : }
2115 688014 : for ( ; i<=dx; i++)
2116 : {
2117 600295 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2118 600295 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2119 600295 : gel(x,i) = c;
2120 : }
2121 102209 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2122 87719 : if (dx < dy) break;
2123 72967 : if (gc_needed(av2,1))
2124 : {
2125 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2126 28 : gc_slice(av2,x,dx+1);
2127 : }
2128 : }
2129 14752 : return gc_bool(av, dx < 0);
2130 : }
2131 :
2132 : /* T either NULL or a t_POL. */
2133 : GEN
2134 138872 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2135 : {
2136 138872 : long vx = varn(x), dx, dy, dz, i, lx, p;
2137 138872 : pari_sp av = avma, av2;
2138 : GEN y_lead;
2139 :
2140 138872 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2141 138872 : dy = degpol(y); y_lead = gel(y,dy+2);
2142 : /* if monic, no point in using pseudo-division */
2143 138872 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2144 125095 : dx = degpol(x);
2145 125095 : if (dx < dy) return RgX_copy(x);
2146 125088 : (void)new_chunk(2);
2147 125088 : x = RgX_recip_i(x)+2;
2148 125088 : y = RgX_recip_i(y)+2;
2149 : /* pay attention to sparse divisors */
2150 535566 : for (i = 1; i <= dy; i++)
2151 410478 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2152 125088 : dz = dx-dy; p = dz+1;
2153 125088 : av2 = avma;
2154 : for (;;)
2155 : {
2156 247161 : gel(x,0) = gneg(gel(x,0)); p--;
2157 1061861 : for (i=1; i<=dy; i++)
2158 : {
2159 814703 : GEN c = gmul(y_lead, gel(x,i));
2160 814690 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2161 814700 : gel(x,i) = rem(c, T);
2162 : }
2163 373910 : for ( ; i<=dx; i++)
2164 : {
2165 126749 : GEN c = gmul(y_lead, gel(x,i));
2166 126749 : gel(x,i) = rem(c, T);
2167 : }
2168 253361 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2169 247161 : if (dx < dy) break;
2170 122073 : if (gc_needed(av2,1))
2171 : {
2172 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2173 0 : gc_slice(av2,x,dx+1);
2174 : }
2175 : }
2176 125088 : if (dx < 0) return pol_0(vx);
2177 124927 : lx = dx+3; x -= 2;
2178 124927 : x[0] = evaltyp(t_POL) | _evallg(lx);
2179 124927 : x[1] = evalsigne(1) | evalvarn(vx);
2180 124927 : x = RgX_recip_i(x);
2181 124927 : if (p)
2182 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2183 3884 : GEN t = y_lead;
2184 3884 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2185 0 : t = RgXQ_powu(t, p, T);
2186 : else
2187 3884 : t = gpowgs(t, p);
2188 11803 : for (i=2; i<lx; i++)
2189 : {
2190 7919 : GEN c = gmul(gel(x,i), t);
2191 7919 : gel(x,i) = rem(c,T);
2192 : }
2193 3884 : if (!T) return gc_upto(av, x);
2194 : }
2195 121043 : return gc_GEN(av, x);
2196 : }
2197 :
2198 : GEN
2199 138872 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2200 :
2201 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2202 : GEN
2203 12063 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2204 : {
2205 12063 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2206 12063 : pari_sp av = avma, av2;
2207 : GEN z, r, ypow, y_lead;
2208 :
2209 12063 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2210 12063 : dy = degpol(y); y_lead = gel(y,dy+2);
2211 12063 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2212 7561 : dx = degpol(x);
2213 7561 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2214 7561 : if (dx == dy)
2215 : {
2216 98 : GEN x_lead = gel(x,lg(x)-1);
2217 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2218 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2219 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2220 98 : *ptr = gc_upto(av, r); return scalarpol(x_lead, vx);
2221 : }
2222 7463 : (void)new_chunk(2);
2223 7463 : x = RgX_recip_i(x)+2;
2224 7463 : y = RgX_recip_i(y)+2;
2225 : /* pay attention to sparse divisors */
2226 39000 : for (i = 1; i <= dy; i++)
2227 31537 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2228 7463 : dz = dx-dy; p = dz+1;
2229 7463 : lz = dz+3;
2230 7463 : z = cgetg(lz, t_POL);
2231 7463 : z[1] = evalsigne(1) | evalvarn(vx);
2232 28640 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2233 7463 : ypow = new_chunk(dz+1);
2234 7463 : gel(ypow,0) = gen_1;
2235 7463 : gel(ypow,1) = y_lead;
2236 13714 : for (i=2; i<=dz; i++)
2237 : {
2238 6251 : GEN c = gmul(gel(ypow,i-1), y_lead);
2239 6251 : gel(ypow,i) = rem(c,T);
2240 : }
2241 7463 : av2 = avma;
2242 7463 : for (iz=2;;)
2243 : {
2244 15605 : p--;
2245 15605 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2246 69915 : for (i=1; i<=dy; i++)
2247 : {
2248 54310 : GEN c = gmul(y_lead, gel(x,i));
2249 54310 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2250 54310 : gel(x,i) = rem(c, T);
2251 : }
2252 41058 : for ( ; i<=dx; i++)
2253 : {
2254 25453 : GEN c = gmul(y_lead, gel(x,i));
2255 25453 : gel(x,i) = rem(c,T);
2256 : }
2257 15605 : x++; dx--;
2258 21177 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2259 15605 : if (dx < dy) break;
2260 8142 : if (gc_needed(av2,1))
2261 : {
2262 0 : GEN X = x-2;
2263 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2264 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2265 0 : (void)gc_all(av2,2, &X, &z); x = X+2;
2266 : }
2267 : }
2268 14008 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2269 7463 : if (dx < 0)
2270 182 : x = pol_0(vx);
2271 : else
2272 : {
2273 7281 : lx = dx+3; x -= 2;
2274 7281 : x[0] = evaltyp(t_POL) | _evallg(lx);
2275 7281 : x[1] = evalsigne(1) | evalvarn(vx);
2276 7281 : x = RgX_recip_i(x);
2277 : }
2278 7463 : z = RgX_recip_i(z);
2279 7463 : r = x;
2280 7463 : if (p)
2281 : {
2282 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2283 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2284 : }
2285 7463 : *ptr = r; return gc_all(av, 2, &z, ptr);
2286 : }
2287 : GEN
2288 11818 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2289 11818 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2290 :
2291 : GEN
2292 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2293 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2294 : GEN
2295 743683538 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2296 : GEN
2297 141351 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2298 : GEN
2299 26320 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2300 : GEN
2301 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2302 : GEN
2303 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2304 :
2305 : GEN
2306 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2307 :
2308 : GEN
2309 0 : RgXQX_powers(GEN P, long n, GEN T)
2310 : {
2311 0 : GEN v = cgetg(n+2, t_VEC);
2312 : long i;
2313 0 : gel(v, 1) = pol_1(varn(T));
2314 0 : if (n==0) return v;
2315 0 : gel(v, 2) = gcopy(P);
2316 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2317 0 : return v;
2318 : }
2319 :
2320 : static GEN
2321 623391 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2322 : static GEN
2323 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2324 : static GEN
2325 74882 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2326 : static GEN
2327 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2328 : static GEN
2329 280996 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2330 : static GEN
2331 1040655 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2332 : static GEN
2333 1019733 : _one(void *data) { return pol_1(varn((GEN)data)); }
2334 : static GEN
2335 1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2336 : static GEN
2337 724293 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2338 :
2339 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2340 : _mul, _sqr, _one, _zero };
2341 :
2342 : GEN
2343 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2344 : {
2345 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2346 : }
2347 :
2348 : GEN
2349 417523 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2350 : {
2351 417523 : int use_sqr = 2*degpol(x) >= degpol(T);
2352 417521 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2353 : }
2354 :
2355 : /* mod X^n */
2356 : struct modXn {
2357 : long v; /* varn(X) */
2358 : long n;
2359 : } ;
2360 : static GEN
2361 11893 : _sqrXn(void *data, GEN x) {
2362 11893 : struct modXn *S = (struct modXn*)data;
2363 11893 : return RgXn_sqr(x, S->n);
2364 : }
2365 : static GEN
2366 4528 : _mulXn(void *data, GEN x, GEN y) {
2367 4528 : struct modXn *S = (struct modXn*)data;
2368 4528 : return RgXn_mul(x,y, S->n);
2369 : }
2370 : static GEN
2371 1939 : _oneXn(void *data) {
2372 1939 : struct modXn *S = (struct modXn*)data;
2373 1939 : return pol_1(S->v);
2374 : }
2375 : static GEN
2376 0 : _zeroXn(void *data) {
2377 0 : struct modXn *S = (struct modXn*)data;
2378 0 : return pol_0(S->v);
2379 : }
2380 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2381 : _oneXn, _zeroXn };
2382 :
2383 : GEN
2384 483 : RgXn_powers(GEN x, long m, long n)
2385 : {
2386 483 : long d = degpol(x);
2387 483 : int use_sqr = (d<<1) >= n;
2388 : struct modXn S;
2389 483 : S.v = varn(x); S.n = n;
2390 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2391 : }
2392 :
2393 : GEN
2394 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2395 : {
2396 : struct modXn S;
2397 : long v;
2398 2286 : if (n == 0) return x;
2399 2286 : v = RgX_valrem(x, &x);
2400 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2401 2265 : S.v = varn(x); S.n = n;
2402 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2403 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2404 2265 : return x;
2405 : }
2406 : GEN
2407 0 : RgXn_powu(GEN x, ulong m, long n)
2408 : {
2409 : pari_sp av;
2410 0 : if (n == 0) return gcopy(x);
2411 0 : av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
2412 : }
2413 :
2414 : GEN
2415 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2416 : {
2417 : struct modXn S;
2418 966 : S.v = varn(gel(x,2)); S.n = n;
2419 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2420 : }
2421 :
2422 : GEN
2423 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2424 : {
2425 0 : int use_sqr = 2*degpol(x) >= n;
2426 : struct modXn S;
2427 0 : S.v = varn(x); S.n = n;
2428 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2429 : }
2430 :
2431 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2432 : GEN
2433 5187 : RgXn_eval(GEN Q, GEN x, long n)
2434 : {
2435 5187 : long d = degpol(x);
2436 : int use_sqr;
2437 : struct modXn S;
2438 5187 : if (d == 1 && isrationalzero(gel(x,2)))
2439 : {
2440 5180 : GEN y = RgX_unscale(Q, gel(x,3));
2441 5180 : setvarn(y, varn(x)); return y;
2442 : }
2443 7 : S.v = varn(x);
2444 7 : S.n = n;
2445 7 : use_sqr = (d<<1) >= n;
2446 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2447 : }
2448 :
2449 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2450 : static GEN
2451 2614799 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2452 : {
2453 2614799 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2454 2614799 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2455 : }
2456 :
2457 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2458 : static GEN
2459 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2460 : {
2461 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2462 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2463 : }
2464 :
2465 : static GEN
2466 2188985 : RgXn_div_gen(GEN g, GEN f, long e)
2467 : {
2468 : pari_sp av;
2469 : ulong mask;
2470 : GEN W, a;
2471 2188985 : long v = varn(f), n = 1;
2472 :
2473 2188985 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2474 2188978 : a = ginv(gel(f,2));
2475 2188977 : if (e == 1 && !g) return scalarpol(a, v);
2476 2184063 : else if (e == 2 && !g)
2477 : {
2478 : GEN b;
2479 836932 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2480 246250 : b = gneg(b);
2481 246248 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2482 246248 : return deg1pol(b, a, v);
2483 : }
2484 1347131 : av = avma;
2485 1347131 : W = scalarpol_shallow(a,v);
2486 1347131 : mask = quadratic_prec_mask(e);
2487 3955973 : while (mask > 1)
2488 : {
2489 : GEN u, fr;
2490 2608842 : long n2 = n;
2491 2608842 : n<<=1; if (mask & 1) n--;
2492 2608842 : mask >>= 1;
2493 2608842 : fr = RgXn_red_shallow(f, n);
2494 2608842 : if (mask>1 || !g)
2495 : {
2496 1324855 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2497 1324855 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2498 : }
2499 : else
2500 : {
2501 1283987 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2502 1283987 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2503 1283987 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2504 : }
2505 2608842 : if (gc_needed(av,2))
2506 : {
2507 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2508 0 : W = gc_upto(av, W);
2509 : }
2510 : }
2511 1347131 : return W;
2512 : }
2513 :
2514 : static GEN
2515 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2516 : {
2517 : GEN r;
2518 119 : if (lgefint(p) == 3)
2519 : {
2520 119 : ulong pp = uel(p, 2);
2521 119 : if (pp == 2)
2522 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2523 : else
2524 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2525 : }
2526 : else
2527 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2528 119 : return FpX_to_mod(r, p);
2529 : }
2530 :
2531 : static GEN
2532 7 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2533 : {
2534 7 : GEN r, T = RgX_to_FpX(pol, p);
2535 7 : if (signe(T) == 0) pari_err_OP("/", x, y);
2536 7 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2537 7 : return FpXQX_to_mod(r, T, p);
2538 : }
2539 :
2540 : static GEN
2541 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2542 : {
2543 : GEN r;
2544 91 : if (lgefint(p) == 3)
2545 : {
2546 91 : ulong pp = uel(p, 2);
2547 91 : if (pp == 2)
2548 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2549 : else
2550 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2551 : }
2552 : else
2553 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2554 91 : return FpX_to_mod(r, p);
2555 : }
2556 :
2557 : static GEN
2558 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2559 : {
2560 0 : GEN r, T = RgX_to_FpX(pol, p);
2561 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2562 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2563 0 : return FpXQX_to_mod(r, T, p);
2564 : }
2565 :
2566 : static GEN
2567 905087 : RgXn_inv_fast(GEN x, long e)
2568 : {
2569 : GEN p, pol;
2570 : long pa;
2571 905087 : long t = RgX_type(x,&p,&pol,&pa);
2572 905089 : switch(t)
2573 : {
2574 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2575 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2576 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2577 904998 : default: return NULL;
2578 : }
2579 : }
2580 :
2581 : static GEN
2582 1284113 : RgXn_div_fast(GEN x, GEN y, long e)
2583 : {
2584 : GEN p, pol;
2585 : long pa;
2586 1284113 : long t = RgX_type2(x,y,&p,&pol,&pa);
2587 1284113 : switch(t)
2588 : {
2589 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2590 7 : case RgX_type_code(t_POLMOD, t_INTMOD):
2591 7 : return RgXn_div_FpXQX(x, y, e, pol, p);
2592 1283987 : default: return NULL;
2593 : }
2594 : }
2595 :
2596 : GEN
2597 1284113 : RgXn_div_i(GEN g, GEN f, long e)
2598 : {
2599 1284113 : GEN h = RgXn_div_fast(g, f, e);
2600 1284113 : if (h) return h;
2601 1283987 : return RgXn_div_gen(g, f, e);
2602 : }
2603 :
2604 : GEN
2605 584450 : RgXn_div(GEN g, GEN f, long e)
2606 : {
2607 584450 : pari_sp av = avma;
2608 584450 : return gc_upto(av, RgXn_div_i(g, f, e));
2609 : }
2610 :
2611 : GEN
2612 905087 : RgXn_inv_i(GEN f, long e)
2613 : {
2614 905087 : GEN h = RgXn_inv_fast(f, e);
2615 905089 : if (h) return h;
2616 904998 : return RgXn_div_gen(NULL, f, e);
2617 : }
2618 :
2619 : GEN
2620 807479 : RgXn_inv(GEN f, long e)
2621 : {
2622 807479 : pari_sp av = avma;
2623 807479 : return gc_upto(av, RgXn_inv_i(f, e));
2624 : }
2625 :
2626 : /* intformal(x^n*S) / x^(n+1) */
2627 : static GEN
2628 56377 : RgX_integXn(GEN x, long n)
2629 114355 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2630 :
2631 : GEN
2632 52918 : RgXn_expint(GEN h, long e)
2633 : {
2634 52918 : pari_sp av = avma, av2;
2635 52918 : long v = varn(h), n;
2636 52918 : GEN f = pol_1(v), g;
2637 : ulong mask;
2638 :
2639 52918 : if (!signe(h)) return f;
2640 50426 : g = pol_1(v);
2641 50426 : n = 1; mask = quadratic_prec_mask(e);
2642 50426 : av2 = avma;
2643 56376 : for (;mask>1;)
2644 : {
2645 : GEN u, w;
2646 56376 : long n2 = n;
2647 56376 : n<<=1; if (mask & 1) n--;
2648 56376 : mask >>= 1;
2649 56376 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2650 56377 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2651 56377 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2652 56375 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2653 56375 : if (mask<=1) break;
2654 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2655 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2656 5950 : if (gc_needed(av2,2))
2657 : {
2658 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2659 0 : (void)gc_all(av2, 2, &f, &g);
2660 : }
2661 : }
2662 50425 : return gc_upto(av, f);
2663 : }
2664 :
2665 : GEN
2666 0 : RgXn_exp(GEN h, long e)
2667 : {
2668 0 : long d = degpol(h);
2669 0 : if (d < 0) return pol_1(varn(h));
2670 0 : if (!d || !gequal0(gel(h,2)))
2671 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2672 0 : return RgXn_expint(RgX_deriv(h), e);
2673 : }
2674 :
2675 : GEN
2676 154 : RgXn_reverse(GEN f, long e)
2677 : {
2678 154 : pari_sp av = avma, av2;
2679 : ulong mask;
2680 : GEN fi, a, df, W, an;
2681 154 : long v = varn(f), n=1;
2682 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2683 0 : pari_err_INV("serreverse",f);
2684 154 : fi = ginv(gel(f,3));
2685 154 : a = deg1pol_shallow(fi,gen_0,v);
2686 154 : if (e <= 2) return gc_GEN(av, a);
2687 133 : W = scalarpol(fi,v);
2688 133 : df = RgX_deriv(f);
2689 133 : mask = quadratic_prec_mask(e);
2690 133 : av2 = avma;
2691 616 : for (;mask>1;)
2692 : {
2693 : GEN u, fa, fr;
2694 483 : long n2 = n, rt;
2695 483 : n<<=1; if (mask & 1) n--;
2696 483 : mask >>= 1;
2697 483 : fr = RgXn_red_shallow(f, n);
2698 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2699 483 : an = RgXn_powers(a, rt, n);
2700 483 : if (n>1)
2701 : {
2702 483 : long n4 = (n2+1)>>1;
2703 483 : GEN dfr = RgXn_red_shallow(df, n2);
2704 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2705 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2706 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2707 : }
2708 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2709 483 : fa = RgX_shift(fa, -n2);
2710 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2711 483 : if (gc_needed(av2,2))
2712 : {
2713 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2714 0 : (void)gc_all(av2, 2, &a, &W);
2715 : }
2716 : }
2717 133 : return gc_upto(av, a);
2718 : }
2719 :
2720 : GEN
2721 7 : RgXn_sqrt(GEN h, long e)
2722 : {
2723 7 : pari_sp av = avma, av2;
2724 7 : long v = varn(h), n = 1;
2725 7 : GEN f = scalarpol(gen_1, v), df = f;
2726 7 : ulong mask = quadratic_prec_mask(e);
2727 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2728 0 : pari_err_SQRTN("RgXn_sqrt",h);
2729 7 : av2 = avma;
2730 : while(1)
2731 7 : {
2732 14 : long n2 = n, m;
2733 : GEN g;
2734 14 : n<<=1; if (mask & 1) n--;
2735 14 : mask >>= 1;
2736 14 : m = n-n2;
2737 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2738 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2739 14 : if (mask==1) return gc_upto(av, f);
2740 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2741 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2742 7 : if (gc_needed(av2,2))
2743 : {
2744 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2745 0 : (void)gc_all(av2, 2, &f, &df);
2746 : }
2747 : }
2748 : }
2749 :
2750 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2751 : GEN
2752 120509 : RgXQ_powu(GEN x, ulong n, GEN T)
2753 : {
2754 120509 : pari_sp av = avma;
2755 :
2756 120509 : if (!n) return pol_1(varn(x));
2757 72249 : if (n == 1) return RgX_copy(x);
2758 23604 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2759 23604 : return gc_GEN(av, x);
2760 : }
2761 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2762 : GEN
2763 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2764 : {
2765 : pari_sp av;
2766 102160 : long s = signe(n);
2767 :
2768 102160 : if (!s) return pol_1(varn(x));
2769 102160 : if (is_pm1(n) == 1)
2770 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2771 13853 : av = avma;
2772 13853 : if (s < 0) x = RgXQ_inv(x, T);
2773 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2774 13853 : return gc_GEN(av, x);
2775 : }
2776 : static GEN
2777 200611 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2778 : static GEN
2779 111195 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2780 :
2781 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2782 : GEN
2783 12145 : ZXQ_powers(GEN x, long l, GEN T)
2784 : {
2785 12145 : int use_sqr = 2*degpol(x) >= degpol(T);
2786 12145 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2787 : }
2788 :
2789 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2790 : GEN
2791 168511 : ZXQ_powu(GEN x, ulong n, GEN T)
2792 : {
2793 168511 : pari_sp av = avma;
2794 :
2795 168511 : if (!n) return pol_1(varn(x));
2796 168511 : if (n == 1) return ZX_copy(x);
2797 113876 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2798 113877 : return gc_GEN(av, x);
2799 : }
2800 :
2801 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2802 : GEN
2803 8610 : RgXQ_powers(GEN x, long l, GEN T)
2804 : {
2805 8610 : int use_sqr = 2*degpol(x) >= degpol(T);
2806 8610 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2807 : }
2808 :
2809 : GEN
2810 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2811 : {
2812 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2813 : }
2814 :
2815 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2816 : GEN
2817 9590 : QXQ_powers(GEN a, long n, GEN T)
2818 : {
2819 : GEN den, v;
2820 9590 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2821 9562 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2822 : /* den*a integral; v[i+1] = (den*a)^i in K */
2823 9562 : if (den)
2824 : { /* restore denominators */
2825 6012 : GEN d = den;
2826 : long i;
2827 6012 : gel(v,2) = a;
2828 27749 : for (i=3; i<=n+1; i++) {
2829 21737 : d = mulii(d,den);
2830 21737 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2831 : }
2832 : }
2833 9562 : return v;
2834 : }
2835 :
2836 : static GEN
2837 3612 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2838 : {
2839 3612 : long l, i, m = 0;
2840 : GEN dz, z;
2841 3612 : GEN V = cgetg_copy(v, &l);
2842 12796 : for (i = imin; i < l; i++)
2843 : {
2844 9184 : GEN c = gel(v, i);
2845 9184 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2846 : }
2847 3612 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2848 3983 : for (i = 1; i < imin; i++) V[i] = v[i];
2849 12796 : for (i = imin; i < l; i++)
2850 : {
2851 9184 : GEN c = gel(v,i);
2852 9184 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2853 9184 : gel(V,i) = c;
2854 : }
2855 3612 : return V;
2856 : }
2857 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2858 : GEN
2859 3241 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2860 3241 : { return do_QXQ_eval(v, 1, a, T); }
2861 :
2862 : GEN
2863 371 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2864 371 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2865 :
2866 : GEN
2867 2940 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2868 : {
2869 2940 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2870 : }
2871 :
2872 : GEN
2873 2680 : RgXQ_norm(GEN x, GEN T)
2874 : {
2875 : pari_sp av;
2876 2680 : long dx = degpol(x);
2877 : GEN L, y;
2878 2680 : if (degpol(T)==0) return gpowgs(x,0);
2879 2673 : av = avma; y = resultant(T, x);
2880 2673 : L = leading_coeff(T);
2881 2673 : if (gequal1(L) || !signe(x)) return y;
2882 0 : return gc_upto(av, gdiv(y, gpowgs(L, dx)));
2883 : }
2884 :
2885 : GEN
2886 476 : RgXQ_trace(GEN x, GEN T)
2887 : {
2888 476 : pari_sp av = avma;
2889 : GEN dT, z;
2890 : long n;
2891 476 : if (degpol(T)==0) return gmulgs(x,0);
2892 469 : dT = RgX_deriv(T); n = degpol(dT);
2893 469 : z = RgXQ_mul(x, dT, T);
2894 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2895 420 : return gc_upto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2896 : }
2897 :
2898 : GEN
2899 2788484 : RgX_blocks(GEN P, long n, long m)
2900 : {
2901 2788484 : GEN z = cgetg(m+1,t_VEC);
2902 2788484 : long i,j, k=2, l = lg(P);
2903 8559030 : for(i=1; i<=m; i++)
2904 : {
2905 5770546 : GEN zi = cgetg(n+2,t_POL);
2906 5770546 : zi[1] = P[1];
2907 5770546 : gel(z,i) = zi;
2908 17024680 : for(j=2; j<n+2; j++)
2909 11254134 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2910 5770546 : zi = RgX_renormalize_lg(zi, n+2);
2911 : }
2912 2788484 : return z;
2913 : }
2914 :
2915 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2916 : void
2917 49770473 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2918 : {
2919 49770473 : long n = degpol(p), v = varn(p), n0, n1, i;
2920 : GEN p0, p1;
2921 :
2922 49771779 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2923 :
2924 49678495 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2925 49678495 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2926 49681267 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2927 146818889 : for (i=0; i<n1; i++)
2928 : {
2929 97137177 : p0[2+i] = p[2+(i<<1)];
2930 97137177 : p1[2+i] = p[3+(i<<1)];
2931 : }
2932 49681712 : if (n1 != n0)
2933 18027402 : p0[2+i] = p[2+(i<<1)];
2934 49681712 : *pe = normalizepol(p0);
2935 49683921 : *po = normalizepol(p1);
2936 : }
2937 :
2938 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2939 : GEN
2940 43631 : RgX_splitting(GEN p, long k)
2941 : {
2942 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
2943 : GEN r;
2944 :
2945 43631 : m = n/k;
2946 43631 : r = cgetg(k+1,t_VEC);
2947 234059 : for(i=1; i<=k; i++)
2948 : {
2949 190428 : gel(r,i) = cgetg(m+3, t_POL);
2950 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2951 : }
2952 571991 : for (j=1, i=0, l=2; i<=n; i++)
2953 : {
2954 528360 : gmael(r,j,l) = gel(p,2+i);
2955 528360 : if (j==k) { j=1; l++; } else j++;
2956 : }
2957 234059 : for(i=1; i<=k; i++)
2958 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2959 43631 : return r;
2960 : }
2961 :
2962 : /*******************************************************************/
2963 : /* */
2964 : /* Kronecker form */
2965 : /* */
2966 : /*******************************************************************/
2967 :
2968 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2969 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2970 : * normalized coefficients */
2971 : GEN
2972 185965 : Kronecker_to_mod(GEN z, GEN T)
2973 : {
2974 185965 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2975 185965 : GEN x, t = cgetg(N,t_POL);
2976 185965 : t[1] = T[1];
2977 185965 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2978 185965 : x[1] = z[1];
2979 185965 : T = RgX_copy(T);
2980 1398524 : for (i=2; i<lx+2; i++, z+= N-2)
2981 : {
2982 5758508 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2983 1212559 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2984 : }
2985 185965 : N = (l-2) % (N-2) + 2;
2986 471775 : for (j=2; j<N; j++) t[j] = z[j];
2987 185965 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2988 185965 : return normalizepol_lg(x, i+1);
2989 : }
2990 :
2991 : /*******************************************************************/
2992 : /* */
2993 : /* Domain detection */
2994 : /* */
2995 : /*******************************************************************/
2996 :
2997 : static GEN
2998 609812 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2999 : {
3000 609812 : pari_sp av = avma;
3001 : GEN r;
3002 609812 : if (lgefint(p) == 3)
3003 : {
3004 558456 : ulong pp = uel(p, 2);
3005 558456 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
3006 : RgX_to_Flx(y, pp), pp));
3007 : }
3008 : else
3009 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3010 609812 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3011 544614 : return gc_upto(av, FpX_to_mod(r, p));
3012 : }
3013 :
3014 : static GEN
3015 2079 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3016 : {
3017 2079 : pari_sp av = avma;
3018 : long dT;
3019 : GEN kx, ky, r;
3020 2079 : GEN T = RgX_to_FpX(pol, p);
3021 2079 : if (signe(T)==0) pari_err_OP("*", x, y);
3022 2072 : dT = degpol(T);
3023 2072 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3024 2072 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3025 2072 : r = FpX_mul(kx, ky, p);
3026 2072 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3027 1463 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3028 : }
3029 :
3030 : static GEN
3031 454903 : RgX_liftred(GEN x, GEN T)
3032 454903 : { return RgXQX_red(liftpol_shallow(x), T); }
3033 :
3034 : static GEN
3035 170888 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3036 : {
3037 170888 : pari_sp av = avma;
3038 170888 : long dT = degpol(T);
3039 170888 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3040 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3041 170888 : return gc_upto(av, Kronecker_to_mod(r, T));
3042 : }
3043 :
3044 : static GEN
3045 1407 : RgX_sqr_FpX(GEN x, GEN p)
3046 : {
3047 1407 : pari_sp av = avma;
3048 : GEN r;
3049 1407 : if (lgefint(p) == 3)
3050 : {
3051 1203 : ulong pp = uel(p, 2);
3052 1203 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3053 : }
3054 : else
3055 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3056 1407 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3057 1274 : return gc_upto(av, FpX_to_mod(r, p));
3058 : }
3059 :
3060 : static GEN
3061 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3062 : {
3063 196 : pari_sp av = avma;
3064 : long dT;
3065 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3066 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3067 189 : dT = degpol(T);
3068 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3069 189 : r = FpX_sqr(kx, p);
3070 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3071 189 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3072 : }
3073 :
3074 : static GEN
3075 13425 : RgX_sqr_QXQX(GEN x, GEN T)
3076 : {
3077 13425 : pari_sp av = avma;
3078 13425 : long dT = degpol(T);
3079 13425 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3080 13425 : return gc_upto(av, Kronecker_to_mod(r, T));
3081 : }
3082 :
3083 : static GEN
3084 68327 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3085 : {
3086 68327 : pari_sp av = avma;
3087 : GEN r;
3088 68327 : if (lgefint(p) == 3)
3089 : {
3090 51691 : ulong pp = uel(p, 2);
3091 51691 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3092 : RgX_to_Flx(y, pp), pp));
3093 : }
3094 : else
3095 16636 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3096 68327 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3097 31031 : return gc_upto(av, FpX_to_mod(r, p));
3098 : }
3099 :
3100 : static GEN
3101 49851 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3102 : {
3103 49851 : pari_sp av = avma;
3104 : GEN r;
3105 49851 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3106 49851 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3107 : }
3108 : static GEN
3109 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3110 : {
3111 70 : pari_sp av = avma;
3112 : GEN r;
3113 70 : GEN T = RgX_to_FpX(pol, p);
3114 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3115 63 : if (lgefint(p) == 3)
3116 : {
3117 55 : ulong pp = uel(p, 2);
3118 55 : GEN Tp = ZX_to_Flx(T, pp);
3119 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3120 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3121 : }
3122 : else
3123 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3124 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3125 56 : return gc_upto(av, FpXQX_to_mod(r, T, p));
3126 : }
3127 :
3128 : static GEN
3129 119205912 : RgX_mul_fast(GEN x, GEN y)
3130 : {
3131 : GEN p, pol;
3132 : long pa;
3133 119205912 : long t = RgX_type2(x,y, &p,&pol,&pa);
3134 119206075 : switch(t)
3135 : {
3136 65247677 : case t_INT: return ZX_mul(x,y);
3137 3360387 : case t_FRAC: return QX_mul(x,y);
3138 105127 : case t_FFELT: return FFX_mul(x, y, pol);
3139 609812 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3140 170888 : case RgX_type_code(t_POLMOD, t_INT):
3141 : case RgX_type_code(t_POLMOD, t_FRAC):
3142 170888 : return RgX_mul_QXQX(x, y, pol);
3143 2069 : case RgX_type_code(t_POLMOD, t_INTMOD):
3144 2069 : return RgX_mul_FpXQX(x, y, pol, p);
3145 49710115 : default: return NULL;
3146 : }
3147 : }
3148 : static GEN
3149 1311361 : RgX_sqr_fast(GEN x)
3150 : {
3151 : GEN p, pol;
3152 : long pa;
3153 1311361 : long t = RgX_type(x,&p,&pol,&pa);
3154 1311360 : switch(t)
3155 : {
3156 1173349 : case t_INT: return ZX_sqr(x);
3157 82524 : case t_FRAC: return QX_sqr(x);
3158 2499 : case t_FFELT: return FFX_sqr(x, pol);
3159 1407 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3160 13425 : case RgX_type_code(t_POLMOD, t_INT):
3161 : case RgX_type_code(t_POLMOD, t_FRAC):
3162 13425 : return RgX_sqr_QXQX(x, pol);
3163 195 : case RgX_type_code(t_POLMOD, t_INTMOD):
3164 195 : return RgX_sqr_FpXQX(x, pol, p);
3165 37961 : default: return NULL;
3166 : }
3167 : }
3168 :
3169 : static GEN
3170 14699691 : RgX_rem_fast(GEN x, GEN y)
3171 : {
3172 : GEN p, pol;
3173 : long pa;
3174 14699691 : long t = RgX_type2(x,y, &p,&pol,&pa);
3175 14699828 : switch(t)
3176 : {
3177 3509594 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3178 1837073 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3179 84 : case t_FFELT: return FFX_rem(x, y, pol);
3180 68327 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3181 49851 : case RgX_type_code(t_POLMOD, t_INT):
3182 : case RgX_type_code(t_POLMOD, t_FRAC):
3183 49851 : return RgX_rem_QXQX(x, y, pol);
3184 61 : case RgX_type_code(t_POLMOD, t_INTMOD):
3185 61 : return RgX_rem_FpXQX(x, y, pol, p);
3186 9234838 : default: return NULL;
3187 : }
3188 : }
3189 :
3190 : GEN
3191 105467026 : RgX_mul(GEN x, GEN y)
3192 : {
3193 105467026 : GEN z = RgX_mul_fast(x,y);
3194 105467103 : if (!z) z = RgX_mul_i(x,y);
3195 105466711 : return z;
3196 : }
3197 :
3198 : GEN
3199 1299209 : RgX_sqr(GEN x)
3200 : {
3201 1299209 : GEN z = RgX_sqr_fast(x);
3202 1299204 : if (!z) z = RgX_sqr_i(x);
3203 1299204 : return z;
3204 : }
3205 :
3206 : GEN
3207 14699699 : RgX_rem(GEN x, GEN y)
3208 : {
3209 14699699 : GEN z = RgX_rem_fast(x, y);
3210 14699811 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3211 14699801 : return z;
3212 : }
3213 :
3214 : static GEN
3215 0 : _RgX_mul(void* E, GEN x, GEN y)
3216 0 : { (void) E; return RgX_mul(x, y); }
3217 :
3218 : GEN
3219 0 : RgXV_prod(GEN V)
3220 0 : { return gen_product(V, NULL, &_RgX_mul); }
3221 :
3222 : GEN
3223 11067686 : RgXn_mul(GEN f, GEN g, long n)
3224 : {
3225 11067686 : pari_sp av = avma;
3226 11067686 : GEN h = RgX_mul_fast(f,g);
3227 11067686 : if (!h) return RgXn_mul2(f,g,n);
3228 1783815 : if (degpol(h) < n) return h;
3229 383270 : return gc_GEN(av, RgXn_red_shallow(h, n));
3230 : }
3231 :
3232 : GEN
3233 12152 : RgXn_sqr(GEN f, long n)
3234 : {
3235 12152 : pari_sp av = avma;
3236 12152 : GEN g = RgX_sqr_fast(f);
3237 12152 : if (!g) return RgXn_sqr2(f,n);
3238 11837 : if (degpol(g) < n) return g;
3239 10640 : return gc_GEN(av, RgXn_red_shallow(g, n));
3240 : }
3241 :
3242 : /* (f*g) \/ x^n */
3243 : GEN
3244 2671189 : RgX_mulhigh_i(GEN f, GEN g, long n)
3245 : {
3246 2671189 : GEN h = RgX_mul_fast(f,g);
3247 2671189 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3248 : }
3249 :
3250 : /* (f*g) \/ x^n */
3251 : GEN
3252 0 : RgX_sqrhigh_i(GEN f, long n)
3253 : {
3254 0 : GEN h = RgX_sqr_fast(f);
3255 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3256 : }
|
