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 19168425 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19168425 : long pold=1, rold=n*(d-1);
35 91184352 : for(p=2; p<=d; p++)
36 : {
37 72015927 : r = m*(p-1) + n*((d-1)/p);
38 72015927 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19168425 : return pold;
41 : }
42 :
43 : static GEN
44 13123292 : 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 13123292 : pari_sp av = avma;
48 : long i;
49 13123292 : GEN z = cmul(E,P,a,ff->one(E));
50 13084861 : if (!z) z = gen_0;
51 67564933 : for (i=1; i<=n; i++)
52 : {
53 54472146 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 54487481 : if (t) {
55 41669933 : z = ff->add(E, z, t);
56 41653736 : if (gc_needed(av,2)) z = gerepileupto(av, z);
57 : }
58 : }
59 13092787 : 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 11671072 : 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 11671072 : pari_sp av = avma;
72 11671072 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 11671072 : if (d < 0) return ff->zero(E);
76 11107971 : if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1007163 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1007163 : 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 1007163 : d -= l;
84 1007163 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2015237 : while (d >= l-1)
86 : {
87 1006971 : d -= l-1;
88 1006971 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1007020 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1006969 : if (gc_needed(av,2))
91 105 : z = gerepileupto(av, z);
92 : }
93 1008266 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1008288 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1008284 : return gerepileupto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 727517 : 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 727517 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 727517 : if (d < 0) return ff->zero(E);
106 726586 : rtd = (long) sqrt((double)d);
107 726586 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 726584 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 726590 : return gerepileupto(av, z);
110 : }
111 :
112 : static GEN
113 1193647 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 15692651 : _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 1040043 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 339860 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 1217201 : _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 225890 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 225890 : pari_sp av = avma;
163 225890 : long i, d = degpol(P);
164 : GEN s;
165 225890 : if (signe(P)==0) return pol_0(varn(P));
166 225890 : s = gel(P, d+2);
167 225890 : if (d == 0) return gcopy(s);
168 1078854 : for (i = d-1; i >= 0; i--)
169 : {
170 856310 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
171 856310 : if (gc_needed(av,1))
172 : {
173 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
174 13 : s = gerepileupto(av, s);
175 : }
176 : }
177 222544 : return gerepileupto(av, s);
178 : }
179 :
180 : GEN
181 1568 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
182 : {
183 1568 : pari_sp av = avma;
184 1568 : long i, d = degpol(P), v = varn(A);
185 : GEN s;
186 1568 : if (signe(P)==0) return pol_0(v);
187 1568 : if (d == 0) return scalarpol(gel(P, d+2), v);
188 1190 : s = scalarpol_shallow(gel(P, d+2), v);
189 4851 : for (i = d-1; i >= 0; i--)
190 : {
191 3661 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
192 3661 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
193 3661 : if (gc_needed(av,1))
194 : {
195 0 : if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
196 0 : s = gerepileupto(av, s);
197 : }
198 : }
199 1190 : return gerepileupto(av, s);
200 : }
201 :
202 : const struct bb_algebra *
203 154026 : get_Rg_algebra(void)
204 : {
205 154026 : return &Rg_algebra;
206 : }
207 :
208 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
209 :
210 : static GEN
211 18004 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
212 : {
213 : (void) E;
214 18004 : return RgX_divrem(x, y, r);
215 : }
216 :
217 : GEN
218 4487 : RgX_digits(GEN x, GEN T)
219 : {
220 4487 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
221 4487 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
222 : }
223 :
224 : /*******************************************************************/
225 : /* */
226 : /* RgX */
227 : /* */
228 : /*******************************************************************/
229 :
230 : long
231 22809689 : RgX_equal(GEN x, GEN y)
232 : {
233 22809689 : long i = lg(x);
234 :
235 22809689 : if (i != lg(y)) return 0;
236 96540029 : for (i--; i > 1; i--)
237 74034648 : if (!gequal(gel(x,i),gel(y,i))) return 0;
238 22505381 : return 1;
239 : }
240 :
241 : /* Returns 1 in the base ring over which x is defined */
242 : /* HACK: this also works for t_SER */
243 : GEN
244 2148740 : Rg_get_1(GEN x)
245 : {
246 : GEN p, T;
247 2148740 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
248 2148739 : if (RgX_type_is_composite(tx))
249 39575 : RgX_type_decode(tx, &i /*junk*/, &tx);
250 2148739 : switch(tx)
251 : {
252 364 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
253 7 : case t_PADIC: return cvtop(gen_1, p, lx);
254 728 : case t_FFELT: return FF_1(T);
255 2147640 : default: return gen_1;
256 : }
257 : }
258 : /* Returns 0 in the base ring over which x is defined */
259 : /* HACK: this also works for t_SER */
260 : GEN
261 6974639 : Rg_get_0(GEN x)
262 : {
263 : GEN p, T;
264 6974639 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
265 6974639 : if (RgX_type_is_composite(tx))
266 52486 : RgX_type_decode(tx, &i /*junk*/, &tx);
267 6974639 : switch(tx)
268 : {
269 441 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
270 42 : case t_PADIC: return zeropadic(p, lx);
271 231 : case t_FFELT: return FF_zero(T);
272 6973925 : default: return gen_0;
273 : }
274 : }
275 :
276 : GEN
277 5292 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
278 : {
279 5292 : long i, n = degpol(P);
280 : GEN z, dz, dP;
281 5292 : if (n < 0) return gen_0;
282 5292 : P = Q_remove_denom(P, &dP);
283 5292 : z = gel(P,2); if (n == 0) return icopy(z);
284 3136 : if (dV) z = mulii(dV, z); /* V[1] = dV */
285 3136 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
286 6146 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
287 3136 : dz = mul_denom(dP, dV);
288 3136 : return dz? RgX_Rg_div(z, dz): z;
289 : }
290 :
291 : /* Return P(h * x), not memory clean */
292 : GEN
293 18281 : RgX_unscale(GEN P, GEN h)
294 : {
295 18281 : long i, l = lg(P);
296 18281 : GEN hi = gen_1, Q = cgetg(l, t_POL);
297 18281 : Q[1] = P[1];
298 18281 : if (l == 2) return Q;
299 18253 : gel(Q,2) = gcopy(gel(P,2));
300 40120 : for (i=3; i<l; i++)
301 : {
302 21867 : hi = gmul(hi,h);
303 21867 : gel(Q,i) = gmul(gel(P,i), hi);
304 : }
305 18253 : return Q;
306 : }
307 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
308 : GEN
309 1448728 : ZX_z_unscale(GEN P, long h)
310 : {
311 1448728 : long i, l = lg(P);
312 1448728 : GEN Q = cgetg(l, t_POL);
313 1448728 : Q[1] = P[1];
314 1448728 : if (l == 2) return Q;
315 1360562 : gel(Q,2) = gel(P,2);
316 1360562 : if (l == 3) return Q;
317 1334144 : if (h == -1)
318 230178 : for (i = 3; i < l; i++)
319 : {
320 191012 : gel(Q,i) = negi(gel(P,i));
321 191015 : if (++i == l) break;
322 141598 : gel(Q,i) = gel(P,i);
323 : }
324 : else
325 : {
326 : GEN hi;
327 1245564 : gel(Q,3) = mulis(gel(P,3), h);
328 1245564 : hi = sqrs(h);
329 5357542 : for (i = 4; i < l; i++)
330 : {
331 4111977 : gel(Q,i) = mulii(gel(P,i), hi);
332 4111976 : if (i != l-1) hi = mulis(hi,h);
333 : }
334 : }
335 1334148 : return Q;
336 : }
337 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
338 : GEN
339 998550 : ZX_unscale(GEN P, GEN h)
340 : {
341 : long i, l;
342 : GEN Q, hi;
343 998550 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
344 804 : l = lg(P); Q = cgetg(l, t_POL);
345 804 : Q[1] = P[1];
346 804 : if (l == 2) return Q;
347 804 : gel(Q,2) = gel(P,2);
348 804 : if (l == 3) return Q;
349 804 : hi = h;
350 804 : gel(Q,3) = mulii(gel(P,3), hi);
351 2531 : for (i = 4; i < l; i++)
352 : {
353 1727 : hi = mulii(hi,h);
354 1727 : gel(Q,i) = mulii(gel(P,i), hi);
355 : }
356 804 : return Q;
357 : }
358 : /* P a ZX. Return P(x << n), not memory clean */
359 : GEN
360 216659 : ZX_unscale2n(GEN P, long n)
361 : {
362 216659 : long i, ni = n, l = lg(P);
363 216659 : GEN Q = cgetg(l, t_POL);
364 216661 : Q[1] = P[1];
365 216661 : if (l == 2) return Q;
366 216661 : gel(Q,2) = gel(P,2);
367 216661 : if (l == 3) return Q;
368 216661 : gel(Q,3) = shifti(gel(P,3), ni);
369 799049 : for (i=4; i<l; i++)
370 : {
371 582426 : ni += n;
372 582426 : gel(Q,i) = shifti(gel(P,i), ni);
373 : }
374 216623 : return Q;
375 : }
376 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
377 : GEN
378 12019 : ZX_unscale_div(GEN P, GEN h)
379 : {
380 12019 : long i, l = lg(P);
381 12019 : GEN hi, Q = cgetg(l, t_POL);
382 12019 : Q[1] = P[1];
383 12019 : if (l == 2) return Q;
384 12019 : gel(Q,2) = diviiexact(gel(P,2), h);
385 12019 : if (l == 3) return Q;
386 12019 : gel(Q,3) = gel(P,3);
387 12019 : if (l == 4) return Q;
388 12019 : hi = h;
389 12019 : gel(Q,4) = mulii(gel(P,4), hi);
390 61187 : for (i=5; i<l; i++)
391 : {
392 49168 : hi = mulii(hi,h);
393 49168 : gel(Q,i) = mulii(gel(P,i), hi);
394 : }
395 12019 : return Q;
396 : }
397 : /* P(h*X) / h^k, assuming the result is a ZX */
398 : GEN
399 966 : ZX_unscale_divpow(GEN P, GEN h, long k)
400 : {
401 966 : long i, j, l = lg(P);
402 966 : GEN H, Q = cgetg(l, t_POL);
403 966 : Q[1] = P[1]; if (l == 2) return Q;
404 966 : H = gpowers(h, maxss(k, l - 3 - k));
405 3864 : for (i = 2, j = k+1; j > 1 && i < l; i++)
406 2898 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
407 966 : if (i == l) return Q;
408 966 : gel(Q, i) = gel(P, i); i++;
409 3542 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
410 966 : return Q;
411 : }
412 :
413 : GEN
414 5558 : RgXV_unscale(GEN x, GEN h)
415 : {
416 5558 : if (isint1(h)) return RgX_copy(x);
417 16202 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
418 : }
419 :
420 : /* Return h^degpol(P) P(x / h), not memory clean */
421 : GEN
422 2439967 : RgX_rescale(GEN P, GEN h)
423 : {
424 2439967 : long i, l = lg(P);
425 2439967 : GEN Q = cgetg(l,t_POL), hi = h;
426 2439948 : gel(Q,l-1) = gel(P,l-1);
427 7641824 : for (i=l-2; i>=2; i--)
428 : {
429 7639121 : gel(Q,i) = gmul(gel(P,i), hi);
430 7638799 : if (i == 2) break;
431 5201333 : hi = gmul(hi,h);
432 : }
433 2440169 : Q[1] = P[1]; return Q;
434 : }
435 :
436 : GEN
437 2289 : RgXV_rescale(GEN x, GEN h)
438 : {
439 2289 : if (isint1(h)) return RgX_copy(x);
440 15862 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
441 : }
442 :
443 : /* A(X^d) --> A(X) */
444 : GEN
445 724325 : RgX_deflate(GEN x0, long d)
446 : {
447 : GEN z, y, x;
448 724325 : long i,id, dy, dx = degpol(x0);
449 724323 : if (d == 1 || dx <= 0) return leafcopy(x0);
450 331644 : dy = dx/d;
451 331644 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
452 331644 : z = y + 2;
453 331644 : x = x0+ 2;
454 1367380 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
455 331644 : return y;
456 : }
457 :
458 : GEN
459 3584 : RgX_homogenize(GEN P, long v)
460 : {
461 : long i, l, d;
462 3584 : GEN Q = cgetg_copy(P, &l);
463 3584 : Q[1] = P[1]; d = l-3;
464 31794 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
465 3584 : return Q;
466 : }
467 :
468 : /* F a t_RFRAC */
469 : long
470 98 : rfrac_deflate_order(GEN F)
471 : {
472 98 : GEN N = gel(F,1), D = gel(F,2);
473 98 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
474 98 : if (m == 1) return 1;
475 35 : if (typ(N) == t_POL && varn(N) == varn(D))
476 21 : m = cgcd(m, RgX_deflate_order(N));
477 35 : return m;
478 : }
479 : /* F a t_RFRAC */
480 : GEN
481 98 : rfrac_deflate_max(GEN F, long *m)
482 : {
483 98 : *m = rfrac_deflate_order(F);
484 98 : return rfrac_deflate(F, *m);
485 : }
486 : /* F a t_RFRAC */
487 : GEN
488 98 : rfrac_deflate(GEN F, long m)
489 : {
490 98 : GEN N = gel(F,1), D = gel(F,2);
491 98 : if (m == 1) return F;
492 35 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
493 35 : D = RgX_deflate(D, m); return mkrfrac(N, D);
494 : }
495 :
496 : /* return x0(X^d) */
497 : GEN
498 794215 : RgX_inflate(GEN x0, long d)
499 : {
500 794215 : long i, id, dy, dx = degpol(x0);
501 794215 : GEN x = x0 + 2, z, y;
502 794215 : if (dx <= 0) return leafcopy(x0);
503 786613 : dy = dx*d;
504 786613 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
505 786612 : z = y + 2;
506 25948517 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
507 11394440 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
508 786612 : return y;
509 : }
510 :
511 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
512 : static GEN
513 4588756 : RgX_translate_basecase(GEN P, GEN c)
514 : {
515 4588756 : pari_sp av = avma;
516 : GEN Q, R;
517 : long i, k, n;
518 :
519 4588756 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
520 4586757 : Q = leafcopy(P);
521 4586754 : R = Q+2; n = degpol(P);
522 4586757 : if (isint1(c))
523 : {
524 26808 : for (i=1; i<=n; i++)
525 : {
526 138991 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
527 22910 : if (gc_needed(av,2))
528 : {
529 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
530 0 : Q = gerepilecopy(av, Q); R = Q+2;
531 : }
532 : }
533 : }
534 4582859 : else if (isintm1(c))
535 : {
536 914875 : for (i=1; i<=n; i++)
537 : {
538 2909487 : for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
539 694385 : if (gc_needed(av,2))
540 : {
541 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
542 0 : Q = gerepilecopy(av, Q); R = Q+2;
543 : }
544 : }
545 : }
546 : else
547 : {
548 15412746 : for (i=1; i<=n; i++)
549 : {
550 38492525 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
551 11050398 : if (gc_needed(av,2))
552 : {
553 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
554 0 : Q = gerepilecopy(av, Q); R = Q+2;
555 : }
556 : }
557 : }
558 4586094 : return gerepilecopy(av, Q);
559 : }
560 : GEN
561 4589188 : RgX_translate(GEN P, GEN c)
562 : {
563 4589188 : pari_sp av = avma;
564 4589188 : long n = degpol(P);
565 4589181 : if (n < 40)
566 4588750 : return RgX_translate_basecase(P, c);
567 : else
568 : {
569 431 : long d = n >> 1;
570 431 : GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
571 432 : GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
572 432 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
573 432 : return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
574 : }
575 : }
576 :
577 : /* P(ax + b) */
578 : GEN
579 0 : RgX_affine(GEN P, GEN a, GEN b)
580 : {
581 0 : if (signe(b)) P = RgX_translate(P, b);
582 0 : return RgX_unscale(P, a);
583 : }
584 :
585 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
586 : GEN
587 22537 : RgXQX_translate(GEN P, GEN c, GEN T)
588 : {
589 22537 : pari_sp av = avma;
590 : GEN Q, R;
591 : long i, k, n;
592 :
593 22537 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
594 22229 : Q = leafcopy(P);
595 22229 : R = Q+2; n = degpol(P);
596 79092 : for (i=1; i<=n; i++)
597 : {
598 253848 : for (k=n-i; k<n; k++)
599 : {
600 196985 : pari_sp av2 = avma;
601 196985 : gel(R,k) = gerepileupto(av2,
602 196985 : RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
603 : }
604 56863 : if (gc_needed(av,2))
605 : {
606 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
607 0 : Q = gerepilecopy(av, Q); R = Q+2;
608 : }
609 : }
610 22229 : return gerepilecopy(av, Q);
611 : }
612 :
613 : /********************************************************************/
614 : /** **/
615 : /** CONVERSIONS **/
616 : /** (not memory clean) **/
617 : /** **/
618 : /********************************************************************/
619 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
620 : * but everything else is */
621 : static GEN
622 160504 : QXQ_to_mod(GEN x, GEN T)
623 : {
624 : long d;
625 160504 : switch(typ(x))
626 : {
627 59807 : case t_INT: return icopy(x);
628 1057 : case t_FRAC: return gcopy(x);
629 99642 : case t_POL:
630 99642 : d = degpol(x);
631 99642 : if (d < 0) return gen_0;
632 97374 : if (d == 0) return gcopy(gel(x,2));
633 95213 : return mkpolmod(RgX_copy(x), T);
634 0 : default: pari_err_TYPE("QXQ_to_mod",x);
635 : return NULL;/* LCOV_EXCL_LINE */
636 : }
637 : }
638 : /* pure shallow version */
639 : GEN
640 695335 : QXQ_to_mod_shallow(GEN x, GEN T)
641 : {
642 : long d;
643 695335 : switch(typ(x))
644 : {
645 444964 : case t_INT:
646 444964 : case t_FRAC: return x;
647 250371 : case t_POL:
648 250371 : d = degpol(x);
649 250371 : if (d < 0) return gen_0;
650 205745 : if (d == 0) return gel(x,2);
651 191560 : return mkpolmod(x, T);
652 0 : default: pari_err_TYPE("QXQ_to_mod",x);
653 : return NULL;/* LCOV_EXCL_LINE */
654 : }
655 : }
656 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
657 : * Not memory clean because T not copied, but everything else is */
658 : static GEN
659 37170 : QXQX_to_mod(GEN z, GEN T)
660 : {
661 37170 : long i,l = lg(z);
662 37170 : GEN x = cgetg(l,t_POL);
663 185750 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
664 37172 : x[1] = z[1]; return normalizepol_lg(x,l);
665 : }
666 : /* pure shallow version */
667 : GEN
668 155611 : QXQX_to_mod_shallow(GEN z, GEN T)
669 : {
670 155611 : long i,l = lg(z);
671 155611 : GEN x = cgetg(l,t_POL);
672 770222 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
673 155611 : x[1] = z[1]; return normalizepol_lg(x,l);
674 : }
675 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
676 : GEN
677 12488 : QXQXV_to_mod(GEN V, GEN T)
678 : {
679 12488 : long i, l = lg(V);
680 12488 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
681 49658 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
682 12488 : return z;
683 : }
684 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
685 : GEN
686 11743 : QXQV_to_mod(GEN V, GEN T)
687 : {
688 11743 : long i, l = lg(V);
689 11743 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
690 23667 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
691 11743 : return z;
692 : }
693 :
694 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
695 : GEN
696 13384 : QXQC_to_mod_shallow(GEN V, GEN T)
697 : {
698 13384 : long i, l = lg(V);
699 13384 : GEN z = cgetg(l, t_COL);
700 94108 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
701 13384 : return z;
702 : }
703 :
704 : GEN
705 6244 : QXQM_to_mod_shallow(GEN V, GEN T)
706 : {
707 6244 : long i, l = lg(V);
708 6244 : GEN z = cgetg(l, t_MAT);
709 19628 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
710 6244 : return z;
711 : }
712 :
713 : GEN
714 6083464 : RgX_renormalize_lg(GEN x, long lx)
715 : {
716 : long i;
717 9081422 : for (i = lx-1; i>1; i--)
718 8635069 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
719 6083464 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
720 6083464 : setlg(x, i+1); setsigne(x, i != 1); return x;
721 : }
722 :
723 : GEN
724 1466982 : RgV_to_RgX(GEN x, long v)
725 : {
726 1466982 : long i, k = lg(x);
727 : GEN p;
728 :
729 4109130 : while (--k && gequal0(gel(x,k)));
730 1466982 : if (!k) return pol_0(v);
731 1460105 : i = k+2; p = cgetg(i,t_POL);
732 1460109 : p[1] = evalsigne(1) | evalvarn(v);
733 10122172 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
734 1460109 : return p;
735 : }
736 : GEN
737 196780 : RgV_to_RgX_reverse(GEN x, long v)
738 : {
739 196780 : long j, k, l = lg(x);
740 : GEN p;
741 :
742 198259 : for (k = 1; k < l; k++)
743 198259 : if (!gequal0(gel(x,k))) break;
744 196780 : if (k == l) return pol_0(v);
745 196780 : k -= 1;
746 196780 : l -= k;
747 196780 : x += k;
748 196780 : p = cgetg(l+1,t_POL);
749 196780 : p[1] = evalsigne(1) | evalvarn(v);
750 1083573 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
751 196780 : return p;
752 : }
753 :
754 : /* return the (N-dimensional) vector of coeffs of p */
755 : GEN
756 12285750 : RgX_to_RgC(GEN x, long N)
757 : {
758 : long i, l;
759 : GEN z;
760 12285750 : l = lg(x)-1; x++;
761 12285750 : if (l > N+1) l = N+1; /* truncate higher degree terms */
762 12285750 : z = cgetg(N+1,t_COL);
763 76564976 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
764 22956846 : for ( ; i<=N; i++) gel(z,i) = gen_0;
765 12285876 : return z;
766 : }
767 : GEN
768 1246164 : Rg_to_RgC(GEN x, long N)
769 : {
770 1246164 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
771 : }
772 :
773 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
774 : GEN
775 301373 : RgM_to_RgXV(GEN x, long v)
776 1285901 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
777 : GEN
778 7091 : RgM_to_RgXV_reverse(GEN x, long v)
779 28364 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
780 :
781 : /* matrix whose entries are given by the coeffs of the polynomials in
782 : * vector v (considered as degree n-1 polynomials) */
783 : GEN
784 287878 : RgV_to_RgM(GEN x, long n)
785 1531835 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
786 :
787 : GEN
788 43650 : RgXV_to_RgM(GEN x, long n)
789 255310 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
790 :
791 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
792 : GEN
793 18035 : RgM_to_RgXX(GEN x, long v,long w)
794 : {
795 18035 : long j, lx = lg(x);
796 18035 : GEN y = cgetg(lx+1, t_POL);
797 18035 : y[1] = evalsigne(1) | evalvarn(v);
798 18035 : y++;
799 103433 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
800 18035 : return normalizepol_lg(--y, lx+1);
801 : }
802 :
803 : /* matrix whose entries are given by the coeffs of the polynomial v in
804 : * two variables (considered as degree n-1 polynomials) */
805 : GEN
806 308 : RgXX_to_RgM(GEN v, long n)
807 : {
808 308 : long j, N = lg(v)-1;
809 308 : GEN y = cgetg(N, t_MAT);
810 1001 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
811 308 : return y;
812 : }
813 :
814 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
815 : GEN
816 29630 : RgXY_swapspec(GEN x, long n, long w, long nx)
817 : {
818 29630 : long j, ly = n+3;
819 29630 : GEN y = cgetg(ly, t_POL);
820 29630 : y[1] = evalsigne(1);
821 375099 : for (j=2; j<ly; j++)
822 : {
823 : long k;
824 345469 : GEN a = cgetg(nx+2,t_POL);
825 345469 : a[1] = evalsigne(1) | evalvarn(w);
826 1786001 : for (k=0; k<nx; k++)
827 : {
828 1440532 : GEN xk = gel(x,k);
829 1440532 : if (typ(xk)==t_POL)
830 1352082 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
831 : else
832 88450 : gel(a,k+2) = j==2 ? xk: gen_0;
833 : }
834 345469 : gel(y,j) = normalizepol_lg(a, nx+2);
835 : }
836 29630 : return normalizepol_lg(y,ly);
837 : }
838 :
839 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
840 : GEN
841 952 : RgXY_swap(GEN x, long n, long w)
842 : {
843 952 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
844 952 : setvarn(z, varn(x)); return z;
845 : }
846 :
847 : long
848 274 : RgXY_degreex(GEN b)
849 : {
850 274 : long deg = 0, i;
851 274 : if (!signe(b)) return -1;
852 1354 : for (i = 2; i < lg(b); ++i)
853 : {
854 1080 : GEN bi = gel(b, i);
855 1080 : if (typ(bi) == t_POL)
856 1016 : deg = maxss(deg, degpol(bi));
857 : }
858 274 : return deg;
859 : }
860 :
861 : GEN
862 16632 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
863 :
864 : /* return (x % X^n). Shallow */
865 : GEN
866 7782825 : RgXn_red_shallow(GEN a, long n)
867 : {
868 7782825 : long i, L = n+2, l = lg(a);
869 : GEN b;
870 7782825 : if (L >= l) return a; /* deg(x) < n */
871 5610358 : b = cgetg(L, t_POL); b[1] = a[1];
872 31326079 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
873 5610359 : return normalizepol_lg(b,L);
874 : }
875 :
876 : GEN
877 483 : RgXnV_red_shallow(GEN x, long n)
878 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
879 :
880 : /* return (x * X^n). Shallow */
881 : GEN
882 127243964 : RgX_shift_shallow(GEN a, long n)
883 : {
884 127243964 : long i, l = lg(a);
885 : GEN b;
886 127243964 : if (l == 2 || !n) return a;
887 86740814 : l += n;
888 86740814 : if (n < 0)
889 : {
890 51767185 : if (l <= 2) return pol_0(varn(a));
891 50315255 : b = cgetg(l, t_POL); b[1] = a[1];
892 50315380 : a -= n;
893 146074675 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
894 : } else {
895 34973629 : b = cgetg(l, t_POL); b[1] = a[1];
896 34974513 : a -= n; n += 2;
897 74266065 : for (i=2; i<n; i++) gel(b,i) = gen_0;
898 147052063 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
899 : }
900 85289893 : return b;
901 : }
902 : /* return (x * X^n). */
903 : GEN
904 1573346 : RgX_shift(GEN a, long n)
905 : {
906 1573346 : long i, l = lg(a);
907 : GEN b;
908 1573346 : if (l == 2 || !n) return RgX_copy(a);
909 1572982 : l += n;
910 1572982 : if (n < 0)
911 : {
912 847 : if (l <= 2) return pol_0(varn(a));
913 777 : b = cgetg(l, t_POL); b[1] = a[1];
914 777 : a -= n;
915 2947 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
916 : } else {
917 1572135 : b = cgetg(l, t_POL); b[1] = a[1];
918 1572135 : a -= n; n += 2;
919 3936789 : for (i=2; i<n; i++) gel(b,i) = gen_0;
920 4379263 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
921 : }
922 1572912 : return b;
923 : }
924 :
925 : GEN
926 317037 : RgX_rotate_shallow(GEN P, long k, long p)
927 : {
928 317037 : long i, l = lgpol(P);
929 : GEN r;
930 317037 : if (signe(P)==0)
931 1365 : return pol_0(varn(P));
932 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
933 2100644 : for(i=0; i<p; i++)
934 : {
935 1784972 : long s = 2+(i+k)%p;
936 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
937 : }
938 315672 : return RgX_renormalize(r);
939 : }
940 :
941 : GEN
942 2979364 : RgX_mulXn(GEN x, long d)
943 : {
944 : pari_sp av;
945 : GEN z;
946 : long v;
947 2979364 : if (d >= 0) return RgX_shift(x, d);
948 1467412 : d = -d;
949 1467412 : v = RgX_val(x);
950 1467412 : if (v >= d) return RgX_shift(x, -d);
951 1467405 : av = avma;
952 1467405 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
953 1467405 : return gerepileupto(av, z);
954 : }
955 :
956 : long
957 588 : RgXV_maxdegree(GEN x)
958 : {
959 588 : long d = -1, i, l = lg(x);
960 4494 : for (i = 1; i < l; i++)
961 3906 : d = maxss(d, degpol(gel(x,i)));
962 588 : return d;
963 : }
964 :
965 : long
966 3599618 : RgX_val(GEN x)
967 : {
968 3599618 : long i, lx = lg(x);
969 3599618 : if (lx == 2) return LONG_MAX;
970 4524864 : for (i = 2; i < lx; i++)
971 4524808 : if (!isexactzero(gel(x,i))) break;
972 3599464 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
973 3599408 : return i - 2;
974 : }
975 : long
976 62373049 : RgX_valrem(GEN x, GEN *Z)
977 : {
978 62373049 : long v, i, lx = lg(x);
979 62373049 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
980 104382522 : for (i = 2; i < lx; i++)
981 104390152 : if (!isexactzero(gel(x,i))) break;
982 : /* possible with nonrational zeros */
983 62391571 : if (i == lx)
984 : {
985 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
986 14 : return LONG_MAX;
987 : }
988 62391557 : v = i - 2;
989 62391557 : *Z = RgX_shift_shallow(x, -v);
990 62405384 : return v;
991 : }
992 : long
993 617940 : RgX_valrem_inexact(GEN x, GEN *Z)
994 : {
995 : long v;
996 617940 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
997 629737 : for (v = 0;; v++)
998 629737 : if (!gequal0(gel(x,2+v))) break;
999 617933 : if (Z) *Z = RgX_shift_shallow(x, -v);
1000 617933 : return v;
1001 : }
1002 :
1003 : GEN
1004 63819 : RgXQC_red(GEN x, GEN T)
1005 394653 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1006 :
1007 : GEN
1008 1204 : RgXQV_red(GEN x, GEN T)
1009 27769 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1010 :
1011 : GEN
1012 12243 : RgXQM_red(GEN x, GEN T)
1013 76062 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1014 :
1015 : GEN
1016 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1017 : {
1018 322 : return RgXQM_red(RgM_mul(P, Q), T);
1019 : }
1020 :
1021 : GEN
1022 447739 : RgXQX_red(GEN P, GEN T)
1023 : {
1024 447739 : long i, l = lg(P);
1025 447739 : GEN Q = cgetg(l, t_POL);
1026 447739 : Q[1] = P[1];
1027 2353601 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1028 447740 : return normalizepol_lg(Q, l);
1029 : }
1030 :
1031 : GEN
1032 574920 : RgX_deriv(GEN x)
1033 : {
1034 574920 : long i,lx = lg(x)-1;
1035 : GEN y;
1036 :
1037 574920 : if (lx<3) return pol_0(varn(x));
1038 573149 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1039 2814424 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1040 573132 : y[1] = x[1]; return normalizepol_lg(y,i);
1041 : }
1042 :
1043 : GEN
1044 1205225 : RgX_recipspec_shallow(GEN x, long l, long n)
1045 : {
1046 : long i;
1047 1205225 : GEN z = cgetg(n+2,t_POL);
1048 1205224 : z[1] = 0; z += 2;
1049 23385909 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1050 1420590 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1051 1205224 : return normalizepol_lg(z-2,n+2);
1052 : }
1053 :
1054 : GEN
1055 642330 : RgXn_recip_shallow(GEN P, long n)
1056 : {
1057 642330 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1058 642341 : setvarn(Q, varn(P));
1059 642341 : return Q;
1060 : }
1061 :
1062 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1063 : GEN
1064 16093 : RgX_recip(GEN x)
1065 : {
1066 : long lx, i, j;
1067 16093 : GEN y = cgetg_copy(x, &lx);
1068 139461 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1069 16093 : return normalizepol_lg(y,lx);
1070 : }
1071 : /* shallow version */
1072 : GEN
1073 57788 : RgX_recip_shallow(GEN x)
1074 : {
1075 : long lx, i, j;
1076 57788 : GEN y = cgetg_copy(x, &lx);
1077 346567 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1078 57788 : return normalizepol_lg(y,lx);
1079 : }
1080 :
1081 : GEN
1082 2314503 : RgX_recip_i(GEN x)
1083 : {
1084 : long lx, i, j;
1085 2314503 : GEN y = cgetg_copy(x, &lx);
1086 13366108 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1087 2314503 : return y;
1088 : }
1089 : /*******************************************************************/
1090 : /* */
1091 : /* ADDITION / SUBTRACTION */
1092 : /* */
1093 : /*******************************************************************/
1094 : /* same variable */
1095 : GEN
1096 74047340 : RgX_add(GEN x, GEN y)
1097 : {
1098 74047340 : long i, lx = lg(x), ly = lg(y);
1099 : GEN z;
1100 74047340 : if (ly <= lx) {
1101 64437803 : z = cgetg(lx,t_POL); z[1] = x[1];
1102 260711478 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1103 124384961 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1104 64386877 : z = normalizepol_lg(z, lx);
1105 : } else {
1106 9609537 : z = cgetg(ly,t_POL); z[1] = y[1];
1107 42297058 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1108 31311827 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1109 9610698 : z = normalizepol_lg(z, ly);
1110 : }
1111 74044100 : return z;
1112 : }
1113 : GEN
1114 43514673 : RgX_sub(GEN x, GEN y)
1115 : {
1116 43514673 : long i, lx = lg(x), ly = lg(y);
1117 : GEN z;
1118 43514673 : if (ly <= lx) {
1119 25570635 : z = cgetg(lx,t_POL); z[1] = x[1];
1120 144203168 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1121 47360092 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1122 25536334 : z = normalizepol_lg(z, lx);
1123 : } else {
1124 17944038 : z = cgetg(ly,t_POL); z[1] = y[1];
1125 72955312 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1126 42493830 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1127 17896451 : z = normalizepol_lg(z, ly);
1128 : }
1129 43509224 : return z;
1130 : }
1131 : GEN
1132 4590564 : RgX_neg(GEN x)
1133 : {
1134 4590564 : long i, lx = lg(x);
1135 4590564 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
1136 38304724 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
1137 4590378 : return y;
1138 : }
1139 :
1140 : GEN
1141 18254249 : RgX_Rg_add(GEN y, GEN x)
1142 : {
1143 : GEN z;
1144 18254249 : long lz = lg(y), i;
1145 18254249 : if (lz == 2) return scalarpol(x,varn(y));
1146 15226814 : z = cgetg(lz,t_POL); z[1] = y[1];
1147 15226773 : gel(z,2) = gadd(gel(y,2),x);
1148 59619252 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1149 : /* probably useless unless lz = 3, but cannot be skipped if y is
1150 : * an inexact 0 */
1151 15226887 : return normalizepol_lg(z,lz);
1152 : }
1153 : GEN
1154 65503 : RgX_Rg_add_shallow(GEN y, GEN x)
1155 : {
1156 : GEN z;
1157 65503 : long lz = lg(y), i;
1158 65503 : if (lz == 2) return scalarpol(x,varn(y));
1159 65503 : z = cgetg(lz,t_POL); z[1] = y[1];
1160 65503 : gel(z,2) = gadd(gel(y,2),x);
1161 131090 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1162 65503 : return normalizepol_lg(z,lz);
1163 : }
1164 : GEN
1165 184650 : RgX_Rg_sub(GEN y, GEN x)
1166 : {
1167 : GEN z;
1168 184650 : long lz = lg(y), i;
1169 184650 : if (lz == 2)
1170 : { /* scalarpol(gneg(x),varn(y)) optimized */
1171 2548 : long v = varn(y);
1172 2548 : if (isrationalzero(x)) return pol_0(v);
1173 14 : z = cgetg(3,t_POL);
1174 14 : z[1] = gequal0(x)? evalvarn(v)
1175 14 : : evalvarn(v) | evalsigne(1);
1176 14 : gel(z,2) = gneg(x); return z;
1177 : }
1178 182102 : z = cgetg(lz,t_POL); z[1] = y[1];
1179 182102 : gel(z,2) = gsub(gel(y,2),x);
1180 475394 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1181 182102 : return normalizepol_lg(z,lz);
1182 : }
1183 : GEN
1184 2975564 : Rg_RgX_sub(GEN x, GEN y)
1185 : {
1186 : GEN z;
1187 2975564 : long lz = lg(y), i;
1188 2975564 : if (lz == 2) return scalarpol(x,varn(y));
1189 2974486 : z = cgetg(lz,t_POL); z[1] = y[1];
1190 2974462 : gel(z,2) = gsub(x, gel(y,2));
1191 5458155 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1192 2974500 : return normalizepol_lg(z,lz);
1193 : }
1194 : /*******************************************************************/
1195 : /* */
1196 : /* KARATSUBA MULTIPLICATION */
1197 : /* */
1198 : /*******************************************************************/
1199 : #if 0
1200 : /* to debug Karatsuba-like routines */
1201 : GEN
1202 : zx_debug_spec(GEN x, long nx)
1203 : {
1204 : GEN z = cgetg(nx+2,t_POL);
1205 : long i;
1206 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1207 : z[1] = evalsigne(1); return z;
1208 : }
1209 :
1210 : GEN
1211 : RgX_debug_spec(GEN x, long nx)
1212 : {
1213 : GEN z = cgetg(nx+2,t_POL);
1214 : long i;
1215 : for (i=0; i<nx; i++) z[i+2] = x[i];
1216 : z[1] = evalsigne(1); return z;
1217 : }
1218 : #endif
1219 :
1220 : /* generic multiplication */
1221 : GEN
1222 7524661 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1223 : {
1224 : GEN z, t;
1225 : long i;
1226 7524661 : if (nx == ny) {
1227 1343413 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1228 4388873 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1229 1343396 : return normalizepol_lg(z, nx+2);
1230 : }
1231 6181248 : if (ny < nx) {
1232 5993153 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1233 23883323 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1234 14842792 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1235 5992556 : return normalizepol_lg(z, nx+2);
1236 : } else {
1237 188095 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1238 3683248 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1239 456278 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1240 188122 : return normalizepol_lg(z, ny+2);
1241 : }
1242 : }
1243 : GEN
1244 234041 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1245 : {
1246 : GEN z, t;
1247 : long i;
1248 234041 : if (nx == ny) {
1249 12810 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1250 2183524 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1251 12810 : return normalizepol_lg(z, nx+2);
1252 : }
1253 221231 : if (ny < nx) {
1254 219453 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1255 3794320 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1256 2526164 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1257 219453 : return normalizepol_lg(z, nx+2);
1258 : } else {
1259 1778 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1260 330904 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1261 12222 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1262 1778 : return normalizepol_lg(z, ny+2);
1263 : }
1264 : }
1265 :
1266 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1267 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1268 : * t_VECSMALL, to hold this, which can be left on stack: gerepile
1269 : * will not crash on it. The returned vector itself is not a proper GEN,
1270 : * we access the coefficients as x[i], i = 0..deg(x) */
1271 : static GEN
1272 48401874 : RgXspec_kill0(GEN x, long lx)
1273 : {
1274 48401874 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1275 : long i;
1276 171658326 : for (i=0; i <lx; i++)
1277 : {
1278 123256536 : GEN c = gel(x,i);
1279 123256536 : z[i] = (long)(isrationalzero(c)? NULL: c);
1280 : }
1281 48401790 : return z;
1282 : }
1283 :
1284 : INLINE GEN
1285 86240009 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1286 : {
1287 86240009 : pari_sp av = avma;
1288 86240009 : GEN s = NULL;
1289 : long i;
1290 :
1291 341877795 : for (i=a; i<b; i++)
1292 255650643 : if (gel(y,i) && gel(x,-i))
1293 : {
1294 177617808 : GEN t = gmul(gel(y,i), gel(x,-i));
1295 177611037 : s = s? gadd(s, t): t;
1296 : }
1297 86227152 : return s? gerepileupto(av, s): gen_0;
1298 : }
1299 :
1300 : /* assume nx >= ny > 0, return x * y * t^v */
1301 : static GEN
1302 17666057 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1303 : {
1304 : long i, lz, nz;
1305 : GEN z;
1306 :
1307 17666057 : x = RgXspec_kill0(x,nx);
1308 17666064 : y = RgXspec_kill0(y,ny);
1309 17666076 : lz = nx + ny + 1; nz = lz-2;
1310 17666076 : lz += v;
1311 17666076 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1312 29402076 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1313 49517329 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1314 33361934 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1315 31851283 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1316 17665498 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1317 : }
1318 :
1319 : /* return (x * X^d) + y. Assume d > 0 */
1320 : GEN
1321 7496228 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1322 : {
1323 : GEN x, y, xd, yd, zd;
1324 : long a, lz, nx, ny;
1325 :
1326 7496228 : if (!signe(x0)) return y0;
1327 7393764 : ny = lgpol(y0);
1328 7393760 : nx = lgpol(x0);
1329 7393932 : zd = (GEN)avma;
1330 7393932 : x = x0 + 2; y = y0 + 2; a = ny-d;
1331 7393932 : if (a <= 0)
1332 : {
1333 848801 : lz = nx+d+2;
1334 848801 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1335 2136063 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1336 848805 : x = zd + a;
1337 864483 : while (zd > x) gel(--zd,0) = gen_0;
1338 : }
1339 : else
1340 : {
1341 6545131 : xd = new_chunk(d); yd = y+d;
1342 6545128 : x = RgX_addspec_shallow(x,yd, nx,a);
1343 6545039 : lz = (a>nx)? ny+2: lg(x)+d;
1344 35200786 : x += 2; while (xd > x) *--zd = *--xd;
1345 : }
1346 19142376 : while (yd > y) *--zd = *--yd;
1347 7393844 : *--zd = x0[1];
1348 7393844 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1349 : }
1350 : GEN
1351 526600 : RgX_addmulXn(GEN x0, GEN y0, long d)
1352 : {
1353 : GEN x, y, xd, yd, zd;
1354 : long a, lz, nx, ny;
1355 :
1356 526600 : if (!signe(x0)) return RgX_copy(y0);
1357 525830 : nx = lgpol(x0);
1358 525830 : ny = lgpol(y0);
1359 525830 : zd = (GEN)avma;
1360 525830 : x = x0 + 2; y = y0 + 2; a = ny-d;
1361 525830 : if (a <= 0)
1362 : {
1363 291789 : lz = nx+d+2;
1364 291789 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1365 4294362 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1366 291789 : x = zd + a;
1367 757955 : while (zd > x) gel(--zd,0) = gen_0;
1368 : }
1369 : else
1370 : {
1371 234041 : xd = new_chunk(d); yd = y+d;
1372 234041 : x = RgX_addspec(x,yd, nx,a);
1373 234041 : lz = (a>nx)? ny+2: lg(x)+d;
1374 8625903 : x += 2; while (xd > x) *--zd = *--xd;
1375 : }
1376 2701771 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1377 525830 : *--zd = x0[1];
1378 525830 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1379 : }
1380 :
1381 : /* return x * y mod t^n */
1382 : static GEN
1383 6530589 : RgXn_mul_basecase(GEN x, GEN y, long n)
1384 : {
1385 6530589 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1386 : GEN z;
1387 6530589 : if (lx < 0) return pol_0(varn(x));
1388 6530589 : if (ly < 0) return pol_0(varn(x));
1389 6530589 : z = cgetg(lz, t_POL) + 2;
1390 6530589 : x+=2; if (lx > n) lx = n;
1391 6530589 : y+=2; if (ly > n) ly = n;
1392 6530589 : z[-1] = x[-1];
1393 6530589 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1394 6530589 : x = RgXspec_kill0(x, lx);
1395 6530589 : y = RgXspec_kill0(y, ly);
1396 : /* x:y:z [i] = term of degree i */
1397 25895550 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1398 11627419 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1399 6576171 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1400 6530589 : return normalizepol_lg(z - 2, lz);
1401 : }
1402 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1403 : static GEN
1404 9285812 : RgXn_mul2(GEN f, GEN g, long n)
1405 : {
1406 9285812 : pari_sp av = avma;
1407 : GEN fe,fo, ge,go, l,h,m;
1408 : long n0, n1;
1409 9285812 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1410 6566009 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1411 35420 : n0 = n>>1; n1 = n-n0;
1412 35420 : RgX_even_odd(f, &fe, &fo);
1413 35420 : RgX_even_odd(g, &ge, &go);
1414 35420 : l = RgXn_mul2(fe,ge,n1);
1415 35420 : h = RgXn_mul2(fo,go,n0);
1416 35420 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1417 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1418 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1419 35420 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1420 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1421 35420 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1422 35420 : m = RgX_inflate(m,2);
1423 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1424 35420 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1425 35420 : h = RgX_inflate(h,2);
1426 35420 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1427 35420 : return gerepileupto(av, h);
1428 : }
1429 : /* (f*g) \/ x^n */
1430 : static GEN
1431 1565360 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1432 : {
1433 1565360 : long d = degpol(f)+degpol(g) + 1 - n;
1434 : GEN h;
1435 1565360 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1436 29955 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1437 : RgX_recip_i(g), d));
1438 29955 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1439 : }
1440 :
1441 : /* (f*g) \/ x^n */
1442 : static GEN
1443 0 : RgX_sqrhigh_i2(GEN f, long n)
1444 : {
1445 0 : long d = 2*degpol(f)+ 1 - n;
1446 : GEN h;
1447 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1448 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1449 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1450 : }
1451 :
1452 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1453 : * b+2 were sent instead. na, nb = number of terms of a, b.
1454 : * Only c, c0, c1, c2 are genuine GEN.
1455 : */
1456 : GEN
1457 18406324 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1458 : {
1459 : GEN a0, c, c0;
1460 18406324 : long n0, n0a, i, v = 0;
1461 : pari_sp av;
1462 :
1463 24792959 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1464 24291050 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1465 18406280 : if (na < nb) swapspec(a,b, na,nb);
1466 18406280 : if (!nb) return pol_0(0);
1467 :
1468 18157156 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1469 491110 : RgX_shift_inplace_init(v);
1470 491141 : i = (na>>1); n0 = na-i; na = i;
1471 491141 : av = avma; a0 = a+n0; n0a = n0;
1472 1346657 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1473 :
1474 491141 : if (nb > n0)
1475 : {
1476 : GEN b0,c1,c2;
1477 : long n0b;
1478 :
1479 489763 : nb -= n0; b0 = b+n0; n0b = n0;
1480 1457742 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1481 489763 : c = RgX_mulspec(a,b,n0a,n0b);
1482 489763 : c0 = RgX_mulspec(a0,b0, na,nb);
1483 :
1484 489763 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1485 489763 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1486 :
1487 489763 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1488 489763 : c2 = RgX_sub(c1, RgX_add(c0,c));
1489 489763 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1490 : }
1491 : else
1492 : {
1493 1378 : c = RgX_mulspec(a,b,n0a,nb);
1494 1378 : c0 = RgX_mulspec(a0,b,na,nb);
1495 : }
1496 491141 : c0 = RgX_addmulXn(c0,c,n0);
1497 491141 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1498 : }
1499 :
1500 : INLINE GEN
1501 50715 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1502 : {
1503 50715 : pari_sp av = avma;
1504 50715 : GEN s = NULL;
1505 50715 : long j, l = (i+1)>>1;
1506 140745 : for (j=a; j<l; j++)
1507 : {
1508 90030 : GEN xj = gel(x,j), xx = gel(x,i-j);
1509 90030 : if (xj && xx)
1510 : {
1511 82386 : GEN t = gmul(xj, xx);
1512 82386 : s = s? gadd(s, t): t;
1513 : }
1514 : }
1515 50715 : if (s) s = gshift(s,1);
1516 50715 : if ((i&1) == 0)
1517 : {
1518 29748 : GEN t = gel(x, i>>1);
1519 29748 : if (t) {
1520 27354 : t = gsqr(t);
1521 27354 : s = s? gadd(s, t): t;
1522 : }
1523 : }
1524 50715 : return s? gerepileupto(av,s): gen_0;
1525 : }
1526 : static GEN
1527 8781 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1528 : {
1529 : long i, lz, nz;
1530 : GEN z;
1531 :
1532 8781 : if (!nx) return pol_0(0);
1533 8781 : x = RgXspec_kill0(x,nx);
1534 8781 : lz = (nx << 1) + 1, nz = lz-2;
1535 8781 : lz += v;
1536 8781 : z = cgetg(lz,t_POL) + 2;
1537 12953 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1538 38529 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1539 29748 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1540 8781 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1541 : }
1542 : /* return x^2 mod t^n */
1543 : static GEN
1544 0 : RgXn_sqr_basecase(GEN x, long n)
1545 : {
1546 0 : long i, lz = n+2, lx = lgpol(x);
1547 : GEN z;
1548 0 : if (lx < 0) return pol_0(varn(x));
1549 0 : z = cgetg(lz, t_POL);
1550 0 : z[1] = x[1];
1551 0 : x+=2; if (lx > n) lx = n;
1552 0 : x = RgXspec_kill0(x,lx);
1553 0 : z+=2;/* x:z [i] = term of degree i */
1554 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1555 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1556 0 : z -= 2; return normalizepol_lg(z, lz);
1557 : }
1558 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1559 : static GEN
1560 315 : RgXn_sqr2(GEN f, long n)
1561 : {
1562 315 : pari_sp av = avma;
1563 : GEN fe,fo, l,h,m;
1564 : long n0, n1;
1565 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1566 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1567 0 : n0 = n>>1; n1 = n-n0;
1568 0 : RgX_even_odd(f, &fe, &fo);
1569 0 : l = RgXn_sqr(fe,n1);
1570 0 : h = RgXn_sqr(fo,n0);
1571 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1572 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1573 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1574 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1575 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1576 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1577 0 : m = RgX_inflate(m,2);
1578 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1579 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1580 0 : h = RgX_inflate(h,2);
1581 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1582 0 : return gerepileupto(av, h);
1583 : }
1584 : GEN
1585 8820 : RgX_sqrspec(GEN a, long na)
1586 : {
1587 : GEN a0, c, c0, c1;
1588 8820 : long n0, n0a, i, v = 0;
1589 : pari_sp av;
1590 :
1591 10906 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1592 8820 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1593 39 : RgX_shift_inplace_init(v);
1594 39 : i = (na>>1); n0 = na-i; na = i;
1595 39 : av = avma; a0 = a+n0; n0a = n0;
1596 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1597 :
1598 39 : c = RgX_sqrspec(a,n0a);
1599 39 : c0 = RgX_sqrspec(a0,na);
1600 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1601 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1602 39 : c0 = RgX_addmulXn(c0,c,n0);
1603 39 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1604 : }
1605 :
1606 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1607 : GEN
1608 1751148 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1609 : {
1610 1751148 : GEN z = RgX_mul(A, B);
1611 1751131 : if (a < b)
1612 10426 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1613 1740705 : else if (a > b)
1614 1096711 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1615 : else
1616 643994 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1617 1751140 : return z;
1618 : }
1619 :
1620 : GEN
1621 16934314 : RgX_mul_i(GEN x, GEN y)
1622 : {
1623 16934314 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1624 16933747 : setvarn(z, varn(x)); return z;
1625 : }
1626 :
1627 : GEN
1628 8742 : RgX_sqr_i(GEN x)
1629 : {
1630 8742 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1631 8742 : setvarn(z,varn(x)); return z;
1632 : }
1633 :
1634 : /*******************************************************************/
1635 : /* */
1636 : /* DIVISION */
1637 : /* */
1638 : /*******************************************************************/
1639 : GEN
1640 3931705 : RgX_Rg_divexact(GEN x, GEN y) {
1641 3931705 : long i, lx = lg(x);
1642 : GEN z;
1643 3931705 : if (lx == 2) return gcopy(x);
1644 3882960 : switch(typ(y))
1645 : {
1646 3853181 : case t_INT:
1647 3853181 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1648 3683435 : break;
1649 4457 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1650 : }
1651 3708757 : z = cgetg(lx, t_POL); z[1] = x[1];
1652 23655081 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1653 3697137 : return z;
1654 : }
1655 : GEN
1656 28847005 : RgX_Rg_div(GEN x, GEN y) {
1657 28847005 : long i, lx = lg(x);
1658 : GEN z;
1659 28847005 : if (lx == 2) return gcopy(x);
1660 28586293 : switch(typ(y))
1661 : {
1662 19439993 : case t_INT:
1663 19439993 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1664 3743293 : break;
1665 4830 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1666 : }
1667 12884763 : z = cgetg(lx, t_POL); z[1] = x[1];
1668 48772653 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1669 12884394 : return normalizepol_lg(z, lx);
1670 : }
1671 : GEN
1672 13860 : RgX_normalize(GEN x)
1673 : {
1674 13860 : GEN z, d = NULL;
1675 13860 : long i, n = lg(x)-1;
1676 13860 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1677 13860 : if (i == 1) return pol_0(varn(x));
1678 13860 : if (i == n && isint1(d)) return x;
1679 5726 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1680 13783 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1681 5726 : gel(z,n) = Rg_get_1(d); return z;
1682 : }
1683 : GEN
1684 10353 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1685 : GEN
1686 230738 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1687 : {
1688 230738 : long l = lg(a), i;
1689 : GEN a0, z0, z;
1690 :
1691 230738 : if (l <= 3)
1692 : {
1693 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1694 0 : return pol_0(0);
1695 : }
1696 230738 : z = cgetg(l-1, t_POL);
1697 230743 : z[1] = a[1];
1698 230743 : a0 = a + l-1;
1699 230743 : z0 = z + l-2; *z0 = *a0--;
1700 2876905 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1701 : {
1702 2646167 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1703 2646162 : gel(z0,0) = t;
1704 : }
1705 230738 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1706 230738 : return z;
1707 : }
1708 : /* Polynomial division x / y:
1709 : * if pr = ONLY_REM return remainder, otherwise return quotient
1710 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1711 : * if pr != NULL set *pr to remainder, as the last object on stack */
1712 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1713 : static GEN
1714 18440172 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1715 : {
1716 : pari_sp avy, av, av1;
1717 : long dx,dy,dz,i,j,sx,lr;
1718 : GEN z,p1,p2,rem,y_lead,mod,p;
1719 : GEN (*f)(GEN,GEN);
1720 :
1721 18440172 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1722 :
1723 18440172 : dy = degpol(y);
1724 18440109 : y_lead = gel(y,dy+2);
1725 18440109 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1726 : {
1727 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1728 0 : for (dy--; dy>=0; dy--)
1729 : {
1730 0 : y_lead = gel(y,dy+2);
1731 0 : if (!gequal0(y_lead)) break;
1732 : }
1733 : }
1734 18440077 : if (!dy) /* y is constant */
1735 : {
1736 3140 : if (pr == ONLY_REM) return pol_0(varn(x));
1737 3133 : z = RgX_Rg_div(x, y_lead);
1738 3133 : if (pr == ONLY_DIVIDES) return z;
1739 1250 : if (pr) *pr = pol_0(varn(x));
1740 1250 : return z;
1741 : }
1742 18436937 : dx = degpol(x);
1743 18436881 : if (dx < dy)
1744 : {
1745 1937844 : if (pr == ONLY_REM) return RgX_copy(x);
1746 375086 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1747 375065 : z = pol_0(varn(x));
1748 375065 : if (pr) *pr = RgX_copy(x);
1749 375065 : return z;
1750 : }
1751 :
1752 : /* x,y in R[X], y non constant */
1753 16499037 : av = avma;
1754 16499037 : p = NULL;
1755 16499037 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1756 : {
1757 226359 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1758 226359 : if (!z) return gc_NULL(av);
1759 226359 : z = FpX_to_mod(z, p);
1760 226359 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1761 121821 : return gerepileupto(av, z);
1762 104538 : *pr = FpX_to_mod(*pr, p);
1763 104538 : return gc_all(av, 2, &z, pr);
1764 : }
1765 16272817 : switch(typ(y_lead))
1766 : {
1767 0 : case t_REAL:
1768 0 : y_lead = ginv(y_lead);
1769 0 : f = gmul; mod = NULL;
1770 0 : break;
1771 2143 : case t_INTMOD:
1772 2143 : case t_POLMOD: y_lead = ginv(y_lead);
1773 2143 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1774 2143 : break;
1775 16270674 : default: if (gequal1(y_lead)) y_lead = NULL;
1776 16270646 : f = gdiv; mod = NULL;
1777 : }
1778 :
1779 16272789 : if (y_lead == NULL)
1780 14449432 : p2 = gel(x,dx+2);
1781 : else {
1782 : for(;;) {
1783 1823357 : p2 = f(gel(x,dx+2),y_lead);
1784 1823363 : p2 = simplify_shallow(p2);
1785 1823363 : if (!isexactzero(p2) || (--dx < 0)) break;
1786 : }
1787 1823363 : if (dx < dy) /* leading coeff of x was in fact zero */
1788 : {
1789 0 : if (pr == ONLY_DIVIDES) {
1790 0 : set_avma(av);
1791 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1792 : }
1793 0 : if (pr == ONLY_REM)
1794 : {
1795 0 : if (dx < 0)
1796 0 : return gerepilecopy(av, scalarpol(p2, varn(x)));
1797 : else
1798 : {
1799 : GEN t;
1800 0 : set_avma(av);
1801 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1802 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1803 0 : return t;
1804 : }
1805 : }
1806 0 : if (pr) /* cf ONLY_REM above */
1807 : {
1808 0 : if (dx < 0)
1809 : {
1810 0 : p2 = gclone(p2);
1811 0 : set_avma(av);
1812 0 : z = pol_0(varn(x));
1813 0 : x = scalarpol(p2, varn(x));
1814 0 : gunclone(p2);
1815 : }
1816 : else
1817 : {
1818 : GEN t;
1819 0 : set_avma(av);
1820 0 : z = pol_0(varn(x));
1821 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1822 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1823 0 : x = t;
1824 : }
1825 0 : *pr = x;
1826 : }
1827 : else
1828 : {
1829 0 : set_avma(av);
1830 0 : z = pol_0(varn(x));
1831 : }
1832 0 : return z;
1833 : }
1834 : }
1835 : /* dx >= dy */
1836 16272795 : avy = avma;
1837 16272795 : dz = dx-dy;
1838 16272795 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1839 16272770 : x += 2;
1840 16272770 : z += 2;
1841 16272770 : y += 2;
1842 16272770 : gel(z,dz) = gcopy(p2);
1843 :
1844 46312396 : for (i=dx-1; i>=dy; i--)
1845 : {
1846 30039660 : av1=avma; p1=gel(x,i);
1847 1100367342 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1848 30033409 : if (y_lead) p1 = simplify(f(p1,y_lead));
1849 :
1850 30033409 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1851 : else
1852 20156816 : p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1853 30039027 : gel(z,i-dy) = p1;
1854 : }
1855 16272736 : if (!pr) return gerepileupto(av,z-2);
1856 :
1857 9478736 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1858 9839498 : for (sx=0; ; i--)
1859 : {
1860 9839498 : p1 = gel(x,i);
1861 : /* we always enter this loop at least once */
1862 26341713 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1863 9838002 : if (mod && avma==av1) p1 = gmul(p1,mod);
1864 9838002 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1865 3503475 : if (!isexactzero(p1)) break;
1866 3425063 : if (!i) break;
1867 360829 : set_avma(av1);
1868 : }
1869 9477777 : if (pr == ONLY_DIVIDES)
1870 : {
1871 1911 : if (sx) return gc_NULL(av);
1872 1890 : set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1873 : }
1874 9475866 : lr=i+3; rem -= lr;
1875 9475866 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1876 9432247 : else p1 = gerepileupto((pari_sp)rem,p1);
1877 9476751 : rem[0] = evaltyp(t_POL) | evallg(lr);
1878 9476726 : rem[1] = z[-1];
1879 9476726 : rem += 2;
1880 9476726 : gel(rem,i) = p1;
1881 14263374 : for (i--; i>=0; i--)
1882 : {
1883 4786680 : av1=avma; p1 = gel(x,i);
1884 13571954 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1885 4786029 : if (mod && avma==av1) p1 = gmul(p1,mod);
1886 4786312 : gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1887 : }
1888 9476694 : rem -= 2;
1889 9476694 : if (!sx) (void)normalizepol_lg(rem, lr);
1890 9476892 : if (pr == ONLY_REM) return gerepileupto(av,rem);
1891 5495263 : z -= 2;
1892 : {
1893 5495263 : GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1894 5495263 : gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1895 : }
1896 : }
1897 :
1898 : GEN
1899 12895783 : RgX_divrem(GEN x, GEN y, GEN *pr)
1900 : {
1901 12895783 : if (pr == ONLY_REM) return RgX_rem(x, y);
1902 12895783 : return RgX_divrem_i(x, y, pr);
1903 : }
1904 :
1905 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1906 : GEN
1907 134919 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1908 : {
1909 : long vx, dx, dy, dz, i, j, sx, lr;
1910 : pari_sp av0, av, tetpil;
1911 : GEN z,p1,rem,lead;
1912 :
1913 134919 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1914 134919 : vx = varn(x);
1915 134919 : dx = degpol(x);
1916 134919 : dy = degpol(y);
1917 134919 : if (dx < dy)
1918 : {
1919 27219 : if (pr)
1920 : {
1921 27219 : av0 = avma; x = RgXQX_red(x, T);
1922 27219 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1923 27219 : if (pr == ONLY_REM) return x;
1924 0 : *pr = x;
1925 : }
1926 0 : return pol_0(vx);
1927 : }
1928 107700 : lead = leading_coeff(y);
1929 107700 : if (!dy) /* y is constant */
1930 : {
1931 546 : if (pr && pr != ONLY_DIVIDES)
1932 : {
1933 0 : if (pr == ONLY_REM) return pol_0(vx);
1934 0 : *pr = pol_0(vx);
1935 : }
1936 546 : if (gequal1(lead)) return RgX_copy(x);
1937 0 : av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
1938 0 : return gerepile(av0,tetpil,RgXQX_red(x,T));
1939 : }
1940 107154 : av0 = avma; dz = dx-dy;
1941 107154 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
1942 107154 : set_avma(av0);
1943 107154 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1944 107154 : x += 2; y += 2; z += 2;
1945 :
1946 107154 : p1 = gel(x,dx); av = avma;
1947 107154 : gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
1948 550100 : for (i=dx-1; i>=dy; i--)
1949 : {
1950 442946 : av=avma; p1=gel(x,i);
1951 2335106 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1952 442922 : if (lead) p1 = gmul(grem(p1, T), lead);
1953 442929 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
1954 : }
1955 107154 : if (!pr) { guncloneNULL(lead); return z-2; }
1956 :
1957 105908 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
1958 178624 : for (sx=0; ; i--)
1959 : {
1960 178624 : p1 = gel(x,i);
1961 615735 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1962 178626 : tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
1963 99784 : if (!i) break;
1964 72716 : set_avma(av);
1965 : }
1966 105907 : if (pr == ONLY_DIVIDES)
1967 : {
1968 27026 : guncloneNULL(lead);
1969 27026 : if (sx) return gc_NULL(av0);
1970 24520 : return gc_const((pari_sp)rem, z-2);
1971 : }
1972 78881 : lr=i+3; rem -= lr;
1973 78881 : rem[0] = evaltyp(t_POL) | evallg(lr);
1974 78881 : rem[1] = z[-1];
1975 78881 : p1 = gerepile((pari_sp)rem,tetpil,p1);
1976 78881 : rem += 2; gel(rem,i) = p1;
1977 169101 : for (i--; i>=0; i--)
1978 : {
1979 90220 : av=avma; p1 = gel(x,i);
1980 305921 : for (j=0; j<=i && j<=dz; j++)
1981 215701 : p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1982 90220 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
1983 : }
1984 78881 : rem -= 2;
1985 78881 : guncloneNULL(lead);
1986 78881 : if (!sx) (void)normalizepol_lg(rem, lr);
1987 78881 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
1988 161 : *pr = rem; return z-2;
1989 : }
1990 :
1991 : /*******************************************************************/
1992 : /* */
1993 : /* PSEUDO-DIVISION */
1994 : /* */
1995 : /*******************************************************************/
1996 : INLINE GEN
1997 1314966 : rem(GEN c, GEN T)
1998 : {
1999 1314966 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2000 1314968 : return c;
2001 : }
2002 :
2003 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2004 : int
2005 12599 : ZXQX_dvd(GEN x, GEN y, GEN T)
2006 : {
2007 : long dx, dy, i, T_ismonic;
2008 12599 : pari_sp av = avma, av2;
2009 : GEN y_lead;
2010 :
2011 12599 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2012 12599 : dy = degpol(y); y_lead = gel(y,dy+2);
2013 12599 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2014 : /* if monic, no point in using pseudo-division */
2015 12599 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2016 10086 : T_ismonic = gequal1(leading_coeff(T));
2017 10086 : dx = degpol(x);
2018 10086 : if (dx < dy) return !signe(x);
2019 10086 : (void)new_chunk(2);
2020 10086 : x = RgX_recip_i(x)+2;
2021 10086 : y = RgX_recip_i(y)+2;
2022 : /* pay attention to sparse divisors */
2023 20668 : for (i = 1; i <= dy; i++)
2024 10582 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2025 10086 : av2 = avma;
2026 : for (;;)
2027 59662 : {
2028 69748 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2029 69748 : x0 = gneg(x0);
2030 69748 : m = gcdii(cx, y0);
2031 69748 : if (!equali1(m))
2032 : {
2033 67868 : x0 = gdiv(x0, m);
2034 67868 : y0 = diviiexact(y0, m);
2035 67866 : if (equali1(y0)) y0 = NULL;
2036 : }
2037 143207 : for (i=1; i<=dy; i++)
2038 : {
2039 73459 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2040 73459 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2041 73461 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2042 73461 : gel(x,i) = c;
2043 : }
2044 626199 : for ( ; i<=dx; i++)
2045 : {
2046 556451 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2047 556451 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2048 556451 : gel(x,i) = c;
2049 : }
2050 79996 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2051 69748 : if (dx < dy) break;
2052 59662 : if (gc_needed(av2,1))
2053 : {
2054 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2055 28 : gerepilecoeffs(av2,x,dx+1);
2056 : }
2057 : }
2058 10086 : return gc_bool(av, dx < 0);
2059 : }
2060 :
2061 : /* T either NULL or a t_POL. */
2062 : GEN
2063 134570 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2064 : {
2065 134570 : long vx = varn(x), dx, dy, dz, i, lx, p;
2066 134570 : pari_sp av = avma, av2;
2067 : GEN y_lead;
2068 :
2069 134570 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2070 134570 : dy = degpol(y); y_lead = gel(y,dy+2);
2071 : /* if monic, no point in using pseudo-division */
2072 134570 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2073 128500 : dx = degpol(x);
2074 128500 : if (dx < dy) return RgX_copy(x);
2075 128500 : (void)new_chunk(2);
2076 128499 : x = RgX_recip_i(x)+2;
2077 128500 : y = RgX_recip_i(y)+2;
2078 : /* pay attention to sparse divisors */
2079 366867 : for (i = 1; i <= dy; i++)
2080 238366 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2081 128501 : dz = dx-dy; p = dz+1;
2082 128501 : av2 = avma;
2083 : for (;;)
2084 : {
2085 357302 : gel(x,0) = gneg(gel(x,0)); p--;
2086 1059708 : for (i=1; i<=dy; i++)
2087 : {
2088 702414 : GEN c = gmul(y_lead, gel(x,i));
2089 702374 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2090 702401 : gel(x,i) = rem(c, T);
2091 : }
2092 848493 : for ( ; i<=dx; i++)
2093 : {
2094 491194 : GEN c = gmul(y_lead, gel(x,i));
2095 491180 : gel(x,i) = rem(c, T);
2096 : }
2097 383783 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2098 357300 : if (dx < dy) break;
2099 228801 : if (gc_needed(av2,1))
2100 : {
2101 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2102 0 : gerepilecoeffs(av2,x,dx+1);
2103 : }
2104 : }
2105 128499 : if (dx < 0) return pol_0(vx);
2106 127316 : lx = dx+3; x -= 2;
2107 127316 : x[0] = evaltyp(t_POL) | evallg(lx);
2108 127316 : x[1] = evalsigne(1) | evalvarn(vx);
2109 127316 : x = RgX_recip_i(x);
2110 127316 : if (p)
2111 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2112 12051 : GEN t = y_lead;
2113 12051 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2114 0 : t = RgXQ_powu(t, p, T);
2115 : else
2116 12051 : t = gpowgs(t, p);
2117 32276 : for (i=2; i<lx; i++)
2118 : {
2119 20226 : GEN c = gmul(gel(x,i), t);
2120 20225 : gel(x,i) = rem(c,T);
2121 : }
2122 12050 : if (!T) return gerepileupto(av, x);
2123 : }
2124 115265 : return gerepilecopy(av, x);
2125 : }
2126 :
2127 : GEN
2128 134570 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2129 :
2130 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2131 : GEN
2132 11955 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2133 : {
2134 11955 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2135 11955 : pari_sp av = avma, av2;
2136 : GEN z, r, ypow, y_lead;
2137 :
2138 11955 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2139 11955 : dy = degpol(y); y_lead = gel(y,dy+2);
2140 11955 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2141 7587 : dx = degpol(x);
2142 7587 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2143 7587 : if (dx == dy)
2144 : {
2145 98 : GEN x_lead = gel(x,lg(x)-1);
2146 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2147 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2148 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2149 98 : *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
2150 : }
2151 7489 : (void)new_chunk(2);
2152 7489 : x = RgX_recip_i(x)+2;
2153 7489 : y = RgX_recip_i(y)+2;
2154 : /* pay attention to sparse divisors */
2155 38804 : for (i = 1; i <= dy; i++)
2156 31315 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2157 7489 : dz = dx-dy; p = dz+1;
2158 7489 : lz = dz+3;
2159 7489 : z = cgetg(lz, t_POL);
2160 7489 : z[1] = evalsigne(1) | evalvarn(vx);
2161 28720 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2162 7489 : ypow = new_chunk(dz+1);
2163 7489 : gel(ypow,0) = gen_1;
2164 7489 : gel(ypow,1) = y_lead;
2165 13742 : for (i=2; i<=dz; i++)
2166 : {
2167 6253 : GEN c = gmul(gel(ypow,i-1), y_lead);
2168 6253 : gel(ypow,i) = rem(c,T);
2169 : }
2170 7489 : av2 = avma;
2171 7489 : for (iz=2;;)
2172 : {
2173 15657 : p--;
2174 15657 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2175 69483 : for (i=1; i<=dy; i++)
2176 : {
2177 53826 : GEN c = gmul(y_lead, gel(x,i));
2178 53826 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2179 53826 : gel(x,i) = rem(c, T);
2180 : }
2181 41117 : for ( ; i<=dx; i++)
2182 : {
2183 25460 : GEN c = gmul(y_lead, gel(x,i));
2184 25460 : gel(x,i) = rem(c,T);
2185 : }
2186 15657 : x++; dx--;
2187 21231 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2188 15657 : if (dx < dy) break;
2189 8168 : if (gc_needed(av2,1))
2190 : {
2191 0 : GEN X = x-2;
2192 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2193 0 : X[0] = evaltyp(t_POL)|evallg(dx+3); X[1] = z[1]; /* hack */
2194 0 : gerepileall(av2,2, &X, &z); x = X+2;
2195 : }
2196 : }
2197 13991 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2198 7489 : if (dx < 0)
2199 182 : x = pol_0(vx);
2200 : else
2201 : {
2202 7307 : lx = dx+3; x -= 2;
2203 7307 : x[0] = evaltyp(t_POL) | evallg(lx);
2204 7307 : x[1] = evalsigne(1) | evalvarn(vx);
2205 7307 : x = RgX_recip_i(x);
2206 : }
2207 7489 : z = RgX_recip_i(z);
2208 7489 : r = x;
2209 7489 : if (p)
2210 : {
2211 3325 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2212 3325 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2213 : }
2214 7489 : *ptr = r; return gc_all(av, 2, &z, ptr);
2215 : }
2216 : GEN
2217 11717 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2218 11717 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2219 :
2220 : GEN
2221 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2222 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2223 : GEN
2224 485345935 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(gmul(y, gel(x,i))); }
2225 : GEN
2226 141344 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2227 : GEN
2228 25963 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2229 : GEN
2230 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2231 : GEN
2232 133 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2233 :
2234 : GEN
2235 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2236 :
2237 : GEN
2238 0 : RgXQX_powers(GEN P, long n, GEN T)
2239 : {
2240 0 : GEN v = cgetg(n+2, t_VEC);
2241 : long i;
2242 0 : gel(v, 1) = pol_1(varn(T));
2243 0 : if (n==0) return v;
2244 0 : gel(v, 2) = gcopy(P);
2245 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2246 0 : return v;
2247 : }
2248 :
2249 : static GEN
2250 600415 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2251 : static GEN
2252 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2253 : static GEN
2254 71958 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2255 : static GEN
2256 87215 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2257 : static GEN
2258 269831 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2259 : static GEN
2260 1000311 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2261 : static GEN
2262 975282 : _one(void *data) { return pol_1(varn((GEN)data)); }
2263 : static GEN
2264 826 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2265 : static GEN
2266 700648 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2267 :
2268 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2269 : _mul, _sqr, _one, _zero };
2270 :
2271 : GEN
2272 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2273 : {
2274 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2275 : }
2276 :
2277 : GEN
2278 399747 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2279 : {
2280 399747 : int use_sqr = 2*degpol(x) >= degpol(T);
2281 399748 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2282 : }
2283 :
2284 : /* mod X^n */
2285 : struct modXn {
2286 : long v; /* varn(X) */
2287 : long n;
2288 : } ;
2289 : static GEN
2290 11893 : _sqrXn(void *data, GEN x) {
2291 11893 : struct modXn *S = (struct modXn*)data;
2292 11893 : return RgXn_sqr(x, S->n);
2293 : }
2294 : static GEN
2295 4528 : _mulXn(void *data, GEN x, GEN y) {
2296 4528 : struct modXn *S = (struct modXn*)data;
2297 4528 : return RgXn_mul(x,y, S->n);
2298 : }
2299 : static GEN
2300 1939 : _oneXn(void *data) {
2301 1939 : struct modXn *S = (struct modXn*)data;
2302 1939 : return pol_1(S->v);
2303 : }
2304 : static GEN
2305 0 : _zeroXn(void *data) {
2306 0 : struct modXn *S = (struct modXn*)data;
2307 0 : return pol_0(S->v);
2308 : }
2309 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2310 : _oneXn, _zeroXn };
2311 :
2312 : GEN
2313 483 : RgXn_powers(GEN x, long m, long n)
2314 : {
2315 483 : long d = degpol(x);
2316 483 : int use_sqr = (d<<1) >= n;
2317 : struct modXn S;
2318 483 : S.v = varn(x); S.n = n;
2319 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2320 : }
2321 :
2322 : GEN
2323 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2324 : {
2325 : struct modXn S;
2326 : long v;
2327 2286 : if (n == 0) return x;
2328 2286 : v = RgX_valrem(x, &x);
2329 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2330 2265 : S.v = varn(x); S.n = n;
2331 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2332 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2333 2265 : return x;
2334 : }
2335 : GEN
2336 0 : RgXn_powu(GEN x, ulong m, long n)
2337 : {
2338 : pari_sp av;
2339 0 : if (n == 0) return gcopy(x);
2340 0 : av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
2341 : }
2342 :
2343 : GEN
2344 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2345 : {
2346 : struct modXn S;
2347 966 : S.v = varn(gel(x,2)); S.n = n;
2348 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2349 : }
2350 :
2351 : GEN
2352 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2353 : {
2354 0 : int use_sqr = 2*degpol(x) >= n;
2355 : struct modXn S;
2356 0 : S.v = varn(x); S.n = n;
2357 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2358 : }
2359 :
2360 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2361 : GEN
2362 4347 : RgXn_eval(GEN Q, GEN x, long n)
2363 : {
2364 4347 : long d = degpol(x);
2365 : int use_sqr;
2366 : struct modXn S;
2367 4347 : if (d == 1 && isrationalzero(gel(x,2)))
2368 : {
2369 4340 : GEN y = RgX_unscale(Q, gel(x,3));
2370 4340 : setvarn(y, varn(x)); return y;
2371 : }
2372 7 : S.v = varn(x);
2373 7 : S.n = n;
2374 7 : use_sqr = (d<<1) >= n;
2375 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2376 : }
2377 :
2378 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2379 : static GEN
2380 2567731 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2381 : {
2382 2567731 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2383 2567731 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2384 : }
2385 :
2386 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2387 : static GEN
2388 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2389 : {
2390 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2391 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2392 : }
2393 :
2394 : static GEN
2395 2166423 : RgXn_div_gen(GEN g, GEN f, long e)
2396 : {
2397 : pari_sp av;
2398 : ulong mask;
2399 : GEN W, a;
2400 2166423 : long v = varn(f), n = 1;
2401 :
2402 2166423 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2403 2166416 : a = ginv(gel(f,2));
2404 2166413 : if (e == 1 && !g) return scalarpol(a, v);
2405 2162388 : else if (e == 2 && !g)
2406 : {
2407 : GEN b;
2408 835129 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2409 244636 : b = gneg(b);
2410 244638 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2411 244638 : return deg1pol(b, a, v);
2412 : }
2413 1327259 : av = avma;
2414 1327259 : W = scalarpol_shallow(a,v);
2415 1327258 : mask = quadratic_prec_mask(e);
2416 3889123 : while (mask > 1)
2417 : {
2418 : GEN u, fr;
2419 2561865 : long n2 = n;
2420 2561865 : n<<=1; if (mask & 1) n--;
2421 2561865 : mask >>= 1;
2422 2561865 : fr = RgXn_red_shallow(f, n);
2423 2561865 : if (mask>1 || !g)
2424 : {
2425 1292837 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2426 1292837 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2427 : }
2428 : else
2429 : {
2430 1269028 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2431 1269028 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2432 1269028 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2433 : }
2434 2561865 : if (gc_needed(av,2))
2435 : {
2436 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2437 0 : W = gerepileupto(av, W);
2438 : }
2439 : }
2440 1327258 : return W;
2441 : }
2442 :
2443 : static GEN
2444 119 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2445 : {
2446 : GEN r;
2447 119 : if (lgefint(p) == 3)
2448 : {
2449 119 : ulong pp = uel(p, 2);
2450 119 : if (pp == 2)
2451 14 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2452 : else
2453 105 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2454 : }
2455 : else
2456 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2457 119 : return FpX_to_mod(r, p);
2458 : }
2459 :
2460 : static GEN
2461 0 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2462 : {
2463 0 : GEN r, T = RgX_to_FpX(pol, p);
2464 0 : if (signe(T) == 0) pari_err_OP("/", x, y);
2465 0 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2466 0 : return FpXQX_to_mod(r, T, p);
2467 : }
2468 :
2469 : static GEN
2470 70 : RgXn_inv_FpX(GEN x, long e, GEN p)
2471 : {
2472 : GEN r;
2473 70 : if (lgefint(p) == 3)
2474 : {
2475 70 : ulong pp = uel(p, 2);
2476 70 : if (pp == 2)
2477 14 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2478 : else
2479 56 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2480 : }
2481 : else
2482 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2483 70 : return FpX_to_mod(r, p);
2484 : }
2485 :
2486 : static GEN
2487 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2488 : {
2489 0 : GEN r, T = RgX_to_FpX(pol, p);
2490 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2491 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2492 0 : return FpXQX_to_mod(r, T, p);
2493 : }
2494 :
2495 : #define code(t1,t2) ((t1 << 6) | t2)
2496 :
2497 : static GEN
2498 897462 : RgXn_inv_fast(GEN x, long e)
2499 : {
2500 : GEN p, pol;
2501 : long pa;
2502 897462 : long t = RgX_type(x,&p,&pol,&pa);
2503 897466 : switch(t)
2504 : {
2505 70 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2506 0 : case code(t_POLMOD, t_INTMOD):
2507 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2508 897396 : default: return NULL;
2509 : }
2510 : }
2511 :
2512 : static GEN
2513 1269147 : RgXn_div_fast(GEN x, GEN y, long e)
2514 : {
2515 : GEN p, pol;
2516 : long pa;
2517 1269147 : long t = RgX_type2(x,y,&p,&pol,&pa);
2518 1269147 : switch(t)
2519 : {
2520 119 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2521 0 : case code(t_POLMOD, t_INTMOD):
2522 0 : return RgXn_div_FpXQX(x, y, e, pol, p);
2523 1269028 : default: return NULL;
2524 : }
2525 : }
2526 : #undef code
2527 :
2528 : GEN
2529 1269147 : RgXn_div_i(GEN g, GEN f, long e)
2530 : {
2531 1269147 : GEN h = RgXn_div_fast(g, f, e);
2532 1269147 : if (h) return h;
2533 1269028 : return RgXn_div_gen(g, f, e);
2534 : }
2535 :
2536 : GEN
2537 578248 : RgXn_div(GEN g, GEN f, long e)
2538 : {
2539 578248 : pari_sp av = avma;
2540 578248 : return gerepileupto(av, RgXn_div_i(g, f, e));
2541 : }
2542 :
2543 : GEN
2544 897462 : RgXn_inv_i(GEN f, long e)
2545 : {
2546 897462 : GEN h = RgXn_inv_fast(f, e);
2547 897465 : if (h) return h;
2548 897395 : return RgXn_div_gen(NULL, f, e);
2549 : }
2550 :
2551 : GEN
2552 806972 : RgXn_inv(GEN f, long e)
2553 : {
2554 806972 : pari_sp av = avma;
2555 806972 : return gerepileupto(av, RgXn_inv_i(f, e));
2556 : }
2557 :
2558 : /* Compute intformal(x^n*S)/x^(n+1) */
2559 : static GEN
2560 58711 : RgX_integXn(GEN x, long n)
2561 : {
2562 58711 : long i, lx = lg(x);
2563 : GEN y;
2564 58711 : if (lx == 2) return RgX_copy(x);
2565 52677 : y = cgetg(lx, t_POL); y[1] = x[1];
2566 110473 : for (i=2; i<lx; i++)
2567 57794 : gel(y,i) = gdivgs(gel(x,i), n+i-1);
2568 52679 : return RgX_renormalize_lg(y, lx);;
2569 : }
2570 :
2571 : GEN
2572 52856 : RgXn_expint(GEN h, long e)
2573 : {
2574 52856 : pari_sp av = avma, av2;
2575 52856 : long v = varn(h), n;
2576 52856 : GEN f = pol_1(v), g;
2577 : ulong mask;
2578 :
2579 52856 : if (!signe(h)) return f;
2580 52856 : g = pol_1(v);
2581 52856 : n = 1; mask = quadratic_prec_mask(e);
2582 52856 : av2 = avma;
2583 58715 : for (;mask>1;)
2584 : {
2585 : GEN u, w;
2586 58715 : long n2 = n;
2587 58715 : n<<=1; if (mask & 1) n--;
2588 58715 : mask >>= 1;
2589 58715 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2590 58713 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2591 58711 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2592 58713 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2593 58713 : if (mask<=1) break;
2594 5859 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2595 5859 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2596 5859 : if (gc_needed(av2,2))
2597 : {
2598 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2599 0 : gerepileall(av2, 2, &f, &g);
2600 : }
2601 : }
2602 52854 : return gerepileupto(av, f);
2603 : }
2604 :
2605 : GEN
2606 0 : RgXn_exp(GEN h, long e)
2607 : {
2608 0 : long d = degpol(h);
2609 0 : if (d < 0) return pol_1(varn(h));
2610 0 : if (!d || !gequal0(gel(h,2)))
2611 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2612 0 : return RgXn_expint(RgX_deriv(h), e);
2613 : }
2614 :
2615 : GEN
2616 154 : RgXn_reverse(GEN f, long e)
2617 : {
2618 154 : pari_sp av = avma, av2;
2619 : ulong mask;
2620 : GEN fi, a, df, W, an;
2621 154 : long v = varn(f), n=1;
2622 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2623 0 : pari_err_INV("serreverse",f);
2624 154 : fi = ginv(gel(f,3));
2625 154 : a = deg1pol_shallow(fi,gen_0,v);
2626 154 : if (e <= 2) return gerepilecopy(av, a);
2627 133 : W = scalarpol(fi,v);
2628 133 : df = RgX_deriv(f);
2629 133 : mask = quadratic_prec_mask(e);
2630 133 : av2 = avma;
2631 616 : for (;mask>1;)
2632 : {
2633 : GEN u, fa, fr;
2634 483 : long n2 = n, rt;
2635 483 : n<<=1; if (mask & 1) n--;
2636 483 : mask >>= 1;
2637 483 : fr = RgXn_red_shallow(f, n);
2638 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2639 483 : an = RgXn_powers(a, rt, n);
2640 483 : if (n>1)
2641 : {
2642 483 : long n4 = (n2+1)>>1;
2643 483 : GEN dfr = RgXn_red_shallow(df, n2);
2644 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2645 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2646 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2647 : }
2648 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2649 483 : fa = RgX_shift(fa, -n2);
2650 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2651 483 : if (gc_needed(av2,2))
2652 : {
2653 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2654 0 : gerepileall(av2, 2, &a, &W);
2655 : }
2656 : }
2657 133 : return gerepileupto(av, a);
2658 : }
2659 :
2660 : GEN
2661 7 : RgXn_sqrt(GEN h, long e)
2662 : {
2663 7 : pari_sp av = avma, av2;
2664 7 : long v = varn(h), n = 1;
2665 7 : GEN f = scalarpol(gen_1, v), df = f;
2666 7 : ulong mask = quadratic_prec_mask(e);
2667 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2668 0 : pari_err_SQRTN("RgXn_sqrt",h);
2669 7 : av2 = avma;
2670 : while(1)
2671 7 : {
2672 14 : long n2 = n, m;
2673 : GEN g;
2674 14 : n<<=1; if (mask & 1) n--;
2675 14 : mask >>= 1;
2676 14 : m = n-n2;
2677 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2678 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2679 14 : if (mask==1) return gerepileupto(av, f);
2680 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2681 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2682 7 : if (gc_needed(av2,2))
2683 : {
2684 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2685 0 : gerepileall(av2, 2, &f, &df);
2686 : }
2687 : }
2688 : }
2689 :
2690 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2691 : GEN
2692 114059 : RgXQ_powu(GEN x, ulong n, GEN T)
2693 : {
2694 114059 : pari_sp av = avma;
2695 :
2696 114059 : if (!n) return pol_1(varn(x));
2697 66367 : if (n == 1) return RgX_copy(x);
2698 18357 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2699 18357 : return gerepilecopy(av, x);
2700 : }
2701 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2702 : GEN
2703 109424 : RgXQ_pow(GEN x, GEN n, GEN T)
2704 : {
2705 : pari_sp av;
2706 109424 : long s = signe(n);
2707 :
2708 109424 : if (!s) return pol_1(varn(x));
2709 109424 : if (is_pm1(n) == 1)
2710 87215 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2711 22209 : av = avma;
2712 22209 : if (s < 0) x = RgXQ_inv(x, T);
2713 22209 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2714 22209 : return gerepilecopy(av, x);
2715 : }
2716 : static GEN
2717 213319 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2718 : static GEN
2719 112816 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2720 :
2721 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2722 : GEN
2723 10658 : ZXQ_powers(GEN x, long l, GEN T)
2724 : {
2725 10658 : int use_sqr = 2*degpol(x) >= degpol(T);
2726 10658 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2727 : }
2728 :
2729 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2730 : GEN
2731 188501 : ZXQ_powu(GEN x, ulong n, GEN T)
2732 : {
2733 188501 : pari_sp av = avma;
2734 :
2735 188501 : if (!n) return pol_1(varn(x));
2736 188501 : if (n == 1) return ZX_copy(x);
2737 125025 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2738 125019 : return gerepilecopy(av, x);
2739 : }
2740 :
2741 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2742 : GEN
2743 4207 : RgXQ_powers(GEN x, long l, GEN T)
2744 : {
2745 4207 : int use_sqr = 2*degpol(x) >= degpol(T);
2746 4207 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2747 : }
2748 :
2749 : GEN
2750 7035 : RgXQV_factorback(GEN L, GEN e, GEN T)
2751 : {
2752 7035 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2753 : }
2754 :
2755 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2756 : GEN
2757 8096 : QXQ_powers(GEN a, long n, GEN T)
2758 : {
2759 : GEN den, v;
2760 8096 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2761 8075 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2762 : /* den*a integral; v[i+1] = (den*a)^i in K */
2763 8075 : if (den)
2764 : { /* restore denominators */
2765 5111 : GEN d = den;
2766 : long i;
2767 5111 : gel(v,2) = a;
2768 25721 : for (i=3; i<=n+1; i++) {
2769 20610 : d = mulii(d,den);
2770 20610 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2771 : }
2772 : }
2773 8075 : return v;
2774 : }
2775 :
2776 : static GEN
2777 2387 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2778 : {
2779 2387 : long l, i, m = 0;
2780 : GEN dz, z;
2781 2387 : GEN V = cgetg_copy(v, &l);
2782 8379 : for (i = imin; i < l; i++)
2783 : {
2784 5992 : GEN c = gel(v, i);
2785 5992 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2786 : }
2787 2387 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2788 2618 : for (i = 1; i < imin; i++) V[i] = v[i];
2789 8379 : for (i = imin; i < l; i++)
2790 : {
2791 5992 : GEN c = gel(v,i);
2792 5992 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2793 5992 : gel(V,i) = c;
2794 : }
2795 2387 : return V;
2796 : }
2797 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2798 : GEN
2799 2156 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2800 2156 : { return do_QXQ_eval(v, 1, a, T); }
2801 :
2802 : GEN
2803 231 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2804 231 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2805 :
2806 : GEN
2807 2891 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2808 : {
2809 2891 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2810 : }
2811 :
2812 : GEN
2813 5096 : RgXQ_norm(GEN x, GEN T)
2814 : {
2815 : pari_sp av;
2816 5096 : long dx = degpol(x);
2817 : GEN L, y;
2818 5096 : if (degpol(T)==0) return gpowgs(x,0);
2819 5089 : av = avma; y = resultant(T, x);
2820 5089 : L = leading_coeff(T);
2821 5089 : if (gequal1(L) || !signe(x)) return y;
2822 0 : return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2823 : }
2824 :
2825 : GEN
2826 441 : RgXQ_trace(GEN x, GEN T)
2827 : {
2828 441 : pari_sp av = avma;
2829 : GEN dT, z;
2830 : long n;
2831 441 : if (degpol(T)==0) return gmulgs(x,0);
2832 434 : dT = RgX_deriv(T); n = degpol(dT);
2833 434 : z = RgXQ_mul(x, dT, T);
2834 434 : if (degpol(z)<n) return gc_const(av, gen_0);
2835 385 : return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2836 : }
2837 :
2838 : GEN
2839 2741768 : RgX_blocks(GEN P, long n, long m)
2840 : {
2841 2741768 : GEN z = cgetg(m+1,t_VEC);
2842 2741768 : long i,j, k=2, l = lg(P);
2843 8418882 : for(i=1; i<=m; i++)
2844 : {
2845 5677114 : GEN zi = cgetg(n+2,t_POL);
2846 5677114 : zi[1] = P[1];
2847 5677114 : gel(z,i) = zi;
2848 16744532 : for(j=2; j<n+2; j++)
2849 11067418 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2850 5677114 : zi = RgX_renormalize_lg(zi, n+2);
2851 : }
2852 2741768 : return z;
2853 : }
2854 :
2855 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2856 : void
2857 29597336 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2858 : {
2859 29597336 : long n = degpol(p), v = varn(p), n0, n1, i;
2860 : GEN p0, p1;
2861 :
2862 29596688 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2863 :
2864 29503759 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2865 29503759 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2866 29507655 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2867 95667304 : for (i=0; i<n1; i++)
2868 : {
2869 66156297 : p0[2+i] = p[2+(i<<1)];
2870 66156297 : p1[2+i] = p[3+(i<<1)];
2871 : }
2872 29511007 : if (n1 != n0)
2873 13899825 : p0[2+i] = p[2+(i<<1)];
2874 29511007 : *pe = normalizepol(p0);
2875 29502790 : *po = normalizepol(p1);
2876 : }
2877 :
2878 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2879 : GEN
2880 43554 : RgX_splitting(GEN p, long k)
2881 : {
2882 43554 : long n = degpol(p), v = varn(p), m, i, j, l;
2883 : GEN r;
2884 :
2885 43554 : m = n/k;
2886 43554 : r = cgetg(k+1,t_VEC);
2887 233828 : for(i=1; i<=k; i++)
2888 : {
2889 190274 : gel(r,i) = cgetg(m+3, t_POL);
2890 190274 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2891 : }
2892 571529 : for (j=1, i=0, l=2; i<=n; i++)
2893 : {
2894 527975 : gmael(r,j,l) = gel(p,2+i);
2895 527975 : if (j==k) { j=1; l++; } else j++;
2896 : }
2897 233828 : for(i=1; i<=k; i++)
2898 190274 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2899 43554 : return r;
2900 : }
2901 :
2902 : /*******************************************************************/
2903 : /* */
2904 : /* Kronecker form */
2905 : /* */
2906 : /*******************************************************************/
2907 :
2908 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2909 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2910 : * normalized coefficients */
2911 : GEN
2912 177355 : Kronecker_to_mod(GEN z, GEN T)
2913 : {
2914 177355 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2915 177355 : GEN x, t = cgetg(N,t_POL);
2916 177355 : t[1] = T[1];
2917 177355 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2918 177355 : x[1] = z[1];
2919 177355 : T = RgX_copy(T);
2920 1314332 : for (i=2; i<lx+2; i++, z+= N-2)
2921 : {
2922 5434342 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2923 1136977 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2924 : }
2925 177355 : N = (l-2) % (N-2) + 2;
2926 456907 : for (j=2; j<N; j++) t[j] = z[j];
2927 177355 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2928 177355 : return normalizepol_lg(x, i+1);
2929 : }
2930 :
2931 : /*******************************************************************/
2932 : /* */
2933 : /* Domain detection */
2934 : /* */
2935 : /*******************************************************************/
2936 :
2937 : static GEN
2938 151600 : zero_FpX_mod(GEN p, long v)
2939 : {
2940 151600 : GEN r = cgetg(3,t_POL);
2941 151600 : r[1] = evalvarn(v);
2942 151600 : gel(r,2) = mkintmod(gen_0, icopy(p));
2943 151600 : return r;
2944 : }
2945 :
2946 : static GEN
2947 797978 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2948 : {
2949 797978 : pari_sp av = avma;
2950 : GEN r;
2951 797978 : if (lgefint(p) == 3)
2952 : {
2953 747530 : ulong pp = uel(p, 2);
2954 747530 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2955 : RgX_to_Flx(y, pp), pp));
2956 : }
2957 : else
2958 50448 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2959 797978 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2960 684094 : return gerepileupto(av, FpX_to_mod(r, p));
2961 : }
2962 :
2963 : static GEN
2964 532 : zero_FpXQX_mod(GEN pol, GEN p, long v)
2965 : {
2966 532 : GEN r = cgetg(3,t_POL);
2967 532 : r[1] = evalvarn(v);
2968 532 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
2969 532 : return r;
2970 : }
2971 :
2972 : static GEN
2973 2044 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
2974 : {
2975 2044 : pari_sp av = avma;
2976 : long dT;
2977 : GEN kx, ky, r;
2978 2044 : GEN T = RgX_to_FpX(pol, p);
2979 2044 : if (signe(T)==0) pari_err_OP("*", x, y);
2980 2037 : dT = degpol(T);
2981 2037 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
2982 2037 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
2983 2037 : r = FpX_mul(kx, ky, p);
2984 2037 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
2985 1512 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
2986 : }
2987 :
2988 : static GEN
2989 412974 : RgX_liftred(GEN x, GEN T)
2990 412974 : { return RgXQX_red(liftpol_shallow(x), T); }
2991 :
2992 : static GEN
2993 162166 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
2994 : {
2995 162166 : pari_sp av = avma;
2996 162166 : long dT = degpol(T);
2997 162166 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
2998 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
2999 162166 : return gerepileupto(av, Kronecker_to_mod(r, T));
3000 : }
3001 :
3002 : static GEN
3003 1365 : RgX_sqr_FpX(GEN x, GEN p)
3004 : {
3005 1365 : pari_sp av = avma;
3006 : GEN r;
3007 1365 : if (lgefint(p) == 3)
3008 : {
3009 1161 : ulong pp = uel(p, 2);
3010 1161 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3011 : }
3012 : else
3013 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3014 1365 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3015 1246 : return gerepileupto(av, FpX_to_mod(r, p));
3016 : }
3017 :
3018 : static GEN
3019 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3020 : {
3021 196 : pari_sp av = avma;
3022 : long dT;
3023 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3024 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3025 189 : dT = degpol(T);
3026 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3027 189 : r = FpX_sqr(kx, p);
3028 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3029 189 : return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3030 : }
3031 :
3032 : static GEN
3033 13488 : RgX_sqr_QXQX(GEN x, GEN T)
3034 : {
3035 13488 : pari_sp av = avma;
3036 13488 : long dT = degpol(T);
3037 13488 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3038 13488 : return gerepileupto(av, Kronecker_to_mod(r, T));
3039 : }
3040 :
3041 : static GEN
3042 67494 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3043 : {
3044 67494 : pari_sp av = avma;
3045 : GEN r;
3046 67494 : if (lgefint(p) == 3)
3047 : {
3048 50704 : ulong pp = uel(p, 2);
3049 50704 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3050 : RgX_to_Flx(y, pp), pp));
3051 : }
3052 : else
3053 16790 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3054 67494 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3055 29897 : return gerepileupto(av, FpX_to_mod(r, p));
3056 : }
3057 :
3058 : static GEN
3059 37577 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3060 : {
3061 37577 : pari_sp av = avma;
3062 : GEN r;
3063 37577 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3064 37577 : return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
3065 : }
3066 : static GEN
3067 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3068 : {
3069 70 : pari_sp av = avma;
3070 : GEN r;
3071 70 : GEN T = RgX_to_FpX(pol, p);
3072 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3073 63 : if (lgefint(p) == 3)
3074 : {
3075 55 : ulong pp = uel(p, 2);
3076 55 : GEN Tp = ZX_to_Flx(T, pp);
3077 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3078 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3079 : }
3080 : else
3081 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3082 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3083 56 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3084 : }
3085 :
3086 : #define code(t1,t2) ((t1 << 6) | t2)
3087 : static GEN
3088 70622596 : RgX_mul_fast(GEN x, GEN y)
3089 : {
3090 : GEN p, pol;
3091 : long pa;
3092 70622596 : long t = RgX_type2(x,y, &p,&pol,&pa);
3093 70622857 : switch(t)
3094 : {
3095 38697129 : case t_INT: return ZX_mul(x,y);
3096 3209594 : case t_FRAC: return QX_mul(x,y);
3097 104665 : case t_FFELT: return FFX_mul(x, y, pol);
3098 797978 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3099 162166 : case code(t_POLMOD, t_INT):
3100 : case code(t_POLMOD, t_FRAC):
3101 162166 : return RgX_mul_QXQX(x, y, pol);
3102 2045 : case code(t_POLMOD, t_INTMOD):
3103 2045 : return RgX_mul_FpXQX(x, y, pol, p);
3104 27649280 : default: return NULL;
3105 : }
3106 : }
3107 : static GEN
3108 1134911 : RgX_sqr_fast(GEN x)
3109 : {
3110 : GEN p, pol;
3111 : long pa;
3112 1134911 : long t = RgX_type(x,&p,&pol,&pa);
3113 1134911 : switch(t)
3114 : {
3115 1094011 : case t_INT: return ZX_sqr(x);
3116 14813 : case t_FRAC: return QX_sqr(x);
3117 2296 : case t_FFELT: return FFX_sqr(x, pol);
3118 1365 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3119 13488 : case code(t_POLMOD, t_INT):
3120 : case code(t_POLMOD, t_FRAC):
3121 13488 : return RgX_sqr_QXQX(x, pol);
3122 196 : case code(t_POLMOD, t_INTMOD):
3123 196 : return RgX_sqr_FpXQX(x, pol, p);
3124 8742 : default: return NULL;
3125 : }
3126 : }
3127 :
3128 : static GEN
3129 10524572 : RgX_rem_fast(GEN x, GEN y)
3130 : {
3131 : GEN p, pol;
3132 : long pa;
3133 10524572 : long t = RgX_type2(x,y, &p,&pol,&pa);
3134 10524903 : switch(t)
3135 : {
3136 3223915 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3137 1666594 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3138 84 : case t_FFELT: return FFX_rem(x, y, pol);
3139 67494 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3140 37577 : case code(t_POLMOD, t_INT):
3141 : case code(t_POLMOD, t_FRAC):
3142 37577 : return RgX_rem_QXQX(x, y, pol);
3143 129 : case code(t_POLMOD, t_INTMOD):
3144 129 : return RgX_rem_FpXQX(x, y, pol, p);
3145 5529110 : default: return NULL;
3146 : }
3147 : }
3148 :
3149 : #undef code
3150 :
3151 : GEN
3152 57103156 : RgX_mul(GEN x, GEN y)
3153 : {
3154 57103156 : GEN z = RgX_mul_fast(x,y);
3155 57103246 : if (!z) z = RgX_mul_i(x,y);
3156 57102635 : return z;
3157 : }
3158 :
3159 : GEN
3160 1122759 : RgX_sqr(GEN x)
3161 : {
3162 1122759 : GEN z = RgX_sqr_fast(x);
3163 1122752 : if (!z) z = RgX_sqr_i(x);
3164 1122752 : return z;
3165 : }
3166 :
3167 : GEN
3168 10524583 : RgX_rem(GEN x, GEN y)
3169 : {
3170 10524583 : GEN z = RgX_rem_fast(x, y);
3171 10524809 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3172 10524734 : return z;
3173 : }
3174 :
3175 : static GEN
3176 0 : _RgX_mul(void* E, GEN x, GEN y)
3177 0 : { (void) E; return RgX_mul(x, y); }
3178 :
3179 : GEN
3180 0 : RgXV_prod(GEN V)
3181 0 : { return gen_product(V, NULL, &_RgX_mul); }
3182 :
3183 : GEN
3184 10893007 : RgXn_mul(GEN f, GEN g, long n)
3185 : {
3186 10893007 : pari_sp av = avma;
3187 10893007 : GEN h = RgX_mul_fast(f,g);
3188 10893008 : if (!h) return RgXn_mul2(f,g,n);
3189 1743411 : if (degpol(h) < n) return h;
3190 362725 : return gerepilecopy(av, RgXn_red_shallow(h, n));
3191 : }
3192 :
3193 : GEN
3194 12152 : RgXn_sqr(GEN f, long n)
3195 : {
3196 12152 : pari_sp av = avma;
3197 12152 : GEN g = RgX_sqr_fast(f);
3198 12152 : if (!g) return RgXn_sqr2(f,n);
3199 11837 : if (degpol(g) < n) return g;
3200 10640 : return gerepilecopy(av, RgXn_red_shallow(g, n));
3201 : }
3202 :
3203 : /* (f*g) \/ x^n */
3204 : GEN
3205 2626459 : RgX_mulhigh_i(GEN f, GEN g, long n)
3206 : {
3207 2626459 : GEN h = RgX_mul_fast(f,g);
3208 2626459 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3209 : }
3210 :
3211 : /* (f*g) \/ x^n */
3212 : GEN
3213 0 : RgX_sqrhigh_i(GEN f, long n)
3214 : {
3215 0 : GEN h = RgX_sqr_fast(f);
3216 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3217 : }
|
