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 19622294 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19622294 : long pold=1, rold=n*(d-1);
35 92318496 : for(p=2; p<=d; p++)
36 : {
37 72696202 : r = m*(p-1) + n*((d-1)/p);
38 72696202 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19622294 : return pold;
41 : }
42 :
43 : static GEN
44 13708713 : 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 13708713 : pari_sp av = avma;
48 : long i;
49 13708713 : GEN z = cmul(E,P,a,ff->one(E));
50 13665393 : if (!z) z = gen_0;
51 69555980 : for (i=1; i<=n; i++)
52 : {
53 55879137 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 55898124 : if (t) {
55 43076386 : z = ff->add(E, z, t);
56 43062680 : if (gc_needed(av,2)) z = gerepileupto(av, z);
57 : }
58 : }
59 13676843 : 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 12260664 : 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 12260664 : pari_sp av = avma;
72 12260664 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12260664 : if (d < 0) return ff->zero(E);
76 11654801 : if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1030490 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1030490 : 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 1030490 : d -= l;
84 1030490 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2051576 : while (d >= l-1)
86 : {
87 1019650 : d -= l-1;
88 1019650 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1019383 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1019660 : if (gc_needed(av,2))
91 105 : z = gerepileupto(av, z);
92 : }
93 1031926 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1031927 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1031928 : return gerepileupto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 731259 : 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 731259 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 731259 : if (d < 0) return ff->zero(E);
106 730314 : rtd = (long) sqrt((double)d);
107 730314 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 730320 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 730319 : return gerepileupto(av, z);
110 : }
111 :
112 : static GEN
113 1194025 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 15702936 : _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 1042127 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 344959 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 1217499 : _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 226006 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
161 : {
162 226006 : pari_sp av = avma;
163 226006 : long i, d = degpol(P);
164 : GEN s;
165 226006 : if (signe(P)==0) return pol_0(varn(P));
166 226006 : s = gel(P, d+2);
167 226006 : if (d == 0) return gcopy(s);
168 1079554 : for (i = d-1; i >= 0; i--)
169 : {
170 856831 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
171 856831 : 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 222723 : return gerepileupto(av, s);
178 : }
179 :
180 : GEN
181 1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
182 : {
183 1652 : pari_sp av = avma;
184 1652 : long i, d = degpol(P), v = varn(A);
185 : GEN s;
186 1652 : if (signe(P)==0) return pol_0(v);
187 1652 : if (d == 0) return scalarpol(gel(P, d+2), v);
188 1232 : s = scalarpol_shallow(gel(P, d+2), v);
189 4963 : for (i = d-1; i >= 0; i--)
190 : {
191 3731 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
192 3731 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
193 3731 : 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 1232 : return gerepileupto(av, s);
200 : }
201 :
202 : const struct bb_algebra *
203 154038 : get_Rg_algebra(void)
204 : {
205 154038 : return &Rg_algebra;
206 : }
207 :
208 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
209 :
210 : static GEN
211 30352 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
212 : {
213 : (void) E;
214 30352 : return RgX_divrem(x, y, r);
215 : }
216 :
217 : GEN
218 6923 : RgX_digits(GEN x, GEN T)
219 : {
220 6923 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
221 6923 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
222 : }
223 :
224 : /*******************************************************************/
225 : /* */
226 : /* RgX */
227 : /* */
228 : /*******************************************************************/
229 :
230 : long
231 25006642 : RgX_equal(GEN x, GEN y)
232 : {
233 25006642 : long i = lg(x);
234 :
235 25006642 : if (i != lg(y)) return 0;
236 107625929 : for (i--; i > 1; i--)
237 82923857 : if (!gequal(gel(x,i),gel(y,i))) return 0;
238 24702072 : 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 156643776 : Rg_get_1(GEN x)
245 : {
246 : GEN p, T;
247 156643776 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
248 156643776 : if (RgX_type_is_composite(tx))
249 11621734 : RgX_type_decode(tx, &i /*junk*/, &tx);
250 156643776 : switch(tx)
251 : {
252 795214 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
253 973 : case t_PADIC: return cvtop(gen_1, p, lx);
254 5481 : case t_FFELT: return FF_1(T);
255 155842108 : 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 6974240 : Rg_get_0(GEN x)
262 : {
263 : GEN p, T;
264 6974240 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
265 6974240 : if (RgX_type_is_composite(tx))
266 51114 : RgX_type_decode(tx, &i /*junk*/, &tx);
267 6974240 : switch(tx)
268 : {
269 525 : 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 6973442 : default: return gen_0;
273 : }
274 : }
275 :
276 : GEN
277 5502 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
278 : {
279 5502 : long i, n = degpol(P);
280 : GEN z, dz, dP;
281 5502 : if (n < 0) return gen_0;
282 5502 : P = Q_remove_denom(P, &dP);
283 5502 : z = gel(P,2); if (n == 0) return icopy(z);
284 3248 : if (dV) z = mulii(dV, z); /* V[1] = dV */
285 3248 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
286 6720 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
287 3248 : dz = mul_denom(dP, dV);
288 3248 : return dz? RgX_Rg_div(z, dz): z;
289 : }
290 :
291 : /* Return P(h * x), not memory clean */
292 : GEN
293 18416 : RgX_unscale(GEN P, GEN h)
294 : {
295 18416 : long i, l = lg(P);
296 18416 : GEN hi = gen_1, Q = cgetg(l, t_POL);
297 18416 : Q[1] = P[1];
298 18416 : if (l == 2) return Q;
299 18388 : gel(Q,2) = gcopy(gel(P,2));
300 39624 : for (i=3; i<l; i++)
301 : {
302 21236 : hi = gmul(hi,h);
303 21236 : gel(Q,i) = gmul(gel(P,i), hi);
304 : }
305 18388 : return Q;
306 : }
307 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
308 : GEN
309 1450481 : ZX_z_unscale(GEN P, long h)
310 : {
311 1450481 : long i, l = lg(P);
312 1450481 : GEN Q = cgetg(l, t_POL);
313 1450481 : Q[1] = P[1];
314 1450481 : if (l == 2) return Q;
315 1362047 : gel(Q,2) = gel(P,2);
316 1362047 : if (l == 3) return Q;
317 1335727 : if (h == -1)
318 228793 : for (i = 3; i < l; i++)
319 : {
320 189946 : gel(Q,i) = negi(gel(P,i));
321 189945 : if (++i == l) break;
322 140337 : gel(Q,i) = gel(P,i);
323 : }
324 : else
325 : {
326 : GEN hi;
327 1247271 : gel(Q,3) = mulis(gel(P,3), h);
328 1247269 : hi = sqrs(h);
329 5360804 : for (i = 4; i < l; i++)
330 : {
331 4113533 : gel(Q,i) = mulii(gel(P,i), hi);
332 4113535 : if (i != l-1) hi = mulis(hi,h);
333 : }
334 : }
335 1335726 : 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 999490 : ZX_unscale(GEN P, GEN h)
340 : {
341 : long i, l;
342 : GEN Q, hi;
343 999490 : 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 216843 : ZX_unscale2n(GEN P, long n)
361 : {
362 216843 : long i, ni = n, l = lg(P);
363 216843 : GEN Q = cgetg(l, t_POL);
364 216843 : Q[1] = P[1];
365 216843 : if (l == 2) return Q;
366 216843 : gel(Q,2) = gel(P,2);
367 216843 : if (l == 3) return Q;
368 216843 : gel(Q,3) = shifti(gel(P,3), ni);
369 798293 : for (i=4; i<l; i++)
370 : {
371 581480 : ni += n;
372 581480 : gel(Q,i) = shifti(gel(P,i), ni);
373 : }
374 216813 : return Q;
375 : }
376 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
377 : GEN
378 12073 : ZX_unscale_div(GEN P, GEN h)
379 : {
380 12073 : long i, l = lg(P);
381 12073 : GEN hi, Q = cgetg(l, t_POL);
382 12073 : Q[1] = P[1];
383 12073 : if (l == 2) return Q;
384 12073 : gel(Q,2) = diviiexact(gel(P,2), h);
385 12073 : if (l == 3) return Q;
386 12073 : gel(Q,3) = gel(P,3);
387 12073 : if (l == 4) return Q;
388 12073 : hi = h;
389 12073 : gel(Q,4) = mulii(gel(P,4), hi);
390 61349 : for (i=5; i<l; i++)
391 : {
392 49276 : hi = mulii(hi,h);
393 49276 : gel(Q,i) = mulii(gel(P,i), hi);
394 : }
395 12073 : 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 5754 : RgXV_unscale(GEN x, GEN h)
415 : {
416 5754 : if (isint1(h)) return RgX_copy(x);
417 16516 : 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 4141244 : RgX_rescale(GEN P, GEN h)
423 : {
424 4141244 : long i, l = lg(P);
425 4141244 : GEN Q = cgetg(l,t_POL), hi = h;
426 4141205 : gel(Q,l-1) = gel(P,l-1);
427 10969669 : for (i=l-2; i>=2; i--)
428 : {
429 10967065 : gel(Q,i) = gmul(gel(P,i), hi);
430 10966872 : if (i == 2) break;
431 6828069 : hi = gmul(hi,h);
432 : }
433 4141407 : 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 959256 : RgX_deflate(GEN x0, long d)
446 : {
447 : GEN z, y, x;
448 959256 : long i,id, dy, dx = degpol(x0);
449 959253 : if (d == 1 || dx <= 0) return leafcopy(x0);
450 371263 : dy = dx/d;
451 371263 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
452 371259 : z = y + 2;
453 371259 : x = x0+ 2;
454 1493924 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
455 371259 : return y;
456 : }
457 :
458 : GEN
459 18732 : RgX_homogenize(GEN P, long v)
460 : {
461 : long i, l, d;
462 18732 : GEN Q = cgetg_copy(P, &l);
463 18732 : Q[1] = P[1]; d = l-3;
464 165998 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
465 18732 : 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 796647 : RgX_inflate(GEN x0, long d)
499 : {
500 796647 : long i, id, dy, dx = degpol(x0);
501 796647 : GEN x = x0 + 2, z, y;
502 796647 : if (dx <= 0) return leafcopy(x0);
503 789045 : dy = dx*d;
504 789045 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
505 789044 : z = y + 2;
506 27876042 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
507 11420033 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
508 789044 : return y;
509 : }
510 :
511 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
512 : static GEN
513 5804061 : RgX_translate_basecase(GEN P, GEN c)
514 : {
515 5804061 : pari_sp av = avma;
516 : GEN Q, R;
517 : long i, k, n;
518 :
519 5804061 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
520 5802017 : Q = leafcopy(P);
521 5802026 : R = Q+2; n = degpol(P);
522 5802027 : if (isint1(c))
523 : {
524 26090 : for (i=1; i<=n; i++)
525 : {
526 134618 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
527 22255 : 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 5798184 : else if (isintm1(c))
535 : {
536 891459 : for (i=1; i<=n; i++)
537 : {
538 2825510 : for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
539 676534 : 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 19103508 : for (i=1; i<=n; i++)
549 : {
550 44897110 : for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
551 13520250 : 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 5801388 : return gerepilecopy(av, Q);
559 : }
560 : GEN
561 5804490 : RgX_translate(GEN P, GEN c)
562 : {
563 5804490 : pari_sp av = avma;
564 5804490 : long n = degpol(P);
565 5804480 : if (n < 40)
566 5804060 : return RgX_translate_basecase(P, c);
567 : else
568 : {
569 420 : long d = n >> 1;
570 420 : GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
571 420 : GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
572 420 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
573 420 : 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 23651 : RgXQX_translate(GEN P, GEN c, GEN T)
588 : {
589 23651 : pari_sp av = avma;
590 : GEN Q, R;
591 : long i, k, n;
592 :
593 23651 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
594 23322 : Q = leafcopy(P);
595 23322 : R = Q+2; n = degpol(P);
596 82367 : for (i=1; i<=n; i++)
597 : {
598 260749 : for (k=n-i; k<n; k++)
599 : {
600 201704 : pari_sp av2 = avma;
601 201704 : gel(R,k) = gerepileupto(av2,
602 201704 : RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
603 : }
604 59045 : 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 23322 : 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 166649 : QXQ_to_mod(GEN x, GEN T)
623 : {
624 : long d;
625 166649 : switch(typ(x))
626 : {
627 62304 : case t_INT: return icopy(x);
628 2079 : case t_FRAC: return gcopy(x);
629 102266 : case t_POL:
630 102266 : d = degpol(x);
631 102266 : if (d < 0) return gen_0;
632 99998 : if (d == 0) return gcopy(gel(x,2));
633 97859 : 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 709351 : QXQ_to_mod_shallow(GEN x, GEN T)
641 : {
642 : long d;
643 709351 : switch(typ(x))
644 : {
645 454645 : case t_INT:
646 454645 : case t_FRAC: return x;
647 254706 : case t_POL:
648 254706 : d = degpol(x);
649 254706 : if (d < 0) return gen_0;
650 209296 : if (d == 0) return gel(x,2);
651 194294 : 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 37184 : QXQX_to_mod(GEN z, GEN T)
660 : {
661 37184 : long i,l = lg(z);
662 37184 : GEN x = cgetg(l,t_POL);
663 185808 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
664 37184 : x[1] = z[1]; return normalizepol_lg(x,l);
665 : }
666 : /* pure shallow version */
667 : GEN
668 161699 : QXQX_to_mod_shallow(GEN z, GEN T)
669 : {
670 161699 : long i,l = lg(z);
671 161699 : GEN x = cgetg(l,t_POL);
672 790270 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
673 161699 : x[1] = z[1]; return normalizepol_lg(x,l);
674 : }
675 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
676 : GEN
677 12495 : QXQXV_to_mod(GEN V, GEN T)
678 : {
679 12495 : long i, l = lg(V);
680 12495 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
681 49679 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
682 12495 : return z;
683 : }
684 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
685 : GEN
686 18177 : QXQV_to_mod(GEN V, GEN T)
687 : {
688 18177 : long i, l = lg(V);
689 18177 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
690 36202 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
691 18177 : return z;
692 : }
693 :
694 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
695 : GEN
696 13412 : QXQC_to_mod_shallow(GEN V, GEN T)
697 : {
698 13412 : long i, l = lg(V);
699 13412 : GEN z = cgetg(l, t_COL);
700 94192 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
701 13412 : return z;
702 : }
703 :
704 : GEN
705 6251 : QXQM_to_mod_shallow(GEN V, GEN T)
706 : {
707 6251 : long i, l = lg(V);
708 6251 : GEN z = cgetg(l, t_MAT);
709 19663 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
710 6251 : return z;
711 : }
712 :
713 : GEN
714 6142809 : RgX_renormalize_lg(GEN x, long lx)
715 : {
716 : long i;
717 9156776 : for (i = lx-1; i>1; i--)
718 8707987 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
719 6142809 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
720 6142809 : setlg(x, i+1); setsigne(x, i != 1); return x;
721 : }
722 :
723 : GEN
724 1488421 : RgV_to_RgX(GEN x, long v)
725 : {
726 1488421 : long i, k = lg(x);
727 : GEN p;
728 :
729 4137070 : while (--k && gequal0(gel(x,k)));
730 1488421 : if (!k) return pol_0(v);
731 1481577 : i = k+2; p = cgetg(i,t_POL);
732 1481579 : p[1] = evalsigne(1) | evalvarn(v);
733 10206062 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
734 1481579 : return p;
735 : }
736 : GEN
737 189375 : RgV_to_RgX_reverse(GEN x, long v)
738 : {
739 189375 : long j, k, l = lg(x);
740 : GEN p;
741 :
742 190852 : for (k = 1; k < l; k++)
743 190852 : if (!gequal0(gel(x,k))) break;
744 189375 : if (k == l) return pol_0(v);
745 189375 : k -= 1;
746 189375 : l -= k;
747 189375 : x += k;
748 189375 : p = cgetg(l+1,t_POL);
749 189375 : p[1] = evalsigne(1) | evalvarn(v);
750 993263 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
751 189375 : return p;
752 : }
753 :
754 : /* return the (N-dimensional) vector of coeffs of p */
755 : GEN
756 13135351 : RgX_to_RgC(GEN x, long N)
757 : {
758 : long i, l;
759 : GEN z;
760 13135351 : l = lg(x)-1; x++;
761 13135351 : if (l > N+1) l = N+1; /* truncate higher degree terms */
762 13135351 : z = cgetg(N+1,t_COL);
763 79858227 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
764 24036810 : for ( ; i<=N; i++) gel(z,i) = gen_0;
765 13135519 : return z;
766 : }
767 : GEN
768 1366249 : Rg_to_RgC(GEN x, long N)
769 : {
770 1366249 : 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 300520 : RgM_to_RgXV(GEN x, long v)
776 1281899 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
777 : GEN
778 7097 : RgM_to_RgXV_reverse(GEN x, long v)
779 28388 : { 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 339721 : RgV_to_RgM(GEN x, long n)
785 1703477 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
786 :
787 : GEN
788 77490 : RgXV_to_RgM(GEN x, long n)
789 392193 : { 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 23116 : RgM_to_RgXX(GEN x, long v,long w)
794 : {
795 23116 : long j, lx = lg(x);
796 23116 : GEN y = cgetg(lx+1, t_POL);
797 23116 : y[1] = evalsigne(1) | evalvarn(v);
798 23116 : y++;
799 127478 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
800 23116 : 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 322 : RgXX_to_RgM(GEN v, long n)
807 : {
808 322 : long j, N = lg(v)-1;
809 322 : GEN y = cgetg(N, t_MAT);
810 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
811 322 : return y;
812 : }
813 :
814 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
815 : GEN
816 31244 : RgXY_swapspec(GEN x, long n, long w, long nx)
817 : {
818 31244 : long j, ly = n+3;
819 31244 : GEN y = cgetg(ly, t_POL);
820 31244 : y[1] = evalsigne(1);
821 384620 : for (j=2; j<ly; j++)
822 : {
823 : long k;
824 353376 : GEN a = cgetg(nx+2,t_POL);
825 353376 : a[1] = evalsigne(1) | evalvarn(w);
826 1809209 : for (k=0; k<nx; k++)
827 : {
828 1455833 : GEN xk = gel(x,k);
829 1455833 : if (typ(xk)==t_POL)
830 1367317 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
831 : else
832 88516 : gel(a,k+2) = j==2 ? xk: gen_0;
833 : }
834 353376 : gel(y,j) = normalizepol_lg(a, nx+2);
835 : }
836 31244 : 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 288 : RgXY_degreex(GEN b)
849 : {
850 288 : long deg = 0, i;
851 288 : if (!signe(b)) return -1;
852 1396 : for (i = 2; i < lg(b); ++i)
853 : {
854 1108 : GEN bi = gel(b, i);
855 1108 : if (typ(bi) == t_POL)
856 1030 : deg = maxss(deg, degpol(bi));
857 : }
858 288 : return deg;
859 : }
860 :
861 : GEN
862 32578 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
863 :
864 : /* return (x % X^n). Shallow */
865 : GEN
866 7899548 : RgXn_red_shallow(GEN a, long n)
867 : {
868 7899548 : long i, L = n+2, l = lg(a);
869 : GEN b;
870 7899548 : if (L >= l) return a; /* deg(x) < n */
871 5715408 : b = cgetg(L, t_POL); b[1] = a[1];
872 36532441 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
873 5715408 : 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 164757571 : RgX_shift_shallow(GEN a, long n)
883 : {
884 164757571 : long i, l = lg(a);
885 : GEN b;
886 164757571 : if (l == 2 || !n) return a;
887 107860569 : l += n;
888 107860569 : if (n < 0)
889 : {
890 55062915 : if (l <= 2) return pol_0(varn(a));
891 53599045 : b = cgetg(l, t_POL); b[1] = a[1];
892 53599260 : a -= n;
893 163687977 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
894 : } else {
895 52797654 : b = cgetg(l, t_POL); b[1] = a[1];
896 52797794 : a -= n; n += 2;
897 114926143 : for (i=2; i<n; i++) gel(b,i) = gen_0;
898 204240900 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
899 : }
900 106397054 : return b;
901 : }
902 : /* return (x * X^n). */
903 : GEN
904 1567501 : RgX_shift(GEN a, long n)
905 : {
906 1567501 : long i, l = lg(a);
907 : GEN b;
908 1567501 : if (l == 2 || !n) return RgX_copy(a);
909 1567137 : l += n;
910 1567137 : 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 1566290 : b = cgetg(l, t_POL); b[1] = a[1];
918 1566290 : a -= n; n += 2;
919 3915502 : for (i=2; i<n; i++) gel(b,i) = gen_0;
920 4318398 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
921 : }
922 1567067 : 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 2984854 : RgX_mulXn(GEN x, long d)
943 : {
944 : pari_sp av;
945 : GEN z;
946 : long v;
947 2984854 : if (d >= 0) return RgX_shift(x, d);
948 1472237 : d = -d;
949 1472237 : v = RgX_val(x);
950 1472237 : if (v >= d) return RgX_shift(x, -d);
951 1472230 : av = avma;
952 1472230 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
953 1472230 : 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 3627123 : RgX_val(GEN x)
967 : {
968 3627123 : long i, lx = lg(x);
969 3627123 : if (lx == 2) return LONG_MAX;
970 4555176 : for (i = 2; i < lx; i++)
971 4555120 : if (!isexactzero(gel(x,i))) break;
972 3626969 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
973 3626913 : return i - 2;
974 : }
975 : long
976 80181178 : RgX_valrem(GEN x, GEN *Z)
977 : {
978 80181178 : long v, i, lx = lg(x);
979 80181178 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
980 125524164 : for (i = 2; i < lx; i++)
981 125530107 : if (!isexactzero(gel(x,i))) break;
982 : /* possible with nonrational zeros */
983 80195083 : if (i == lx)
984 : {
985 14 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
986 14 : return LONG_MAX;
987 : }
988 80195069 : v = i - 2;
989 80195069 : *Z = RgX_shift_shallow(x, -v);
990 80209413 : return v;
991 : }
992 : long
993 831798 : RgX_valrem_inexact(GEN x, GEN *Z)
994 : {
995 : long v;
996 831798 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
997 847559 : for (v = 0;; v++)
998 847559 : if (!gequal0(gel(x,2+v))) break;
999 831791 : if (Z) *Z = RgX_shift_shallow(x, -v);
1000 831791 : return v;
1001 : }
1002 :
1003 : GEN
1004 63875 : RgXQC_red(GEN x, GEN T)
1005 394877 : { 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 12257 : RgXQM_red(GEN x, GEN T)
1013 76132 : { 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 456154 : RgXQX_red(GEN P, GEN T)
1023 : {
1024 456154 : long i, l = lg(P);
1025 456154 : GEN Q = cgetg(l, t_POL);
1026 456154 : Q[1] = P[1];
1027 2383598 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1028 456154 : return normalizepol_lg(Q, l);
1029 : }
1030 :
1031 : GEN
1032 793024 : RgX_deriv(GEN x)
1033 : {
1034 793024 : long i,lx = lg(x)-1;
1035 : GEN y;
1036 :
1037 793024 : if (lx<3) return pol_0(varn(x));
1038 790329 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1039 3516843 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1040 790326 : y[1] = x[1]; return normalizepol_lg(y,i);
1041 : }
1042 :
1043 : GEN
1044 1378109 : RgX_recipspec_shallow(GEN x, long l, long n)
1045 : {
1046 : long i;
1047 1378109 : GEN z = cgetg(n+2,t_POL);
1048 1378117 : z[1] = 0; z += 2;
1049 36196216 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1050 1601221 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1051 1378117 : return normalizepol_lg(z-2,n+2);
1052 : }
1053 :
1054 : GEN
1055 636818 : RgXn_recip_shallow(GEN P, long n)
1056 : {
1057 636818 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1058 636839 : setvarn(Q, varn(P));
1059 636839 : return Q;
1060 : }
1061 :
1062 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1063 : GEN
1064 31150 : RgX_recip(GEN x)
1065 : {
1066 : long lx, i, j;
1067 31150 : GEN y = cgetg_copy(x, &lx);
1068 272104 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1069 31150 : return normalizepol_lg(y,lx);
1070 : }
1071 : /* shallow version */
1072 : GEN
1073 57526 : RgX_recip_shallow(GEN x)
1074 : {
1075 : long lx, i, j;
1076 57526 : GEN y = cgetg_copy(x, &lx);
1077 344995 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1078 57526 : return normalizepol_lg(y,lx);
1079 : }
1080 :
1081 : GEN
1082 3626492 : RgX_recip_i(GEN x)
1083 : {
1084 : long lx, i, j;
1085 3626492 : GEN y = cgetg_copy(x, &lx);
1086 19192861 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1087 3626503 : return y;
1088 : }
1089 : /*******************************************************************/
1090 : /* */
1091 : /* ADDITION / SUBTRACTION */
1092 : /* */
1093 : /*******************************************************************/
1094 : /* same variable */
1095 : GEN
1096 76776833 : RgX_add(GEN x, GEN y)
1097 : {
1098 76776833 : long i, lx = lg(x), ly = lg(y);
1099 : GEN z;
1100 76776833 : if (ly <= lx) {
1101 66743712 : z = cgetg(lx,t_POL); z[1] = x[1];
1102 264583035 : for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1103 118575028 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1104 66696893 : z = normalizepol_lg(z, lx);
1105 : } else {
1106 10033121 : z = cgetg(ly,t_POL); z[1] = y[1];
1107 43108089 : for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1108 31062328 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1109 10034408 : z = normalizepol_lg(z, ly);
1110 : }
1111 76774159 : return z;
1112 : }
1113 : GEN
1114 61869497 : RgX_sub(GEN x, GEN y)
1115 : {
1116 61869497 : long i, lx = lg(x), ly = lg(y);
1117 : GEN z;
1118 61869497 : if (ly <= lx) {
1119 29665602 : z = cgetg(lx,t_POL); z[1] = x[1];
1120 155391063 : for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1121 55307596 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1122 29634924 : z = normalizepol_lg(z, lx);
1123 : } else {
1124 32203895 : z = cgetg(ly,t_POL); z[1] = y[1];
1125 120188190 : for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1126 71515410 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1127 32159029 : z = normalizepol_lg(z, ly);
1128 : }
1129 61865645 : return z;
1130 : }
1131 : GEN
1132 5171531 : RgX_neg(GEN x)
1133 : {
1134 5171531 : long i, lx = lg(x);
1135 5171531 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
1136 40222294 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
1137 5171278 : return y;
1138 : }
1139 :
1140 : GEN
1141 18863746 : RgX_Rg_add(GEN y, GEN x)
1142 : {
1143 : GEN z;
1144 18863746 : long lz = lg(y), i;
1145 18863746 : if (lz == 2) return scalarpol(x,varn(y));
1146 15825029 : z = cgetg(lz,t_POL); z[1] = y[1];
1147 15824973 : gel(z,2) = gadd(gel(y,2),x);
1148 61630080 : 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 15825110 : return normalizepol_lg(z,lz);
1152 : }
1153 : GEN
1154 64688 : RgX_Rg_add_shallow(GEN y, GEN x)
1155 : {
1156 : GEN z;
1157 64688 : long lz = lg(y), i;
1158 64688 : if (lz == 2) return scalarpol(x,varn(y));
1159 64688 : z = cgetg(lz,t_POL); z[1] = y[1];
1160 64688 : gel(z,2) = gadd(gel(y,2),x);
1161 129525 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1162 64688 : return normalizepol_lg(z,lz);
1163 : }
1164 : GEN
1165 182785 : RgX_Rg_sub(GEN y, GEN x)
1166 : {
1167 : GEN z;
1168 182785 : long lz = lg(y), i;
1169 182785 : 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 180237 : z = cgetg(lz,t_POL); z[1] = y[1];
1179 180234 : gel(z,2) = gsub(gel(y,2),x);
1180 470291 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1181 180240 : return normalizepol_lg(z,lz);
1182 : }
1183 : GEN
1184 4674666 : Rg_RgX_sub(GEN x, GEN y)
1185 : {
1186 : GEN z;
1187 4674666 : long lz = lg(y), i;
1188 4674666 : if (lz == 2) return scalarpol(x,varn(y));
1189 4673588 : z = cgetg(lz,t_POL); z[1] = y[1];
1190 4673548 : gel(z,2) = gsub(x, gel(y,2));
1191 7063026 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1192 4673634 : 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 8504418 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1223 : {
1224 : GEN z, t;
1225 : long i;
1226 8504418 : if (nx == ny) {
1227 1326018 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1228 4238244 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1229 1326019 : return normalizepol_lg(z, nx+2);
1230 : }
1231 7178400 : if (ny < nx) {
1232 6991877 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1233 25831291 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1234 16784971 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1235 6991380 : return normalizepol_lg(z, nx+2);
1236 : } else {
1237 186523 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1238 3674609 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1239 451009 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1240 186552 : return normalizepol_lg(z, ny+2);
1241 : }
1242 : }
1243 : GEN
1244 222771 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1245 : {
1246 : GEN z, t;
1247 : long i;
1248 222771 : 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 209961 : if (ny < nx) {
1254 208183 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1255 3723489 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1256 2377381 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1257 208183 : 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 55567727 : RgXspec_kill0(GEN x, long lx)
1273 : {
1274 55567727 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1275 : long i;
1276 189226316 : for (i=0; i <lx; i++)
1277 : {
1278 133658653 : GEN c = gel(x,i);
1279 133658653 : z[i] = (long)(isrationalzero(c)? NULL: c);
1280 : }
1281 55567663 : return z;
1282 : }
1283 :
1284 : INLINE GEN
1285 92976915 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1286 : {
1287 92976915 : pari_sp av = avma;
1288 92976915 : GEN s = NULL;
1289 : long i;
1290 :
1291 354091966 : for (i=a; i<b; i++)
1292 261124603 : if (gel(y,i) && gel(x,-i))
1293 : {
1294 183091710 : GEN t = gmul(gel(y,i), gel(x,-i));
1295 183086250 : s = s? gadd(s, t): t;
1296 : }
1297 92967363 : return s? gerepileupto(av, s): gen_0;
1298 : }
1299 :
1300 : /* assume nx >= ny > 0, return x * y * t^v */
1301 : static GEN
1302 21192777 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1303 : {
1304 : long i, lz, nz;
1305 : GEN z;
1306 :
1307 21192777 : x = RgXspec_kill0(x,nx);
1308 21192712 : y = RgXspec_kill0(y,ny);
1309 21192760 : lz = nx + ny + 1; nz = lz-2;
1310 21192760 : lz += v;
1311 21192760 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1312 32926118 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1313 56298652 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1314 40420568 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1315 35105951 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1316 21192391 : 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 9065474 : 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 9065474 : if (!signe(x0)) return y0;
1327 8959891 : ny = lgpol(y0);
1328 8959884 : nx = lgpol(x0);
1329 8960003 : zd = (GEN)avma;
1330 8960003 : x = x0 + 2; y = y0 + 2; a = ny-d;
1331 8960003 : if (a <= 0)
1332 : {
1333 1412744 : lz = nx+d+2;
1334 1412744 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1335 3264367 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1336 1412749 : x = zd + a;
1337 1430338 : while (zd > x) gel(--zd,0) = gen_0;
1338 : }
1339 : else
1340 : {
1341 7547259 : xd = new_chunk(d); yd = y+d;
1342 7547272 : x = RgX_addspec_shallow(x,yd, nx,a);
1343 7547251 : lz = (a>nx)? ny+2: lg(x)+d;
1344 38091865 : x += 2; while (xd > x) *--zd = *--xd;
1345 : }
1346 21980998 : while (yd > y) *--zd = *--yd;
1347 8960000 : *--zd = x0[1];
1348 8960000 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1349 : }
1350 : GEN
1351 515330 : RgX_addmulXn(GEN x0, GEN y0, long d)
1352 : {
1353 : GEN x, y, xd, yd, zd;
1354 : long a, lz, nx, ny;
1355 :
1356 515330 : if (!signe(x0)) return RgX_copy(y0);
1357 514560 : nx = lgpol(x0);
1358 514560 : ny = lgpol(y0);
1359 514560 : zd = (GEN)avma;
1360 514560 : x = x0 + 2; y = y0 + 2; a = ny-d;
1361 514560 : 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 222771 : xd = new_chunk(d); yd = y+d;
1372 222771 : x = RgX_addspec(x,yd, nx,a);
1373 222771 : lz = (a>nx)? ny+2: lg(x)+d;
1374 8417559 : x += 2; while (xd > x) *--zd = *--xd;
1375 : }
1376 2619546 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1377 514560 : *--zd = x0[1];
1378 514560 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1379 : }
1380 :
1381 : /* return x * y mod t^n */
1382 : static GEN
1383 6586820 : RgXn_mul_basecase(GEN x, GEN y, long n)
1384 : {
1385 6586820 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1386 : GEN z;
1387 6586820 : if (lx < 0) return pol_0(varn(x));
1388 6586820 : if (ly < 0) return pol_0(varn(x));
1389 6586820 : z = cgetg(lz, t_POL) + 2;
1390 6586820 : x+=2; if (lx > n) lx = n;
1391 6586820 : y+=2; if (ly > n) ly = n;
1392 6586820 : z[-1] = x[-1];
1393 6586820 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1394 6586820 : x = RgXspec_kill0(x, lx);
1395 6586820 : y = RgXspec_kill0(y, ly);
1396 : /* x:y:z [i] = term of degree i */
1397 26091200 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1398 11766558 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1399 6632787 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1400 6586820 : return normalizepol_lg(z - 2, lz);
1401 : }
1402 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1403 : static GEN
1404 9357149 : RgXn_mul2(GEN f, GEN g, long n)
1405 : {
1406 9357149 : pari_sp av = avma;
1407 : GEN fe,fo, ge,go, l,h,m;
1408 : long n0, n1;
1409 9357149 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1410 6622240 : 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 1579843 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1432 : {
1433 1579843 : long d = degpol(f)+degpol(g) + 1 - n;
1434 : GEN h;
1435 1579843 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1436 30067 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1437 : RgX_recip_i(g), d));
1438 30067 : 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 21920781 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1458 : {
1459 : GEN a0, c, c0;
1460 21920781 : long n0, n0a, i, v = 0;
1461 : pari_sp av;
1462 :
1463 28304301 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1464 27804743 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1465 21920760 : if (na < nb) swapspec(a,b, na,nb);
1466 21920760 : if (!nb) return pol_0(0);
1467 :
1468 21672602 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1469 479827 : RgX_shift_inplace_init(v);
1470 479871 : i = (na>>1); n0 = na-i; na = i;
1471 479871 : av = avma; a0 = a+n0; n0a = n0;
1472 1335387 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1473 :
1474 479871 : if (nb > n0)
1475 : {
1476 : GEN b0,c1,c2;
1477 : long n0b;
1478 :
1479 478562 : nb -= n0; b0 = b+n0; n0b = n0;
1480 1446541 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1481 478562 : c = RgX_mulspec(a,b,n0a,n0b);
1482 478562 : c0 = RgX_mulspec(a0,b0, na,nb);
1483 :
1484 478562 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1485 478562 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1486 :
1487 478562 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1488 478562 : c2 = RgX_sub(c1, RgX_add(c0,c));
1489 478562 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1490 : }
1491 : else
1492 : {
1493 1309 : c = RgX_mulspec(a,b,n0a,nb);
1494 1309 : c0 = RgX_mulspec(a0,b,na,nb);
1495 : }
1496 479871 : c0 = RgX_addmulXn(c0,c,n0);
1497 479871 : return RgX_shift_inplace(gerepileupto(av,c0), v);
1498 : }
1499 :
1500 : INLINE GEN
1501 50354 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1502 : {
1503 50354 : pari_sp av = avma;
1504 50354 : GEN s = NULL;
1505 50354 : long j, l = (i+1)>>1;
1506 139496 : for (j=a; j<l; j++)
1507 : {
1508 89142 : GEN xj = gel(x,j), xx = gel(x,i-j);
1509 89142 : if (xj && xx)
1510 : {
1511 81498 : GEN t = gmul(xj, xx);
1512 81498 : s = s? gadd(s, t): t;
1513 : }
1514 : }
1515 50354 : if (s) s = gshift(s,1);
1516 50354 : if ((i&1) == 0)
1517 : {
1518 29603 : GEN t = gel(x, i>>1);
1519 29603 : if (t) {
1520 27209 : t = gsqr(t);
1521 27209 : s = s? gadd(s, t): t;
1522 : }
1523 : }
1524 50354 : return s? gerepileupto(av,s): gen_0;
1525 : }
1526 : static GEN
1527 8852 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1528 : {
1529 : long i, lz, nz;
1530 : GEN z;
1531 :
1532 8852 : if (!nx) return pol_0(0);
1533 8852 : x = RgXspec_kill0(x,nx);
1534 8852 : lz = (nx << 1) + 1, nz = lz-2;
1535 8852 : lz += v;
1536 8852 : z = cgetg(lz,t_POL) + 2;
1537 13024 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1538 38455 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1539 29603 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1540 8852 : 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 8891 : RgX_sqrspec(GEN a, long na)
1586 : {
1587 : GEN a0, c, c0, c1;
1588 8891 : long n0, n0a, i, v = 0;
1589 : pari_sp av;
1590 :
1591 10977 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1592 8891 : 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 1691304 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1609 : {
1610 1691304 : GEN z = RgX_mul(A, B);
1611 1691283 : if (a < b)
1612 10413 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1613 1680870 : else if (a > b)
1614 1099477 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1615 : else
1616 581393 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1617 1691297 : return z;
1618 : }
1619 :
1620 : GEN
1621 20482508 : RgX_mul_i(GEN x, GEN y)
1622 : {
1623 20482508 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1624 20482093 : setvarn(z, varn(x)); return z;
1625 : }
1626 :
1627 : GEN
1628 8813 : RgX_sqr_i(GEN x)
1629 : {
1630 8813 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1631 8813 : setvarn(z,varn(x)); return z;
1632 : }
1633 :
1634 : /*******************************************************************/
1635 : /* */
1636 : /* DIVISION */
1637 : /* */
1638 : /*******************************************************************/
1639 : GEN
1640 3950903 : RgX_Rg_divexact(GEN x, GEN y) {
1641 3950903 : long i, lx = lg(x);
1642 : GEN z;
1643 3950903 : if (lx == 2) return gcopy(x);
1644 3901143 : switch(typ(y))
1645 : {
1646 3871332 : case t_INT:
1647 3871332 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1648 3696617 : break;
1649 4655 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1650 : }
1651 3721773 : z = cgetg(lx, t_POL); z[1] = x[1];
1652 23785845 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1653 3711811 : return z;
1654 : }
1655 : GEN
1656 29122000 : RgX_Rg_div(GEN x, GEN y) {
1657 29122000 : long i, lx = lg(x);
1658 : GEN z;
1659 29122000 : if (lx == 2) return gcopy(x);
1660 28861963 : switch(typ(y))
1661 : {
1662 19672401 : case t_INT:
1663 19672401 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1664 3709091 : break;
1665 12929 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1666 : }
1667 12885724 : z = cgetg(lx, t_POL); z[1] = x[1];
1668 47748208 : for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1669 12885399 : return normalizepol_lg(z, lx);
1670 : }
1671 : GEN
1672 18998 : RgX_normalize(GEN x)
1673 : {
1674 18998 : GEN z, d = NULL;
1675 18998 : long i, n = lg(x)-1;
1676 18998 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1677 18998 : if (i == 1) return pol_0(varn(x));
1678 18998 : if (i == n && isint1(d)) return x;
1679 7322 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1680 16177 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1681 7322 : gel(z,n) = Rg_get_1(d); return z;
1682 : }
1683 : GEN
1684 10423 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1685 : GEN
1686 243161 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1687 : {
1688 243161 : long l = lg(a), i;
1689 : GEN a0, z0, z;
1690 :
1691 243161 : if (l <= 3)
1692 : {
1693 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1694 0 : return pol_0(varn(a));
1695 : }
1696 243161 : z = cgetg(l-1, t_POL);
1697 243164 : z[1] = a[1];
1698 243164 : a0 = a + l-1;
1699 243164 : z0 = z + l-2; *z0 = *a0--;
1700 2928674 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1701 : {
1702 2685513 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1703 2685510 : gel(z0,0) = t;
1704 : }
1705 243161 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1706 243161 : 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 22877887 : 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 22877887 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1722 :
1723 22877887 : dy = degpol(y);
1724 22877849 : y_lead = gel(y,dy+2);
1725 22877849 : 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 22877816 : if (!dy) /* y is constant */
1735 : {
1736 5373 : if (pr == ONLY_REM) return pol_0(varn(x));
1737 5366 : z = RgX_Rg_div(x, y_lead);
1738 5366 : if (pr == ONLY_DIVIDES) return z;
1739 2580 : if (pr) *pr = pol_0(varn(x));
1740 2580 : return z;
1741 : }
1742 22872443 : dx = degpol(x);
1743 22872402 : if (dx < dy)
1744 : {
1745 3696540 : if (pr == ONLY_REM) return RgX_copy(x);
1746 374610 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1747 374589 : z = pol_0(varn(x));
1748 374589 : if (pr) *pr = RgX_copy(x);
1749 374589 : return z;
1750 : }
1751 :
1752 : /* x,y in R[X], y non constant */
1753 19175862 : av = avma;
1754 19175862 : p = NULL;
1755 19175862 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1756 : {
1757 238693 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1758 238693 : if (!z) return gc_NULL(av);
1759 238693 : z = FpX_to_mod(z, p);
1760 238693 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1761 122458 : return gerepileupto(av, z);
1762 116235 : *pr = FpX_to_mod(*pr, p);
1763 116235 : return gc_all(av, 2, &z, pr);
1764 : }
1765 18937335 : 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 2137 : case t_INTMOD:
1772 2137 : case t_POLMOD: y_lead = ginv(y_lead);
1773 2137 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1774 2137 : break;
1775 18935198 : default: if (gequal1(y_lead)) y_lead = NULL;
1776 18935189 : f = gdiv; mod = NULL;
1777 : }
1778 :
1779 18937326 : if (y_lead == NULL)
1780 17113060 : p2 = gel(x,dx+2);
1781 : else {
1782 : for(;;) {
1783 1824266 : p2 = f(gel(x,dx+2),y_lead);
1784 1824268 : p2 = simplify_shallow(p2);
1785 1824268 : if (!isexactzero(p2) || (--dx < 0)) break;
1786 : }
1787 1824268 : 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 18937328 : avy = avma;
1837 18937328 : dz = dx-dy;
1838 18937328 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1839 18937289 : x += 2;
1840 18937289 : z += 2;
1841 18937289 : y += 2;
1842 18937289 : gel(z,dz) = gcopy(p2);
1843 :
1844 51803562 : for (i=dx-1; i>=dy; i--)
1845 : {
1846 32866354 : av1=avma; p1=gel(x,i);
1847 1135066020 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1848 32861681 : if (y_lead) p1 = simplify(f(p1,y_lead));
1849 :
1850 32861682 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1851 : else
1852 22893541 : p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1853 32865835 : gel(z,i-dy) = p1;
1854 : }
1855 18937208 : if (!pr) return gerepileupto(av,z-2);
1856 :
1857 12148409 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1858 12509447 : for (sx=0; ; i--)
1859 : {
1860 12509447 : p1 = gel(x,i);
1861 : /* we always enter this loop at least once */
1862 31483655 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1863 12508354 : if (mod && avma==av1) p1 = gmul(p1,mod);
1864 12508354 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1865 3571938 : if (!isexactzero(p1)) break;
1866 3428297 : if (!i) break;
1867 361034 : set_avma(av1);
1868 : }
1869 12147710 : if (pr == ONLY_DIVIDES)
1870 : {
1871 2709 : if (sx) return gc_NULL(av);
1872 2688 : set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1873 : }
1874 12145001 : lr=i+3; rem -= lr;
1875 12145001 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1876 12101382 : else p1 = gerepileupto((pari_sp)rem,p1);
1877 12145615 : rem[0] = evaltyp(t_POL) | evallg(lr);
1878 12145584 : rem[1] = z[-1];
1879 12145584 : rem += 2;
1880 12145584 : gel(rem,i) = p1;
1881 16676266 : for (i--; i>=0; i--)
1882 : {
1883 4530696 : av1=avma; p1 = gel(x,i);
1884 12890097 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1885 4530235 : if (mod && avma==av1) p1 = gmul(p1,mod);
1886 4530480 : gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1887 : }
1888 12145570 : rem -= 2;
1889 12145570 : if (!sx) (void)normalizepol_lg(rem, lr);
1890 12145764 : if (pr == ONLY_REM) return gerepileupto(av,rem);
1891 6633472 : z -= 2;
1892 : {
1893 6633472 : GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1894 6633472 : gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1895 : }
1896 : }
1897 :
1898 : GEN
1899 14043693 : RgX_divrem(GEN x, GEN y, GEN *pr)
1900 : {
1901 14043693 : if (pr == ONLY_REM) return RgX_rem(x, y);
1902 14043693 : 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 136541 : 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 136541 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1914 136541 : vx = varn(x);
1915 136541 : dx = degpol(x);
1916 136541 : dy = degpol(y);
1917 136541 : if (dx < dy)
1918 : {
1919 28316 : if (pr)
1920 : {
1921 28288 : av0 = avma; x = RgXQX_red(x, T);
1922 28288 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1923 28288 : if (pr == ONLY_REM) return x;
1924 0 : *pr = x;
1925 : }
1926 28 : return pol_0(vx);
1927 : }
1928 108225 : lead = leading_coeff(y);
1929 108225 : if (!dy) /* y is constant */
1930 : {
1931 602 : if (pr && pr != ONLY_DIVIDES)
1932 : {
1933 0 : if (pr == ONLY_REM) return pol_0(vx);
1934 0 : *pr = pol_0(vx);
1935 : }
1936 602 : 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 107623 : av0 = avma; dz = dx-dy;
1941 107623 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
1942 107623 : set_avma(av0);
1943 107623 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1944 107623 : x += 2; y += 2; z += 2;
1945 :
1946 107623 : p1 = gel(x,dx); av = avma;
1947 107623 : gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
1948 556143 : for (i=dx-1; i>=dy; i--)
1949 : {
1950 448520 : av=avma; p1=gel(x,i);
1951 2352606 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1952 448501 : if (lead) p1 = gmul(grem(p1, T), lead);
1953 448504 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
1954 : }
1955 107623 : if (!pr) { guncloneNULL(lead); return z-2; }
1956 :
1957 106349 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
1958 179408 : for (sx=0; ; i--)
1959 : {
1960 179408 : p1 = gel(x,i);
1961 618427 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1962 179409 : tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
1963 100135 : if (!i) break;
1964 73059 : set_avma(av);
1965 : }
1966 106349 : if (pr == ONLY_DIVIDES)
1967 : {
1968 27034 : guncloneNULL(lead);
1969 27034 : if (sx) return gc_NULL(av0);
1970 24528 : return gc_const((pari_sp)rem, z-2);
1971 : }
1972 79315 : lr=i+3; rem -= lr;
1973 79315 : rem[0] = evaltyp(t_POL) | evallg(lr);
1974 79315 : rem[1] = z[-1];
1975 79315 : p1 = gerepile((pari_sp)rem,tetpil,p1);
1976 79315 : rem += 2; gel(rem,i) = p1;
1977 169687 : for (i--; i>=0; i--)
1978 : {
1979 90372 : av=avma; p1 = gel(x,i);
1980 306377 : for (j=0; j<=i && j<=dz; j++)
1981 216005 : p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1982 90372 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
1983 : }
1984 79315 : rem -= 2;
1985 79315 : guncloneNULL(lead);
1986 79315 : if (!sx) (void)normalizepol_lg(rem, lr);
1987 79315 : 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 1318047 : rem(GEN c, GEN T)
1998 : {
1999 1318047 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2000 1318048 : return c;
2001 : }
2002 :
2003 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2004 : int
2005 17083 : ZXQX_dvd(GEN x, GEN y, GEN T)
2006 : {
2007 : long dx, dy, i, T_ismonic;
2008 17083 : pari_sp av = avma, av2;
2009 : GEN y_lead;
2010 :
2011 17083 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2012 17083 : dy = degpol(y); y_lead = gel(y,dy+2);
2013 17083 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2014 : /* if monic, no point in using pseudo-division */
2015 17083 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2016 14570 : T_ismonic = gequal1(leading_coeff(T));
2017 14570 : dx = degpol(x);
2018 14570 : if (dx < dy) return !signe(x);
2019 14570 : (void)new_chunk(2);
2020 14570 : x = RgX_recip_i(x)+2;
2021 14570 : y = RgX_recip_i(y)+2;
2022 : /* pay attention to sparse divisors */
2023 29641 : for (i = 1; i <= dy; i++)
2024 15071 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2025 14570 : av2 = avma;
2026 : for (;;)
2027 72617 : {
2028 87187 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2029 87187 : x0 = gneg(x0);
2030 87187 : m = gcdii(cx, y0);
2031 87187 : if (!equali1(m))
2032 : {
2033 84692 : x0 = gdiv(x0, m);
2034 84692 : y0 = diviiexact(y0, m);
2035 84692 : if (equali1(y0)) y0 = NULL;
2036 : }
2037 178147 : for (i=1; i<=dy; i++)
2038 : {
2039 90960 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2040 90960 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2041 90960 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2042 90960 : gel(x,i) = c;
2043 : }
2044 686726 : for ( ; i<=dx; i++)
2045 : {
2046 599539 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2047 599539 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2048 599539 : gel(x,i) = c;
2049 : }
2050 101523 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2051 87187 : if (dx < dy) break;
2052 72617 : 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 14570 : return gc_bool(av, dx < 0);
2059 : }
2060 :
2061 : /* T either NULL or a t_POL. */
2062 : GEN
2063 138459 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2064 : {
2065 138459 : long vx = varn(x), dx, dy, dz, i, lx, p;
2066 138459 : pari_sp av = avma, av2;
2067 : GEN y_lead;
2068 :
2069 138459 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2070 138459 : dy = degpol(y); y_lead = gel(y,dy+2);
2071 : /* if monic, no point in using pseudo-division */
2072 138459 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2073 131518 : dx = degpol(x);
2074 131518 : if (dx < dy) return RgX_copy(x);
2075 131518 : (void)new_chunk(2);
2076 131518 : x = RgX_recip_i(x)+2;
2077 131519 : y = RgX_recip_i(y)+2;
2078 : /* pay attention to sparse divisors */
2079 376624 : for (i = 1; i <= dy; i++)
2080 245106 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2081 131518 : dz = dx-dy; p = dz+1;
2082 131518 : av2 = avma;
2083 : for (;;)
2084 : {
2085 362128 : gel(x,0) = gneg(gel(x,0)); p--;
2086 1075744 : for (i=1; i<=dy; i++)
2087 : {
2088 713618 : GEN c = gmul(y_lead, gel(x,i));
2089 713573 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2090 713597 : gel(x,i) = rem(c, T);
2091 : }
2092 843279 : for ( ; i<=dx; i++)
2093 : {
2094 481153 : GEN c = gmul(y_lead, gel(x,i));
2095 481133 : gel(x,i) = rem(c, T);
2096 : }
2097 391051 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2098 362127 : if (dx < dy) break;
2099 230610 : 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 131517 : if (dx < 0) return pol_0(vx);
2106 130334 : lx = dx+3; x -= 2;
2107 130334 : x[0] = evaltyp(t_POL) | evallg(lx);
2108 130334 : x[1] = evalsigne(1) | evalvarn(vx);
2109 130334 : x = RgX_recip_i(x);
2110 130334 : if (p)
2111 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2112 13432 : GEN t = y_lead;
2113 13432 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2114 0 : t = RgXQ_powu(t, p, T);
2115 : else
2116 13432 : t = gpowgs(t, p);
2117 35917 : for (i=2; i<lx; i++)
2118 : {
2119 22485 : GEN c = gmul(gel(x,i), t);
2120 22485 : gel(x,i) = rem(c,T);
2121 : }
2122 13432 : if (!T) return gerepileupto(av, x);
2123 : }
2124 116902 : return gerepilecopy(av, x);
2125 : }
2126 :
2127 : GEN
2128 138459 : 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 11888 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2133 : {
2134 11888 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2135 11888 : pari_sp av = avma, av2;
2136 : GEN z, r, ypow, y_lead;
2137 :
2138 11888 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2139 11888 : dy = degpol(y); y_lead = gel(y,dy+2);
2140 11888 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2141 7526 : dx = degpol(x);
2142 7526 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2143 7526 : 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 7428 : (void)new_chunk(2);
2152 7428 : x = RgX_recip_i(x)+2;
2153 7428 : y = RgX_recip_i(y)+2;
2154 : /* pay attention to sparse divisors */
2155 38664 : for (i = 1; i <= dy; i++)
2156 31236 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2157 7428 : dz = dx-dy; p = dz+1;
2158 7428 : lz = dz+3;
2159 7428 : z = cgetg(lz, t_POL);
2160 7428 : z[1] = evalsigne(1) | evalvarn(vx);
2161 28521 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2162 7428 : ypow = new_chunk(dz+1);
2163 7428 : gel(ypow,0) = gen_1;
2164 7428 : gel(ypow,1) = y_lead;
2165 13665 : for (i=2; i<=dz; i++)
2166 : {
2167 6237 : GEN c = gmul(gel(ypow,i-1), y_lead);
2168 6237 : gel(ypow,i) = rem(c,T);
2169 : }
2170 7428 : av2 = avma;
2171 7428 : for (iz=2;;)
2172 : {
2173 15549 : p--;
2174 15549 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2175 69229 : for (i=1; i<=dy; i++)
2176 : {
2177 53680 : GEN c = gmul(y_lead, gel(x,i));
2178 53680 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2179 53680 : gel(x,i) = rem(c, T);
2180 : }
2181 40932 : for ( ; i<=dx; i++)
2182 : {
2183 25383 : GEN c = gmul(y_lead, gel(x,i));
2184 25383 : gel(x,i) = rem(c,T);
2185 : }
2186 15549 : x++; dx--;
2187 21093 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2188 15549 : if (dx < dy) break;
2189 8121 : 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 13938 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2198 7428 : if (dx < 0)
2199 182 : x = pol_0(vx);
2200 : else
2201 : {
2202 7246 : lx = dx+3; x -= 2;
2203 7246 : x[0] = evaltyp(t_POL) | evallg(lx);
2204 7246 : x[1] = evalsigne(1) | evalvarn(vx);
2205 7246 : x = RgX_recip_i(x);
2206 : }
2207 7428 : z = RgX_recip_i(z);
2208 7428 : r = x;
2209 7428 : if (p)
2210 : {
2211 3311 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2212 3311 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2213 : }
2214 7428 : *ptr = r; return gc_all(av, 2, &z, ptr);
2215 : }
2216 : GEN
2217 11650 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2218 11650 : { 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 496829432 : 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 26145 : 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 603748 : _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 65138 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2255 : static GEN
2256 86452 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2257 : static GEN
2258 270579 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2259 : static GEN
2260 1004402 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2261 : static GEN
2262 978168 : _one(void *data) { return pol_1(varn((GEN)data)); }
2263 : static GEN
2264 840 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2265 : static GEN
2266 703243 : _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 400515 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2279 : {
2280 400515 : int use_sqr = 2*degpol(x) >= degpol(T);
2281 400515 : 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 2598322 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2381 : {
2382 2598322 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2383 2598322 : 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 2178703 : RgXn_div_gen(GEN g, GEN f, long e)
2396 : {
2397 : pari_sp av;
2398 : ulong mask;
2399 : GEN W, a;
2400 2178703 : long v = varn(f), n = 1;
2401 :
2402 2178703 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2403 2178696 : a = ginv(gel(f,2));
2404 2178695 : if (e == 1 && !g) return scalarpol(a, v);
2405 2174684 : else if (e == 2 && !g)
2406 : {
2407 : GEN b;
2408 835134 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2409 244643 : b = gneg(b);
2410 244646 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2411 244646 : return deg1pol(b, a, v);
2412 : }
2413 1339550 : av = avma;
2414 1339550 : W = scalarpol_shallow(a,v);
2415 1339550 : mask = quadratic_prec_mask(e);
2416 3932006 : while (mask > 1)
2417 : {
2418 : GEN u, fr;
2419 2592456 : long n2 = n;
2420 2592456 : n<<=1; if (mask & 1) n--;
2421 2592456 : mask >>= 1;
2422 2592456 : fr = RgXn_red_shallow(f, n);
2423 2592456 : if (mask>1 || !g)
2424 : {
2425 1313964 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2426 1313964 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2427 : }
2428 : else
2429 : {
2430 1278492 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2431 1278492 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2432 1278492 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2433 : }
2434 2592456 : 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 1339550 : 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 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2471 : {
2472 : GEN r;
2473 91 : if (lgefint(p) == 3)
2474 : {
2475 91 : ulong pp = uel(p, 2);
2476 91 : if (pp == 2)
2477 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2478 : else
2479 63 : 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 91 : 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 900299 : RgXn_inv_fast(GEN x, long e)
2499 : {
2500 : GEN p, pol;
2501 : long pa;
2502 900299 : long t = RgX_type(x,&p,&pol,&pa);
2503 900300 : switch(t)
2504 : {
2505 91 : 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 900211 : default: return NULL;
2509 : }
2510 : }
2511 :
2512 : static GEN
2513 1278611 : RgXn_div_fast(GEN x, GEN y, long e)
2514 : {
2515 : GEN p, pol;
2516 : long pa;
2517 1278611 : long t = RgX_type2(x,y,&p,&pol,&pa);
2518 1278611 : 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 1278492 : default: return NULL;
2524 : }
2525 : }
2526 : #undef code
2527 :
2528 : GEN
2529 1278611 : RgXn_div_i(GEN g, GEN f, long e)
2530 : {
2531 1278611 : GEN h = RgXn_div_fast(g, f, e);
2532 1278611 : if (h) return h;
2533 1278492 : return RgXn_div_gen(g, f, e);
2534 : }
2535 :
2536 : GEN
2537 583120 : RgXn_div(GEN g, GEN f, long e)
2538 : {
2539 583120 : pari_sp av = avma;
2540 583120 : return gerepileupto(av, RgXn_div_i(g, f, e));
2541 : }
2542 :
2543 : GEN
2544 900298 : RgXn_inv_i(GEN f, long e)
2545 : {
2546 900298 : GEN h = RgXn_inv_fast(f, e);
2547 900301 : if (h) return h;
2548 900210 : return RgXn_div_gen(NULL, f, e);
2549 : }
2550 :
2551 : GEN
2552 806976 : RgXn_inv(GEN f, long e)
2553 : {
2554 806976 : pari_sp av = avma;
2555 806976 : return gerepileupto(av, RgXn_inv_i(f, e));
2556 : }
2557 :
2558 : /* Compute intformal(x^n*S)/x^(n+1) */
2559 : static GEN
2560 58714 : RgX_integXn(GEN x, long n)
2561 : {
2562 58714 : long i, lx = lg(x);
2563 : GEN y;
2564 58714 : if (lx == 2) return RgX_copy(x);
2565 52680 : y = cgetg(lx, t_POL); y[1] = x[1];
2566 110477 : for (i=2; i<lx; i++)
2567 57797 : gel(y,i) = gdivgs(gel(x,i), n+i-1);
2568 52680 : return RgX_renormalize_lg(y, lx);;
2569 : }
2570 :
2571 : GEN
2572 52857 : RgXn_expint(GEN h, long e)
2573 : {
2574 52857 : pari_sp av = avma, av2;
2575 52857 : long v = varn(h), n;
2576 52857 : GEN f = pol_1(v), g;
2577 : ulong mask;
2578 :
2579 52857 : if (!signe(h)) return f;
2580 52857 : g = pol_1(v);
2581 52857 : n = 1; mask = quadratic_prec_mask(e);
2582 52857 : av2 = avma;
2583 58716 : for (;mask>1;)
2584 : {
2585 : GEN u, w;
2586 58716 : long n2 = n;
2587 58716 : n<<=1; if (mask & 1) n--;
2588 58716 : mask >>= 1;
2589 58716 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2590 58715 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2591 58714 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2592 58713 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2593 58714 : 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 52855 : 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 115507 : RgXQ_powu(GEN x, ulong n, GEN T)
2693 : {
2694 115507 : pari_sp av = avma;
2695 :
2696 115507 : if (!n) return pol_1(varn(x));
2697 67767 : if (n == 1) return RgX_copy(x);
2698 19670 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2699 19670 : return gerepilecopy(av, x);
2700 : }
2701 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2702 : GEN
2703 100165 : RgXQ_pow(GEN x, GEN n, GEN T)
2704 : {
2705 : pari_sp av;
2706 100165 : long s = signe(n);
2707 :
2708 100165 : if (!s) return pol_1(varn(x));
2709 100165 : if (is_pm1(n) == 1)
2710 86452 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2711 13713 : av = avma;
2712 13713 : if (s < 0) x = RgXQ_inv(x, T);
2713 13713 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2714 13713 : return gerepilecopy(av, x);
2715 : }
2716 : static GEN
2717 213715 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2718 : static GEN
2719 113840 : _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 10906 : ZXQ_powers(GEN x, long l, GEN T)
2724 : {
2725 10906 : int use_sqr = 2*degpol(x) >= degpol(T);
2726 10906 : 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 188545 : ZXQ_powu(GEN x, ulong n, GEN T)
2732 : {
2733 188545 : pari_sp av = avma;
2734 :
2735 188545 : if (!n) return pol_1(varn(x));
2736 188545 : if (n == 1) return ZX_copy(x);
2737 125035 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2738 125041 : return gerepilecopy(av, x);
2739 : }
2740 :
2741 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2742 : GEN
2743 4214 : RgXQ_powers(GEN x, long l, GEN T)
2744 : {
2745 4214 : int use_sqr = 2*degpol(x) >= degpol(T);
2746 4214 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2747 : }
2748 :
2749 : GEN
2750 7013 : RgXQV_factorback(GEN L, GEN e, GEN T)
2751 : {
2752 7013 : 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 8344 : QXQ_powers(GEN a, long n, GEN T)
2758 : {
2759 : GEN den, v;
2760 8344 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2761 8323 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2762 : /* den*a integral; v[i+1] = (den*a)^i in K */
2763 8323 : if (den)
2764 : { /* restore denominators */
2765 5222 : GEN d = den;
2766 : long i;
2767 5222 : gel(v,2) = a;
2768 26502 : for (i=3; i<=n+1; i++) {
2769 21280 : d = mulii(d,den);
2770 21280 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2771 : }
2772 : }
2773 8323 : return v;
2774 : }
2775 :
2776 : static GEN
2777 2492 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2778 : {
2779 2492 : long l, i, m = 0;
2780 : GEN dz, z;
2781 2492 : GEN V = cgetg_copy(v, &l);
2782 8701 : for (i = imin; i < l; i++)
2783 : {
2784 6209 : GEN c = gel(v, i);
2785 6209 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2786 : }
2787 2492 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2788 2723 : for (i = 1; i < imin; i++) V[i] = v[i];
2789 8701 : for (i = imin; i < l; i++)
2790 : {
2791 6209 : GEN c = gel(v,i);
2792 6209 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2793 6209 : gel(V,i) = c;
2794 : }
2795 2492 : return V;
2796 : }
2797 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2798 : GEN
2799 2261 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2800 2261 : { 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 5179 : RgXQ_norm(GEN x, GEN T)
2814 : {
2815 : pari_sp av;
2816 5179 : long dx = degpol(x);
2817 : GEN L, y;
2818 5179 : if (degpol(T)==0) return gpowgs(x,0);
2819 5172 : av = avma; y = resultant(T, x);
2820 5172 : L = leading_coeff(T);
2821 5172 : if (gequal1(L) || !signe(x)) return y;
2822 0 : return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2823 : }
2824 :
2825 : GEN
2826 455 : RgXQ_trace(GEN x, GEN T)
2827 : {
2828 455 : pari_sp av = avma;
2829 : GEN dT, z;
2830 : long n;
2831 455 : if (degpol(T)==0) return gmulgs(x,0);
2832 448 : dT = RgX_deriv(T); n = degpol(dT);
2833 448 : z = RgXQ_mul(x, dT, T);
2834 448 : if (degpol(z)<n) return gc_const(av, gen_0);
2835 399 : return gerepileupto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2836 : }
2837 :
2838 : GEN
2839 2771440 : RgX_blocks(GEN P, long n, long m)
2840 : {
2841 2771440 : GEN z = cgetg(m+1,t_VEC);
2842 2771440 : long i,j, k=2, l = lg(P);
2843 8507898 : for(i=1; i<=m; i++)
2844 : {
2845 5736458 : GEN zi = cgetg(n+2,t_POL);
2846 5736458 : zi[1] = P[1];
2847 5736458 : gel(z,i) = zi;
2848 16905856 : for(j=2; j<n+2; j++)
2849 11169398 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2850 5736458 : zi = RgX_renormalize_lg(zi, n+2);
2851 : }
2852 2771440 : return z;
2853 : }
2854 :
2855 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2856 : void
2857 47587561 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
2858 : {
2859 47587561 : long n = degpol(p), v = varn(p), n0, n1, i;
2860 : GEN p0, p1;
2861 :
2862 47587440 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2863 :
2864 47494491 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2865 47494491 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2866 47496490 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2867 140951297 : for (i=0; i<n1; i++)
2868 : {
2869 93450674 : p0[2+i] = p[2+(i<<1)];
2870 93450674 : p1[2+i] = p[3+(i<<1)];
2871 : }
2872 47500623 : if (n1 != n0)
2873 17416030 : p0[2+i] = p[2+(i<<1)];
2874 47500623 : *pe = normalizepol(p0);
2875 47497867 : *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 43547 : RgX_splitting(GEN p, long k)
2881 : {
2882 43547 : long n = degpol(p), v = varn(p), m, i, j, l;
2883 : GEN r;
2884 :
2885 43547 : m = n/k;
2886 43547 : r = cgetg(k+1,t_VEC);
2887 233807 : for(i=1; i<=k; i++)
2888 : {
2889 190260 : gel(r,i) = cgetg(m+3, t_POL);
2890 190260 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
2891 : }
2892 571487 : for (j=1, i=0, l=2; i<=n; i++)
2893 : {
2894 527940 : gmael(r,j,l) = gel(p,2+i);
2895 527940 : if (j==k) { j=1; l++; } else j++;
2896 : }
2897 233807 : for(i=1; i<=k; i++)
2898 190260 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2899 43547 : 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 178559 : Kronecker_to_mod(GEN z, GEN T)
2913 : {
2914 178559 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2915 178559 : GEN x, t = cgetg(N,t_POL);
2916 178559 : t[1] = T[1];
2917 178559 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2918 178559 : x[1] = z[1];
2919 178559 : T = RgX_copy(T);
2920 1318040 : for (i=2; i<lx+2; i++, z+= N-2)
2921 : {
2922 5462080 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2923 1139481 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2924 : }
2925 178559 : N = (l-2) % (N-2) + 2;
2926 459076 : for (j=2; j<N; j++) t[j] = z[j];
2927 178559 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2928 178559 : return normalizepol_lg(x, i+1);
2929 : }
2930 :
2931 : /*******************************************************************/
2932 : /* */
2933 : /* Domain detection */
2934 : /* */
2935 : /*******************************************************************/
2936 :
2937 : static GEN
2938 152783 : zero_FpX_mod(GEN p, long v)
2939 : {
2940 152783 : GEN r = cgetg(3,t_POL);
2941 152783 : r[1] = evalvarn(v);
2942 152783 : gel(r,2) = mkintmod(gen_0, icopy(p));
2943 152783 : return r;
2944 : }
2945 :
2946 : static GEN
2947 818061 : RgX_mul_FpX(GEN x, GEN y, GEN p)
2948 : {
2949 818061 : pari_sp av = avma;
2950 : GEN r;
2951 818061 : if (lgefint(p) == 3)
2952 : {
2953 766705 : ulong pp = uel(p, 2);
2954 766705 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2955 : RgX_to_Flx(y, pp), pp));
2956 : }
2957 : else
2958 51356 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2959 818061 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2960 703001 : 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 416946 : RgX_liftred(GEN x, GEN T)
2990 416946 : { return RgXQX_red(liftpol_shallow(x), T); }
2991 :
2992 : static GEN
2993 163468 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
2994 : {
2995 163468 : pari_sp av = avma;
2996 163468 : long dT = degpol(T);
2997 163468 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
2998 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
2999 163468 : 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 1239 : 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 13390 : RgX_sqr_QXQX(GEN x, GEN T)
3034 : {
3035 13390 : pari_sp av = avma;
3036 13390 : long dT = degpol(T);
3037 13390 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3038 13390 : return gerepileupto(av, Kronecker_to_mod(r, T));
3039 : }
3040 :
3041 : static GEN
3042 68670 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3043 : {
3044 68670 : pari_sp av = avma;
3045 : GEN r;
3046 68670 : if (lgefint(p) == 3)
3047 : {
3048 51880 : ulong pp = uel(p, 2);
3049 51880 : 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 68670 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3055 31073 : return gerepileupto(av, FpX_to_mod(r, p));
3056 : }
3057 :
3058 : static GEN
3059 38310 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3060 : {
3061 38310 : pari_sp av = avma;
3062 : GEN r;
3063 38310 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3064 38310 : 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 74477847 : RgX_mul_fast(GEN x, GEN y)
3089 : {
3090 : GEN p, pol;
3091 : long pa;
3092 74477847 : long t = RgX_type2(x,y, &p,&pol,&pa);
3093 74477980 : switch(t)
3094 : {
3095 38827319 : case t_INT: return ZX_mul(x,y);
3096 3279171 : case t_FRAC: return QX_mul(x,y);
3097 104742 : case t_FFELT: return FFX_mul(x, y, pol);
3098 818061 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3099 163468 : case code(t_POLMOD, t_INT):
3100 : case code(t_POLMOD, t_FRAC):
3101 163468 : return RgX_mul_QXQX(x, y, pol);
3102 2028 : case code(t_POLMOD, t_INTMOD):
3103 2028 : return RgX_mul_FpXQX(x, y, pol, p);
3104 31283191 : default: return NULL;
3105 : }
3106 : }
3107 : static GEN
3108 1206302 : RgX_sqr_fast(GEN x)
3109 : {
3110 : GEN p, pol;
3111 : long pa;
3112 1206302 : long t = RgX_type(x,&p,&pol,&pa);
3113 1206302 : switch(t)
3114 : {
3115 1149295 : case t_INT: return ZX_sqr(x);
3116 30709 : case t_FRAC: return QX_sqr(x);
3117 2534 : case t_FFELT: return FFX_sqr(x, pol);
3118 1365 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3119 13390 : case code(t_POLMOD, t_INT):
3120 : case code(t_POLMOD, t_FRAC):
3121 13390 : return RgX_sqr_QXQX(x, pol);
3122 196 : case code(t_POLMOD, t_INTMOD):
3123 196 : return RgX_sqr_FpXQX(x, pol, p);
3124 8813 : default: return NULL;
3125 : }
3126 : }
3127 :
3128 : static GEN
3129 13858621 : RgX_rem_fast(GEN x, GEN y)
3130 : {
3131 : GEN p, pol;
3132 : long pa;
3133 13858621 : long t = RgX_type2(x,y, &p,&pol,&pa);
3134 13858812 : switch(t)
3135 : {
3136 3244455 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3137 1688320 : 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 68670 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3140 38310 : case code(t_POLMOD, t_INT):
3141 : case code(t_POLMOD, t_FRAC):
3142 38310 : return RgX_rem_QXQX(x, y, pol);
3143 72 : case code(t_POLMOD, t_INTMOD):
3144 72 : return RgX_rem_FpXQX(x, y, pol, p);
3145 8818901 : default: return NULL;
3146 : }
3147 : }
3148 :
3149 : #undef code
3150 :
3151 : GEN
3152 60825964 : RgX_mul(GEN x, GEN y)
3153 : {
3154 60825964 : GEN z = RgX_mul_fast(x,y);
3155 60825948 : if (!z) z = RgX_mul_i(x,y);
3156 60825477 : return z;
3157 : }
3158 :
3159 : GEN
3160 1194150 : RgX_sqr(GEN x)
3161 : {
3162 1194150 : GEN z = RgX_sqr_fast(x);
3163 1194143 : if (!z) z = RgX_sqr_i(x);
3164 1194143 : return z;
3165 : }
3166 :
3167 : GEN
3168 13858635 : RgX_rem(GEN x, GEN y)
3169 : {
3170 13858635 : GEN z = RgX_rem_fast(x, y);
3171 13858776 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3172 13858667 : 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 10994840 : RgXn_mul(GEN f, GEN g, long n)
3185 : {
3186 10994840 : pari_sp av = avma;
3187 10994840 : GEN h = RgX_mul_fast(f,g);
3188 10994841 : if (!h) return RgXn_mul2(f,g,n);
3189 1774019 : if (degpol(h) < n) return h;
3190 375838 : 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 2657051 : RgX_mulhigh_i(GEN f, GEN g, long n)
3206 : {
3207 2657051 : GEN h = RgX_mul_fast(f,g);
3208 2657050 : 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 : }
|
