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 19609657 : brent_kung_optpow(long d, long n, long m)
32 : {
33 : long p, r;
34 19609657 : long pold=1, rold=n*(d-1);
35 91826747 : for(p=2; p<=d; p++)
36 : {
37 72217090 : r = m*(p-1) + n*((d-1)/p);
38 72217090 : if (r<rold) { pold=p; rold=r; }
39 : }
40 19609657 : return pold;
41 : }
42 :
43 : static GEN
44 14037253 : 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 14037253 : pari_sp av = avma;
48 : long i;
49 14037253 : GEN z = cmul(E,P,a,ff->one(E));
50 14013303 : if (!z) z = gen_0;
51 74523397 : for (i=1; i<=n; i++)
52 : {
53 60510289 : GEN t = cmul(E,P,a+i,gel(V,i+1));
54 60514951 : if (t) {
55 44602380 : z = ff->add(E, z, t);
56 44567690 : if (gc_needed(av,2)) z = gc_upto(av, z);
57 : }
58 : }
59 14013108 : 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 12120596 : 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 12120596 : pari_sp av = avma;
72 12120596 : long l = lg(V)-1;
73 : GEN z, u;
74 :
75 12120596 : if (d < 0) return ff->zero(E);
76 11535053 : if (d < l) return gc_upto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
77 1005073 : if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
78 1005073 : 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 1005073 : d -= l;
84 1005073 : z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
85 2501241 : while (d >= l-1)
86 : {
87 1494898 : d -= l-1;
88 1494898 : u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
89 1495193 : z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
90 1494848 : if (gc_needed(av,2))
91 91 : z = gc_upto(av, z);
92 : }
93 1006343 : u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
94 1006354 : z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
95 1006341 : return gc_upto(av, ff->red(E,z));
96 : }
97 :
98 : GEN
99 881599 : 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 881599 : pari_sp av = avma;
103 : GEN z, V;
104 : long rtd;
105 881599 : if (d < 0) return ff->zero(E);
106 880325 : rtd = (long) sqrt((double)d);
107 880325 : V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
108 880338 : z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
109 880331 : return gc_upto(av, z);
110 : }
111 :
112 : static GEN
113 2150489 : _gen_nored(void *E, GEN x) { (void)E; return x; }
114 : static GEN
115 19637444 : _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 1974339 : _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
120 : static GEN
121 635449 : _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
122 : static GEN
123 2193292 : _gen_one(void *E) { (void)E; return gen_1; }
124 : static GEN
125 301 : _gen_zero(void *E) { (void)E; return gen_0; }
126 :
127 : static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
128 : _gen_mul, _gen_sqr,_gen_one,_gen_zero };
129 :
130 : static GEN
131 512818 : _gen_cmul(void *E, GEN P, long a, GEN x)
132 512818 : {(void)E; return gmul(gel(P,a+2), x);}
133 :
134 : GEN
135 142885 : RgX_RgV_eval(GEN Q, GEN x)
136 : {
137 142885 : return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
138 : }
139 :
140 : GEN
141 0 : RgX_Rg_eval_bk(GEN Q, GEN x)
142 : {
143 0 : return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
144 : }
145 :
146 : GEN
147 2947 : RgXV_RgV_eval(GEN Q, GEN x)
148 : {
149 2947 : long i, l = lg(Q), vQ = gvar(Q);
150 2947 : GEN v = cgetg(l, t_VEC);
151 248311 : for (i = 1; i < l; i++)
152 : {
153 245364 : GEN Qi = gel(Q, i);
154 245364 : gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
155 : }
156 2947 : return v;
157 : }
158 :
159 : GEN
160 629213 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B, long d)
161 : {
162 629213 : pari_sp av = avma, btop;
163 629213 : long i, o, dP = degpol(P);
164 : GEN s;
165 629213 : if (signe(P)==0) return gen_0;
166 504754 : s = dP >=d ? gel(P, d+2): gen_0;
167 504754 : if (d == 0) return gcopy(s);
168 504033 : o = RgX_deflate_order(P); if (d % o) o = 1;
169 504041 : if (o > 1) A = gpowgs(A, o);
170 504045 : btop = avma;
171 2015531 : for (i = d-o; i >= 0; i-=o)
172 1511487 : if (i <= dP)
173 : {
174 1502191 : s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
175 1502190 : if (gc_needed(btop,1))
176 : {
177 13 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
178 13 : s = gc_upto(btop, s);
179 : }
180 : }
181 504044 : return gc_upto(av, s);
182 : }
183 :
184 : GEN
185 16604 : RgX_homogenous_eval(GEN P, GEN A, GEN B, long d)
186 : {
187 16604 : pari_sp av = avma;
188 16604 : return gc_upto(av, RgX_homogenous_evalpow(P, A, gpowers(B,d), d));
189 : }
190 :
191 : /*FIXME : Actually QXQX_QXQX_QXQXV_homogenous_evalpow */
192 : GEN
193 1652 : QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
194 : {
195 1652 : pari_sp av = avma;
196 1652 : long i, d = degpol(P), v = varn(A);
197 : GEN s;
198 1652 : if (signe(P)==0) return pol_0(v);
199 1652 : if (d == 0) return scalarpol(gel(P, d+2), v);
200 1232 : s = scalarpol_shallow(gel(P, d+2), v);
201 4963 : for (i = d-1; i >= 0; i--)
202 : {
203 3731 : GEN c = gel(P,i+2), b = gel(B,d+1-i);
204 3731 : s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
205 3731 : if (gc_needed(av,1))
206 : {
207 0 : if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
208 0 : s = gc_upto(av, s);
209 : }
210 : }
211 1232 : return gc_upto(av, s);
212 : }
213 :
214 : const struct bb_algebra *
215 294988 : get_Rg_algebra(void)
216 : {
217 294988 : return &Rg_algebra;
218 : }
219 :
220 : static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
221 :
222 : static GEN
223 11809 : _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
224 : {
225 : (void) E;
226 11809 : return RgX_divrem(x, y, r);
227 : }
228 :
229 : GEN
230 3234 : RgX_digits(GEN x, GEN T)
231 : {
232 3234 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
233 3234 : if (signe(x)==0) return(cgetg(1, t_VEC));
234 3234 : return gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
235 : }
236 :
237 : /*******************************************************************/
238 : /* */
239 : /* RgX */
240 : /* */
241 : /*******************************************************************/
242 :
243 : long
244 24796437 : RgX_equal(GEN x, GEN y)
245 : {
246 24796437 : long i = lg(x);
247 :
248 24796437 : if (i != lg(y)) return 0;
249 106702741 : for (i--; i > 1; i--)
250 82234298 : if (!gequal(gel(x,i),gel(y,i))) return 0;
251 24468443 : return 1;
252 : }
253 :
254 : /* Returns 1 in the base ring over which x is defined */
255 : /* HACK: this also works for t_SER */
256 : GEN
257 156849065 : Rg_get_1(GEN x)
258 : {
259 : GEN p, T;
260 156849065 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
261 156849066 : if (RgX_type_is_composite(tx))
262 11194785 : RgX_type_decode(tx, &i /*junk*/, &tx);
263 156849066 : switch(tx)
264 : {
265 422919 : case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
266 1015 : case t_PADIC: return cvtop(gen_1, p, lx);
267 7077 : case t_FFELT: return FF_1(T);
268 156418055 : default: return gen_1;
269 : }
270 : }
271 : /* Returns 0 in the base ring over which x is defined */
272 : /* HACK: this also works for t_SER */
273 : GEN
274 1862498 : Rg_get_0(GEN x)
275 : {
276 : GEN p, T;
277 1862498 : long i, lx, tx = Rg_type(x, &p, &T, &lx);
278 1862498 : if (RgX_type_is_composite(tx))
279 25053 : RgX_type_decode(tx, &i /*junk*/, &tx);
280 1862498 : switch(tx)
281 : {
282 504 : case t_INTMOD: retmkintmod(gen_0, icopy(p));
283 42 : case t_PADIC: return zeropadic(p, lx);
284 210 : case t_FFELT: return FF_zero(T);
285 1861742 : default: return gen_0;
286 : }
287 : }
288 :
289 : GEN
290 7455 : QX_ZXQV_eval(GEN P, GEN V, GEN dV)
291 : {
292 7455 : long i, n = degpol(P);
293 : GEN z, dz, dP;
294 7455 : if (n < 0) return gen_0;
295 7455 : P = Q_remove_denom(P, &dP);
296 7455 : z = gel(P,2); if (n == 0) return icopy(z);
297 4214 : if (dV) z = mulii(dV, z); /* V[1] = dV */
298 4214 : z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
299 7931 : for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
300 4214 : dz = mul_denom(dP, dV);
301 4214 : return dz? RgX_Rg_div(z, dz): z;
302 : }
303 :
304 : /* Return P(h * x), not memory clean */
305 : GEN
306 52062 : RgX_unscale(GEN P, GEN h)
307 : {
308 52062 : long i, l = lg(P);
309 52062 : GEN hi = gen_1, Q = cgetg(l, t_POL);
310 52062 : Q[1] = P[1];
311 52062 : if (l == 2) return Q;
312 40820 : gel(Q,2) = gcopy(gel(P,2));
313 95249 : for (i=3; i<l; i++)
314 : {
315 54431 : hi = gmul(hi,h);
316 54430 : gel(Q,i) = gmul(gel(P,i), hi);
317 : }
318 40818 : return Q;
319 : }
320 : /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 and 2^k */
321 : GEN
322 1179248 : ZX_z_unscale(GEN P, long h)
323 : {
324 1179248 : long i, l = lg(P);
325 1179248 : GEN Q = cgetg(l, t_POL);
326 1179248 : Q[1] = P[1];
327 1179248 : if (l == 2) return Q;
328 1136303 : gel(Q,2) = gel(P,2);
329 1136303 : if (l == 3) return Q;
330 1115709 : if (h == -1)
331 242040 : for (i = 3; i < l; i++)
332 : {
333 199257 : gel(Q,i) = negi(gel(P,i));
334 199257 : if (++i == l) break;
335 148226 : gel(Q,i) = gel(P,i);
336 : }
337 1021895 : else if (h > 0 && !(h & (h-1))) return ZX_unscale2n(P, vals(h));
338 : else
339 : {
340 : GEN hi;
341 493193 : gel(Q,3) = mulis(gel(P,3), h);
342 493193 : hi = sqrs(h);
343 2458897 : for (i = 4; i < l; i++)
344 : {
345 1965703 : gel(Q,i) = mulii(gel(P,i), hi);
346 1965704 : if (i != l-1) hi = mulis(hi,h);
347 : }
348 : }
349 587008 : return Q;
350 : }
351 : /* P a ZX, h a t_INT. Return P(h * x), not memory clean */
352 : GEN
353 681402 : ZX_unscale(GEN P, GEN h)
354 : {
355 : long i, l;
356 : GEN Q, hi;
357 681402 : i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
358 888 : l = lg(P); Q = cgetg(l, t_POL);
359 888 : Q[1] = P[1];
360 888 : if (l == 2) return Q;
361 888 : gel(Q,2) = gel(P,2);
362 888 : if (l == 3) return Q;
363 888 : hi = h;
364 888 : gel(Q,3) = mulii(gel(P,3), hi);
365 2783 : for (i = 4; i < l; i++)
366 : {
367 1895 : hi = mulii(hi,h);
368 1895 : gel(Q,i) = mulii(gel(P,i), hi);
369 : }
370 888 : return Q;
371 : }
372 : /* P a ZX. Return P(x << n), not memory clean */
373 : GEN
374 1988707 : ZX_unscale2n(GEN P, long n)
375 : {
376 1988707 : long i, ni = n, l = lg(P);
377 1988707 : GEN Q = cgetg(l, t_POL);
378 1988704 : Q[1] = P[1];
379 1988704 : if (l == 2) return Q;
380 1942907 : gel(Q,2) = gel(P,2);
381 1942907 : if (l == 3) return Q;
382 1936852 : gel(Q,3) = shifti(gel(P,3), ni);
383 6775079 : for (i=4; i<l; i++)
384 : {
385 4838273 : ni += n;
386 4838273 : gel(Q,i) = shifti(gel(P,i), ni);
387 : }
388 1936806 : return Q;
389 : }
390 : /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
391 : GEN
392 12997 : ZX_unscale_div(GEN P, GEN h)
393 : {
394 12997 : long i, l = lg(P);
395 12997 : GEN hi, Q = cgetg(l, t_POL);
396 12997 : Q[1] = P[1];
397 12997 : if (l == 2) return Q;
398 12997 : gel(Q,2) = diviiexact(gel(P,2), h);
399 12997 : if (l == 3) return Q;
400 12997 : gel(Q,3) = gel(P,3);
401 12997 : if (l == 4) return Q;
402 12997 : hi = h;
403 12997 : gel(Q,4) = mulii(gel(P,4), hi);
404 64121 : for (i=5; i<l; i++)
405 : {
406 51124 : hi = mulii(hi,h);
407 51124 : gel(Q,i) = mulii(gel(P,i), hi);
408 : }
409 12997 : return Q;
410 : }
411 : /* P(h*X) / h^k, assuming the result is a ZX */
412 : GEN
413 1393 : ZX_unscale_divpow(GEN P, GEN h, long k)
414 : {
415 1393 : long i, j, l = lg(P);
416 1393 : GEN H, Q = cgetg(l, t_POL);
417 1393 : Q[1] = P[1]; if (l == 2) return Q;
418 1393 : H = gpowers(h, maxss(k, l - 3 - k));
419 5572 : for (i = 2, j = k+1; j > 1 && i < l; i++)
420 4179 : gel(Q, i) = diviiexact(gel(P, i), gel(H, j--));
421 1393 : if (i == l) return Q;
422 1393 : gel(Q, i) = gel(P, i); i++;
423 5082 : for (j = 2; i < l; i++) gel(Q, i) = mulii(gel(P, i), gel(H, j++));
424 1393 : return Q;
425 : }
426 :
427 : GEN
428 6489 : RgXV_unscale(GEN x, GEN h)
429 : {
430 6489 : if (isint1(h)) return gcopy(x);
431 18423 : pari_APPLY_same(RgX_unscale(gel(x,i), h));
432 : }
433 :
434 : /* Return h^degpol(P) P(x / h), not memory clean */
435 : GEN
436 4352770 : RgX_rescale(GEN P, GEN h)
437 : {
438 4352770 : long i, l = lg(P);
439 4352770 : GEN Q = cgetg(l,t_POL), hi = h;
440 4352760 : gel(Q,l-1) = gel(P,l-1);
441 11455753 : for (i=l-2; i>=2; i--)
442 : {
443 11453125 : gel(Q,i) = gmul(gel(P,i), hi);
444 11453039 : if (i == 2) break;
445 7102706 : hi = gmul(hi,h);
446 : }
447 4352961 : Q[1] = P[1]; return Q;
448 : }
449 :
450 : GEN
451 2401 : RgXV_rescale(GEN x, GEN h)
452 : {
453 2401 : if (isint1(h)) return RgX_copy(x);
454 16086 : pari_APPLY_same(RgX_rescale(gel(x,i), h));
455 : }
456 :
457 : /* A(X^d) --> A(X) */
458 : GEN
459 1187346 : RgX_deflate(GEN x0, long d)
460 : {
461 : GEN z, y, x;
462 1187346 : long i,id, dy, dx = degpol(x0);
463 1187342 : if (d == 1 || dx <= 0) return x0;
464 463421 : dy = dx/d;
465 463421 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
466 463416 : z = y + 2;
467 463416 : x = x0+ 2;
468 1775995 : for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
469 463416 : return y;
470 : }
471 :
472 : GEN
473 1260 : RgX_homogenize_deg(GEN P, long d, long v)
474 : {
475 : long i, l;
476 1260 : GEN Q = cgetg_copy(P, &l);
477 1260 : Q[1] = P[1];
478 4018 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
479 1260 : return Q;
480 : }
481 :
482 : GEN
483 18102 : RgX_homogenize(GEN P, long v)
484 : {
485 : long i, l, d;
486 18102 : GEN Q = cgetg_copy(P, &l);
487 18102 : Q[1] = P[1]; d = l-3;
488 159348 : for (i = 2; i < l; i++) gel(Q,i) = monomial(gel(P,i), d--, v);
489 18102 : return Q;
490 : }
491 :
492 : /* F a t_RFRAC */
493 : long
494 140 : rfrac_deflate_order(GEN F)
495 : {
496 140 : GEN N = gel(F,1), D = gel(F,2);
497 140 : long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
498 140 : if (m == 1) return 1;
499 49 : if (typ(N) == t_POL && varn(N) == varn(D))
500 28 : m = cgcd(m, RgX_deflate_order(N));
501 49 : return m;
502 : }
503 : /* F a t_RFRAC */
504 : GEN
505 140 : rfrac_deflate_max(GEN F, long *m)
506 : {
507 140 : *m = rfrac_deflate_order(F);
508 140 : return rfrac_deflate(F, *m);
509 : }
510 : /* F a t_RFRAC */
511 : GEN
512 140 : rfrac_deflate(GEN F, long m)
513 : {
514 140 : GEN N = gel(F,1), D = gel(F,2);
515 140 : if (m == 1) return F;
516 49 : if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
517 49 : D = RgX_deflate(D, m); return mkrfrac(N, D);
518 : }
519 :
520 : /* return x0(X^d) */
521 : GEN
522 885259 : RgX_inflate(GEN x0, long d)
523 : {
524 885259 : long i, id, dy, dx = degpol(x0);
525 885259 : GEN x = x0 + 2, z, y;
526 885259 : if (dx <= 0) return leafcopy(x0);
527 807969 : dy = dx*d;
528 807969 : y = cgetg(dy+3, t_POL); y[1] = x0[1];
529 807969 : z = y + 2;
530 28115967 : for (i=0; i<=dy; i++) gel(z,i) = gen_0;
531 11504581 : for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
532 807969 : return y;
533 : }
534 :
535 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
536 : static GEN
537 6084975 : RgX_Rg_translate_basecase(GEN P, GEN c)
538 : {
539 6084975 : pari_sp av = avma;
540 : GEN Q;
541 : long i, k, n;
542 :
543 6084975 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
544 6083046 : Q = leafcopy(P); n = degpol(P);
545 6083047 : if (isint1(c))
546 : {
547 6867 : for (i=1; i<=n; i++)
548 : {
549 20265 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gel(Q,2+k+1));
550 5390 : if (gc_needed(av,2))
551 : {
552 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(1), i = %ld/%ld", i,n);
553 0 : Q = gc_GEN(av, Q);
554 : }
555 : }
556 : }
557 6081569 : else if (isintm1(c))
558 : {
559 15246 : for (i=1; i<=n; i++)
560 : {
561 49630 : for (k=n-i; k<n; k++) gel(Q,2+k) = gsub(gel(Q,2+k), gel(Q,2+k+1));
562 12299 : if (gc_needed(av,2))
563 : {
564 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate(-1), i = %ld/%ld", i,n);
565 0 : Q = gc_GEN(av, Q);
566 : }
567 : }
568 : }
569 : else
570 : {
571 20578640 : for (i=1; i<=n; i++)
572 : {
573 47487591 : for (k=n-i; k<n; k++) gel(Q,2+k) = gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1)));
574 14500004 : if (gc_needed(av,2))
575 : {
576 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_Rg_translate, i = %ld/%ld", i,n);
577 0 : Q = gc_GEN(av, Q);
578 : }
579 : }
580 : }
581 6082811 : return gc_GEN(av, Q);
582 : }
583 :
584 : static GEN
585 83272 : zero_FpX_mod(GEN p, long v)
586 : {
587 83272 : GEN r = cgetg(3,t_POL);
588 83272 : r[1] = evalvarn(v);
589 83272 : gel(r,2) = mkintmod(gen_0, icopy(p));
590 83272 : return r;
591 : }
592 :
593 : static GEN
594 91 : zero_FpXQX_mod(GEN pol, GEN p, long v)
595 : {
596 91 : GEN r = cgetg(3,t_POL);
597 91 : r[1] = evalvarn(v);
598 91 : gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
599 91 : return r;
600 : }
601 :
602 : static GEN
603 0 : RgX_Rg_translate_FpX(GEN P, GEN c, GEN p)
604 : {
605 0 : pari_sp av = avma;
606 : GEN r;
607 : #if 0
608 : /* 'divide by 0' error if p is not prime and c not invertible */
609 : if (lgefint(p) == 3)
610 : {
611 : ulong pp = uel(p, 2);
612 : r = Flx_to_ZX_inplace(Flx_Fl_translate(RgX_to_Flx(x, pp), Rg_to_Fl(c, pp), pp));
613 : }
614 : else
615 : #endif
616 0 : r = FpX_Fp_translate(RgX_to_FpX(P, p), Rg_to_Fp(c, p), p);
617 0 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(P)); }
618 0 : return gc_upto(av, FpX_to_mod(r, p));
619 : }
620 :
621 : static GEN
622 7 : RgX_Rg_translate_FpXQX(GEN x, GEN c, GEN pol, GEN p)
623 : {
624 7 : pari_sp av = avma;
625 7 : GEN r, T = RgX_to_FpX(pol, p);
626 7 : if (signe(T)==0) pari_err_OP("subst", x, c);
627 7 : r = FpXQX_FpXQ_translate(RgX_to_FpXQX(x, T, p), Rg_to_FpXQ(c, T, p), T, p);
628 7 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
629 7 : return gc_upto(av, FpXQX_to_mod(r, T, p));
630 : }
631 :
632 : static GEN
633 6445694 : RgX_Rg_translate_fast(GEN P, GEN c)
634 : {
635 : GEN p, pol;
636 : long pa;
637 6445694 : long t = RgX_Rg_type(P, c, &p,&pol,&pa);
638 6445685 : switch(t)
639 : {
640 360892 : case t_INT: return ZX_Z_translate(P, c);
641 0 : case t_INTMOD: return RgX_Rg_translate_FpX(P, c, p);
642 6 : case RgX_type_code(t_POLMOD, t_INTMOD):
643 6 : return RgX_Rg_translate_FpXQX(P, c, pol, p);
644 6084788 : default: return NULL;
645 : }
646 : }
647 :
648 : static GEN
649 6085164 : RgX_Rg_translate_i(GEN P, GEN c)
650 : {
651 6085164 : pari_sp av = avma;
652 : long n;
653 6085164 : n = degpol(P);
654 6085163 : if (n < 40)
655 6084974 : return RgX_Rg_translate_basecase(P, c);
656 : else
657 : {
658 189 : long d = n >> 1;
659 189 : GEN Q = RgX_Rg_translate_i(RgX_shift_shallow(P, -d), c);
660 189 : GEN R = RgX_Rg_translate_i(RgXn_red_shallow(P, d), c);
661 189 : GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
662 189 : return gc_upto(av, RgX_add(RgX_mul(Q, S), R));
663 : }
664 : }
665 :
666 : GEN
667 6445696 : RgX_Rg_translate(GEN P, GEN c)
668 : {
669 6445696 : GEN R = RgX_Rg_translate_fast(P, c);
670 6445684 : return R ? R: RgX_Rg_translate_i(P,c);
671 : }
672 : /* P(ax + b) */
673 : GEN
674 30226 : RgX_affine(GEN P, GEN a, GEN b)
675 : {
676 30226 : if (!gequal0(b)) P = RgX_Rg_translate(P, b);
677 30226 : return RgX_unscale(P, a);
678 : }
679 :
680 : /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
681 : GEN
682 33584 : RgXQX_RgXQ_translate(GEN P, GEN c, GEN T)
683 : {
684 33584 : pari_sp av = avma;
685 : GEN Q;
686 : long i, k, n;
687 :
688 33584 : if (!signe(P) || gequal0(c)) return RgX_copy(P);
689 33241 : Q = leafcopy(P); n = degpol(P);
690 105778 : for (i=1; i<=n; i++)
691 : {
692 303686 : for (k=n-i; k<n; k++)
693 : {
694 231149 : pari_sp av2 = avma;
695 231149 : gel(Q,2+k) = gc_upto(av2,
696 231149 : RgX_rem(gadd(gel(Q,2+k), gmul(c, gel(Q,2+k+1))), T));
697 : }
698 72537 : if (gc_needed(av,2))
699 : {
700 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_RgXQ_translate, i = %ld/%ld", i,n);
701 0 : Q = gc_GEN(av, Q);
702 : }
703 : }
704 33241 : return gc_GEN(av, Q);
705 : }
706 :
707 : /********************************************************************/
708 : /** **/
709 : /** CONVERSIONS **/
710 : /** (not memory clean) **/
711 : /** **/
712 : /********************************************************************/
713 : /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
714 : * but everything else is */
715 : static GEN
716 168329 : QXQ_to_mod(GEN x, GEN T)
717 : {
718 : long d;
719 168329 : switch(typ(x))
720 : {
721 68391 : case t_INT: return icopy(x);
722 2079 : case t_FRAC: return gcopy(x);
723 97859 : case t_POL:
724 97859 : d = degpol(x);
725 97859 : if (d < 0) return gen_0;
726 97831 : if (d == 0) return gcopy(gel(x,2));
727 97761 : return mkpolmod(RgX_copy(x), T);
728 0 : default: pari_err_TYPE("QXQ_to_mod",x);
729 : return NULL;/* LCOV_EXCL_LINE */
730 : }
731 : }
732 : /* pure shallow version */
733 : GEN
734 840695 : QXQ_to_mod_shallow(GEN x, GEN T)
735 : {
736 : long d;
737 840695 : switch(typ(x))
738 : {
739 546003 : case t_INT:
740 546003 : case t_FRAC: return x;
741 294692 : case t_POL:
742 294692 : d = degpol(x);
743 294692 : if (d < 0) return gen_0;
744 249702 : if (d == 0) return gel(x,2);
745 233361 : return mkpolmod(x, T);
746 0 : default: pari_err_TYPE("QXQ_to_mod",x);
747 : return NULL;/* LCOV_EXCL_LINE */
748 : }
749 : }
750 : /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
751 : * Not memory clean because T not copied, but everything else is */
752 : static GEN
753 37212 : QXQX_to_mod(GEN z, GEN T)
754 : {
755 37212 : long i,l = lg(z);
756 37212 : GEN x = cgetg(l,t_POL);
757 186046 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
758 37212 : x[1] = z[1]; return normalizepol_lg(x,l);
759 : }
760 : /* pure shallow version */
761 : GEN
762 204096 : QXQX_to_mod_shallow(GEN z, GEN T)
763 : {
764 204096 : long i,l = lg(z);
765 204096 : GEN x = cgetg(l,t_POL);
766 961337 : for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
767 204096 : x[1] = z[1]; return normalizepol_lg(x,l);
768 : }
769 : /* Apply QXQX_to_mod to all entries. Memory-clean ! */
770 : GEN
771 12572 : QXQXV_to_mod(GEN V, GEN T)
772 : {
773 12572 : long i, l = lg(V);
774 12572 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
775 49784 : for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
776 12572 : return z;
777 : }
778 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
779 : GEN
780 22343 : QXQV_to_mod(GEN V, GEN T)
781 : {
782 22343 : long i, l = lg(V);
783 22343 : GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
784 41838 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
785 22343 : return z;
786 : }
787 :
788 : /* Apply QXQ_to_mod to all entries. Memory-clean ! */
789 : GEN
790 14854 : QXQC_to_mod_shallow(GEN V, GEN T)
791 : {
792 14854 : long i, l = lg(V);
793 14854 : GEN z = cgetg(l, t_COL);
794 98308 : for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
795 14854 : return z;
796 : }
797 :
798 : GEN
799 6720 : QXQM_to_mod_shallow(GEN V, GEN T)
800 : {
801 6720 : long i, l = lg(V);
802 6720 : GEN z = cgetg(l, t_MAT);
803 21574 : for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
804 6720 : return z;
805 : }
806 :
807 : GEN
808 8301660 : RgX_renormalize_lg(GEN x, long lx)
809 : {
810 : long i;
811 11350965 : for (i = lx-1; i>1; i--)
812 10946769 : if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
813 8301661 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
814 8301661 : setlg(x, i+1); setsigne(x, i != 1); return x;
815 : }
816 :
817 : GEN
818 2283429 : RgV_to_RgX(GEN x, long v)
819 : {
820 2283429 : long i, k = lg(x);
821 : GEN p;
822 :
823 6357476 : while (--k && gequal0(gel(x,k)));
824 2283429 : if (!k) return pol_0(v);
825 2250420 : i = k+2; p = cgetg(i,t_POL);
826 2250412 : p[1] = evalsigne(1) | evalvarn(v);
827 24386140 : x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
828 2250412 : return p;
829 : }
830 : GEN
831 188283 : RgV_to_RgX_reverse(GEN x, long v)
832 : {
833 188283 : long j, k, l = lg(x);
834 : GEN p;
835 :
836 189851 : for (k = 1; k < l; k++)
837 189851 : if (!gequal0(gel(x,k))) break;
838 188283 : if (k == l) return pol_0(v);
839 188283 : k -= 1;
840 188283 : l -= k;
841 188283 : x += k;
842 188283 : p = cgetg(l+1,t_POL);
843 188283 : p[1] = evalsigne(1) | evalvarn(v);
844 985038 : for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
845 188283 : return p;
846 : }
847 :
848 : /* return the (N-dimensional) vector of coeffs of p */
849 : GEN
850 15338902 : RgX_to_RgC(GEN x, long N)
851 : {
852 : long i, l;
853 : GEN z;
854 15338902 : l = lg(x)-1; x++;
855 15338902 : if (l > N+1) l = N+1; /* truncate higher degree terms */
856 15338902 : z = cgetg(N+1,t_COL);
857 112252165 : for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
858 31410792 : for ( ; i<=N; i++) gel(z,i) = gen_0;
859 15339031 : return z;
860 : }
861 : GEN
862 1360199 : Rg_to_RgC(GEN x, long N)
863 : {
864 1360199 : return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
865 : }
866 :
867 : /* vector of polynomials (in v) whose coefs are given by the columns of x */
868 : GEN
869 302956 : RgM_to_RgXV(GEN x, long v)
870 1289731 : { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
871 : GEN
872 7202 : RgM_to_RgXV_reverse(GEN x, long v)
873 28808 : { pari_APPLY_type(t_VEC, RgV_to_RgX_reverse(gel(x,i), v)) }
874 :
875 : /* matrix whose entries are given by the coeffs of the polynomials in
876 : * vector v (considered as degree n-1 polynomials) */
877 : GEN
878 336530 : RgV_to_RgM(GEN x, long n)
879 1691536 : { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
880 :
881 : GEN
882 79011 : RgXV_to_RgM(GEN x, long n)
883 401767 : { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
884 :
885 : /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
886 : GEN
887 23714 : RgM_to_RgXX(GEN x, long v,long w)
888 : {
889 23714 : long j, lx = lg(x);
890 23714 : GEN y = cgetg(lx+1, t_POL);
891 23714 : y[1] = evalsigne(1) | evalvarn(v);
892 23714 : y++;
893 132085 : for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
894 23714 : return normalizepol_lg(--y, lx+1);
895 : }
896 :
897 : /* matrix whose entries are given by the coeffs of the polynomial v in
898 : * two variables (considered as degree n-1 polynomials) */
899 : GEN
900 322 : RgXX_to_RgM(GEN v, long n)
901 : {
902 322 : long j, N = lg(v)-1;
903 322 : GEN y = cgetg(N, t_MAT);
904 1043 : for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
905 322 : return y;
906 : }
907 :
908 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
909 : GEN
910 33141 : RgXY_swapspec(GEN x, long n, long w, long nx)
911 : {
912 33141 : long j, ly = n+3;
913 33141 : GEN y = cgetg(ly, t_POL);
914 33141 : y[1] = evalsigne(1);
915 401591 : for (j=2; j<ly; j++)
916 : {
917 : long k;
918 368450 : GEN a = cgetg(nx+2,t_POL);
919 368450 : a[1] = evalsigne(1) | evalvarn(w);
920 2026528 : for (k=0; k<nx; k++)
921 : {
922 1658078 : GEN xk = gel(x,k);
923 1658078 : if (typ(xk)==t_POL && varn(xk)==w)
924 1558059 : gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
925 : else
926 100019 : gel(a,k+2) = j==2 ? xk: gen_0;
927 : }
928 368450 : gel(y,j) = normalizepol_lg(a, nx+2);
929 : }
930 33141 : return normalizepol_lg(y,ly);
931 : }
932 :
933 : /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
934 : GEN
935 2051 : RgXY_swap(GEN x, long n, long w)
936 : {
937 2051 : GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
938 2051 : setvarn(z, varn(x)); return z;
939 : }
940 :
941 : long
942 2022605 : RgXY_degreex(GEN b)
943 : {
944 2022605 : long deg = 0, i;
945 2022605 : if (!signe(b)) return -1;
946 8497731 : for (i = 2; i < lg(b); ++i)
947 : {
948 6475127 : GEN bi = gel(b, i);
949 6475127 : if (typ(bi) == t_POL)
950 1023536 : deg = maxss(deg, degpol(bi));
951 : }
952 2022604 : return deg;
953 : }
954 :
955 : GEN
956 38738 : RgXY_derivx(GEN x) { pari_APPLY_pol(RgX_deriv(gel(x,i))); }
957 :
958 : /* return (x % X^n). Shallow */
959 : GEN
960 8301966 : RgXn_red_shallow(GEN a, long n)
961 : {
962 8301966 : long i, L = n+2, l = lg(a);
963 : GEN b;
964 8301966 : if (L >= l) return a; /* deg(x) < n */
965 5999117 : b = cgetg(L, t_POL); b[1] = a[1];
966 38226093 : for (i=2; i<L; i++) gel(b,i) = gel(a,i);
967 5999117 : return normalizepol_lg(b,L);
968 : }
969 :
970 : GEN
971 483 : RgXnV_red_shallow(GEN x, long n)
972 2268 : { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
973 :
974 : /* return (x * X^n). Shallow */
975 : GEN
976 175797338 : RgX_shift_shallow(GEN a, long n)
977 : {
978 175797338 : long i, l = lg(a);
979 : GEN b;
980 175797338 : if (l == 2 || !n) return a;
981 109489658 : l += n;
982 109489658 : if (n < 0)
983 : {
984 54102544 : if (l <= 2) return pol_0(varn(a));
985 52538558 : b = cgetg(l, t_POL); b[1] = a[1];
986 52539350 : a -= n;
987 163535186 : for (i=2; i<l; i++) gel(b,i) = gel(a,i);
988 : } else {
989 55387114 : b = cgetg(l, t_POL); b[1] = a[1];
990 55390942 : a -= n; n += 2;
991 120452302 : for (i=2; i<n; i++) gel(b,i) = gen_0;
992 213289773 : for ( ; i<l; i++) gel(b,i) = gel(a,i);
993 : }
994 107930292 : return b;
995 : }
996 : /* return (x * X^n). */
997 : GEN
998 1577735 : RgX_shift(GEN a, long n)
999 : {
1000 1577735 : long i, l = lg(a);
1001 : GEN b;
1002 1577735 : if (l == 2 || !n) return RgX_copy(a);
1003 1576783 : l += n;
1004 1576783 : if (n < 0)
1005 : {
1006 1442 : if (l <= 2) return pol_0(varn(a));
1007 1372 : b = cgetg(l, t_POL); b[1] = a[1];
1008 1372 : a -= n;
1009 9534 : for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1010 : } else {
1011 1575341 : b = cgetg(l, t_POL); b[1] = a[1];
1012 1575341 : a -= n; n += 2;
1013 3976367 : for (i=2; i<n; i++) gel(b,i) = gen_0;
1014 4344046 : for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
1015 : }
1016 1576713 : return b;
1017 : }
1018 :
1019 : GEN
1020 317037 : RgX_rotate_shallow(GEN P, long k, long p)
1021 : {
1022 317037 : long i, l = lgpol(P);
1023 : GEN r;
1024 317037 : if (signe(P)==0)
1025 1365 : return pol_0(varn(P));
1026 315672 : r = cgetg(p+2,t_POL); r[1] = P[1];
1027 2100644 : for(i=0; i<p; i++)
1028 : {
1029 1784972 : long s = 2+(i+k)%p;
1030 1784972 : gel(r,s) = i<l? gel(P,2+i): gen_0;
1031 : }
1032 315672 : return RgX_renormalize(r);
1033 : }
1034 :
1035 : GEN
1036 2997517 : RgX_mulXn(GEN x, long d)
1037 : {
1038 : pari_sp av;
1039 : GEN z;
1040 : long v;
1041 2997517 : if (d >= 0) return RgX_shift(x, d);
1042 1472762 : d = -d;
1043 1472762 : v = RgX_val(x);
1044 1472762 : if (v >= d) return RgX_shift(x, -d);
1045 1472748 : av = avma;
1046 1472748 : z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
1047 1472748 : return gc_upto(av, z);
1048 : }
1049 :
1050 : long
1051 588 : RgXV_maxdegree(GEN x)
1052 : {
1053 588 : long d = -1, i, l = lg(x);
1054 4494 : for (i = 1; i < l; i++)
1055 3906 : d = maxss(d, degpol(gel(x,i)));
1056 588 : return d;
1057 : }
1058 :
1059 : long
1060 3984033 : RgX_val(GEN x)
1061 : {
1062 3984033 : long i, lx = lg(x);
1063 3984033 : if (lx == 2) return LONG_MAX;
1064 4948367 : for (i = 2; i < lx; i++)
1065 4948325 : if (!isexactzero(gel(x,i))) break;
1066 3983879 : if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
1067 3983837 : return i - 2;
1068 : }
1069 : long
1070 80780249 : RgX_valrem(GEN x, GEN *Z)
1071 : {
1072 80780249 : long v, i, lx = lg(x);
1073 80780249 : if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
1074 125154277 : for (i = 2; i < lx; i++)
1075 125155476 : if (!isexactzero(gel(x,i))) break;
1076 : /* possible with nonrational zeros */
1077 80780132 : if (i == lx)
1078 : {
1079 21 : *Z = scalarpol_shallow(Rg_get_0(x), varn(x));
1080 21 : return LONG_MAX;
1081 : }
1082 80780111 : v = i - 2;
1083 80780111 : *Z = RgX_shift_shallow(x, -v);
1084 80791940 : return v;
1085 : }
1086 : long
1087 882625 : RgX_valrem_inexact(GEN x, GEN *Z)
1088 : {
1089 : long v;
1090 882625 : if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
1091 901089 : for (v = 0;; v++)
1092 901089 : if (!gequal0(gel(x,2+v))) break;
1093 882618 : if (Z) *Z = RgX_shift_shallow(x, -v);
1094 882618 : return v;
1095 : }
1096 :
1097 : GEN
1098 68187 : RgXQC_red(GEN x, GEN T)
1099 414505 : { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
1100 :
1101 : GEN
1102 1449 : RgXQV_red(GEN x, GEN T)
1103 35007 : { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
1104 :
1105 : GEN
1106 13195 : RgXQM_red(GEN x, GEN T)
1107 81382 : { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
1108 :
1109 : GEN
1110 322 : RgXQM_mul(GEN P, GEN Q, GEN T)
1111 : {
1112 322 : return RgXQM_red(RgM_mul(P, Q), T);
1113 : }
1114 :
1115 : GEN
1116 497492 : RgXQX_red(GEN P, GEN T)
1117 : {
1118 497492 : long i, l = lg(P);
1119 497492 : GEN Q = cgetg(l, t_POL);
1120 497492 : Q[1] = P[1];
1121 2590665 : for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
1122 497492 : return normalizepol_lg(Q, l);
1123 : }
1124 :
1125 : GEN
1126 831446 : RgX_deriv(GEN x)
1127 : {
1128 831446 : long i,lx = lg(x)-1;
1129 : GEN y;
1130 :
1131 831446 : if (lx<3) return pol_0(varn(x));
1132 828464 : y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
1133 3656483 : for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
1134 828456 : y[1] = x[1]; return normalizepol_lg(y,i);
1135 : }
1136 :
1137 : GEN
1138 798 : RgX_derivn(GEN x, long n)
1139 : {
1140 798 : long i, vx = varn(x), lx = lg(x);
1141 : GEN y;
1142 798 : if (lx <= n+2) return pol_0(vx);
1143 749 : lx -= n;
1144 749 : y = cgetg(lx,t_POL);
1145 749 : y[1] = evalsigne(1)| evalvarn(vx);
1146 50904 : for (i=0; i<lx-2; i++)
1147 50155 : gel(y,i+2) = gmul(mulu_interval(i+1,i+n),gel(x,i+2+n));
1148 749 : return normalizepol_lg(y, lx);
1149 : }
1150 :
1151 : GEN
1152 0 : RgX_homogenous_derivn(GEN x, long n, long d)
1153 : {
1154 0 : long i, vx = varn(x), lx = minss(lg(x), d+3-n);
1155 : GEN y;
1156 0 : if (d < n) return pol_0(vx);
1157 0 : y = cgetg(lx, t_POL);
1158 0 : y[1] = evalsigne(1)| evalvarn(vx);
1159 0 : for (i = 0; i < lx-2; i++)
1160 0 : gel(y,2+i) = gmul(mulu_interval(d-i+1-n,d-i), gel(x,2+i));
1161 0 : return normalizepol_lg(y, lx);
1162 : }
1163 :
1164 : GEN
1165 2586560 : RgX_recipspec_shallow(GEN x, long l, long n)
1166 : {
1167 : long i;
1168 2586560 : GEN z = cgetg(n+2,t_POL);
1169 2586567 : z[1] = 0; z += 2;
1170 144317999 : for(i=0; i<l; i++) gel(z,n-i-1) = gel(x,i);
1171 2835219 : for( ; i<n; i++) gel(z, n-i-1) = gen_0;
1172 2586567 : return normalizepol_lg(z-2,n+2);
1173 : }
1174 :
1175 : GEN
1176 645781 : RgXn_recip_shallow(GEN P, long n)
1177 : {
1178 645781 : GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1179 645793 : setvarn(Q, varn(P));
1180 645793 : return Q;
1181 : }
1182 :
1183 : /* return coefficients s.t x = x_0 X^n + ... + x_n */
1184 : GEN
1185 35301 : RgX_recip(GEN x)
1186 : {
1187 : long lx, i, j;
1188 35301 : GEN y = cgetg_copy(x, &lx);
1189 295071 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1190 35301 : return normalizepol_lg(y,lx);
1191 : }
1192 : /* shallow version */
1193 : GEN
1194 59388 : RgX_recip_shallow(GEN x)
1195 : {
1196 : long lx, i, j;
1197 59388 : GEN y = cgetg_copy(x, &lx);
1198 356118 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1199 59388 : return normalizepol_lg(y,lx);
1200 : }
1201 :
1202 : GEN
1203 3756394 : RgX_recip_i(GEN x)
1204 : {
1205 : long lx, i, j;
1206 3756394 : GEN y = cgetg_copy(x, &lx);
1207 20259541 : y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1208 3756386 : return y;
1209 : }
1210 : /*******************************************************************/
1211 : /* */
1212 : /* ADDITION / SUBTRACTION */
1213 : /* */
1214 : /*******************************************************************/
1215 : /* cf RgX_coeff_simplify */
1216 : INLINE GEN
1217 441813691 : RgX_coeff_add(GEN a, GEN b)
1218 : {
1219 441813691 : pari_sp av = avma;
1220 441813691 : GEN s = gadd(a, b);
1221 441774263 : if (typ(s) != t_POL) return s;
1222 50947327 : switch(lg(s))
1223 : {
1224 1505140 : case 2: return gc_const(av, gen_0);
1225 131919 : case 3: return gc_GEN(av, gel(s,2));
1226 49310268 : default: return s;
1227 : }
1228 : }
1229 : INLINE GEN
1230 247257640 : RgX_coeff_sub(GEN a, GEN b)
1231 : {
1232 247257640 : pari_sp av = avma;
1233 247257640 : GEN s = gsub(a, b);
1234 247131015 : if (typ(s) != t_POL) return s;
1235 9753414 : switch(lg(s))
1236 : {
1237 2275051 : case 2: return gc_const(av, gen_0);
1238 250071 : case 3: return gc_GEN(av, gel(s,2));
1239 7228292 : default: return s;
1240 : }
1241 : }
1242 : /* same variable */
1243 : GEN
1244 118065957 : RgX_add(GEN x, GEN y)
1245 : {
1246 118065957 : long i, lx = lg(x), ly = lg(y);
1247 : GEN z;
1248 118065957 : if (ly <= lx) {
1249 104828385 : z = cgetg(lx,t_POL); z[1] = x[1];
1250 481623352 : for (i=2; i < ly; i++) gel(z,i) = RgX_coeff_add(gel(x,i),gel(y,i));
1251 162112784 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1252 104803950 : z = normalizepol_lg(z, lx);
1253 : } else {
1254 13237572 : z = cgetg(ly,t_POL); z[1] = y[1];
1255 52608135 : for (i=2; i < lx; i++) gel(z,i) = RgX_coeff_add(gel(x,i),gel(y,i));
1256 38506155 : for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1257 13238825 : z = normalizepol_lg(z, ly);
1258 : }
1259 118057521 : return z;
1260 : }
1261 : GEN
1262 69191393 : RgX_sub(GEN x, GEN y)
1263 : {
1264 69191393 : long i, lx = lg(x), ly = lg(y);
1265 : GEN z;
1266 69191393 : if (ly <= lx) {
1267 35280878 : z = cgetg(lx,t_POL); z[1] = x[1];
1268 185323267 : for (i=2; i < ly; i++) gel(z,i) = RgX_coeff_sub(gel(x,i),gel(y,i));
1269 61465307 : for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1270 35259079 : z = normalizepol_lg(z, lx);
1271 : } else {
1272 33910515 : z = cgetg(ly,t_POL); z[1] = y[1];
1273 126093800 : for (i=2; i < lx; i++) gel(z,i) = RgX_coeff_sub(gel(x,i),gel(y,i));
1274 74706421 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1275 33878680 : z = normalizepol_lg(z, ly);
1276 : }
1277 69189525 : return z;
1278 : }
1279 : GEN
1280 6522081 : RgX_neg(GEN x)
1281 46307327 : { pari_APPLY_pol_normalized(gneg(gel(x,i))); }
1282 :
1283 : GEN
1284 26293531 : RgX_Rg_add(GEN y, GEN x)
1285 : {
1286 : GEN z;
1287 26293531 : long lz = lg(y), i;
1288 26293531 : if (lz == 2) return scalarpol(x,varn(y));
1289 25578197 : z = cgetg(lz,t_POL); z[1] = y[1];
1290 25578143 : gel(z,2) = RgX_coeff_add(gel(y,2),x);
1291 86286702 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1292 : /* probably useless unless lz = 3, but cannot be skipped if y is
1293 : * an inexact 0 */
1294 25578125 : return normalizepol_lg(z,lz);
1295 : }
1296 : GEN
1297 65086 : RgX_Rg_add_shallow(GEN y, GEN x)
1298 : {
1299 : GEN z;
1300 65086 : long lz = lg(y), i;
1301 65086 : if (lz == 2) return scalarpol(x,varn(y));
1302 65086 : z = cgetg(lz,t_POL); z[1] = y[1];
1303 65086 : gel(z,2) = RgX_coeff_add(gel(y,2),x);
1304 130316 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1305 65086 : return normalizepol_lg(z,lz);
1306 : }
1307 : GEN
1308 223144 : RgX_Rg_sub(GEN y, GEN x)
1309 : {
1310 : GEN z;
1311 223144 : long lz = lg(y), i;
1312 223144 : if (lz == 2)
1313 : { /* scalarpol(gneg(x),varn(y)) optimized */
1314 133 : long v = varn(y);
1315 133 : if (isrationalzero(x)) return pol_0(v);
1316 7 : z = cgetg(3,t_POL);
1317 7 : z[1] = gequal0(x)? evalvarn(v)
1318 7 : : evalvarn(v) | evalsigne(1);
1319 7 : gel(z,2) = gneg(x); return z;
1320 : }
1321 223011 : z = cgetg(lz,t_POL); z[1] = y[1];
1322 223011 : gel(z,2) = RgX_coeff_sub(gel(y,2),x);
1323 589558 : for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1324 223010 : return normalizepol_lg(z,lz);
1325 : }
1326 : GEN
1327 4886741 : Rg_RgX_sub(GEN x, GEN y)
1328 : {
1329 : GEN z;
1330 4886741 : long lz = lg(y), i;
1331 4886741 : if (lz == 2) return scalarpol(x,varn(y));
1332 4886601 : z = cgetg(lz,t_POL); z[1] = y[1];
1333 4886587 : gel(z,2) = RgX_coeff_sub(x, gel(y,2));
1334 7341645 : for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1335 4886563 : return normalizepol_lg(z,lz);
1336 : }
1337 : /*******************************************************************/
1338 : /* */
1339 : /* KARATSUBA MULTIPLICATION */
1340 : /* */
1341 : /*******************************************************************/
1342 : #if 0
1343 : /* to debug Karatsuba-like routines */
1344 : GEN
1345 : zx_debug_spec(GEN x, long nx)
1346 : {
1347 : GEN z = cgetg(nx+2,t_POL);
1348 : long i;
1349 : for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1350 : z[1] = evalsigne(1); return z;
1351 : }
1352 :
1353 : GEN
1354 : RgX_debug_spec(GEN x, long nx)
1355 : {
1356 : GEN z = cgetg(nx+2,t_POL);
1357 : long i;
1358 : for (i=0; i<nx; i++) z[i+2] = x[i];
1359 : z[1] = evalsigne(1); return z;
1360 : }
1361 : #endif
1362 :
1363 : /* generic multiplication */
1364 : GEN
1365 9061967 : RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1366 : {
1367 : GEN z, t;
1368 : long i;
1369 9061967 : if (nx == ny) {
1370 1579319 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1371 4997398 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1372 1579311 : return normalizepol_lg(z, nx+2);
1373 : }
1374 7482648 : if (ny < nx) {
1375 7297293 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1376 26682533 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1377 17394368 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1378 7296762 : return normalizepol_lg(z, nx+2);
1379 : } else {
1380 185355 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1381 3672730 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1382 447632 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1383 185366 : return normalizepol_lg(z, ny+2);
1384 : }
1385 : }
1386 : GEN
1387 222463 : RgX_addspec(GEN x, GEN y, long nx, long ny)
1388 : {
1389 : GEN z, t;
1390 : long i;
1391 222463 : if (nx == ny) {
1392 12824 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1393 2185778 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1394 12824 : return normalizepol_lg(z, nx+2);
1395 : }
1396 209639 : if (ny < nx) {
1397 207854 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1398 3725567 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1399 2373010 : for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1400 207853 : return normalizepol_lg(z, nx+2);
1401 : } else {
1402 1785 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1403 331478 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1404 12236 : for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1405 1785 : return normalizepol_lg(z, ny+2);
1406 : }
1407 : }
1408 :
1409 : /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1410 : * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1411 : * t_VECSMALL, to hold this, which can be left on stack: GC functions
1412 : * will not crash on it. The returned vector itself is not a proper GEN,
1413 : * we access the coefficients as x[i], i = 0..deg(x) */
1414 : static GEN
1415 82439150 : RgXspec_kill0(GEN x, long lx)
1416 : {
1417 82439150 : GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit GC-wise */
1418 : long i;
1419 251921309 : for (i=0; i <lx; i++)
1420 : {
1421 169482108 : GEN c = gel(x,i);
1422 169482108 : z[i] = (long)(isrationalzero(c)? NULL: c);
1423 : }
1424 82439201 : return z;
1425 : }
1426 :
1427 : /* For consistency we want operations on t_POL of the same variable v
1428 : * to return a t_POL in v; but this becomes a nuisance for coefficients,
1429 : * e.g., multivariate t_POL, which should have simplest possible types */
1430 : INLINE GEN
1431 113647662 : RgX_coeff_simplify(pari_sp av, GEN s)
1432 : {
1433 113647662 : if (!s) return gen_0;
1434 106597599 : if (typ(s) != t_POL) return gc_upto(av, s);
1435 24893087 : switch(lg(s))
1436 : {
1437 12222 : case 2: return gen_0;
1438 2152039 : case 3: return gc_GEN(av, gel(s,2));
1439 22728826 : default: return gc_upto(av, s);
1440 : }
1441 : }
1442 : INLINE GEN
1443 113557764 : RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1444 : {
1445 113557764 : pari_sp av = avma;
1446 113557764 : GEN s = NULL;
1447 : long i;
1448 :
1449 403067515 : for (i=a; i<b; i++)
1450 289519215 : if (gel(y,i) && gel(x,-i))
1451 : {
1452 211023808 : GEN t = gmul(gel(y,i), gel(x,-i));
1453 211019319 : s = s? gadd(s, t): t;
1454 : }
1455 113548300 : return RgX_coeff_simplify(av, s);
1456 : }
1457 :
1458 : /* assume nx >= ny > 0, return x * y * t^v */
1459 : static GEN
1460 33611288 : RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1461 : {
1462 : long i, lz, nz;
1463 : GEN z;
1464 :
1465 33611288 : x = RgXspec_kill0(x,nx);
1466 33611277 : y = RgXspec_kill0(y,ny);
1467 33611257 : lz = nx + ny + 1; nz = lz-2;
1468 33611257 : lz += v;
1469 33611257 : z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1470 71505326 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1471 82098969 : for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1472 55411136 : for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1473 48487796 : for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1474 33610778 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1475 : }
1476 :
1477 : /* return (x * X^d) + y. Assume d > 0 */
1478 : GEN
1479 9743626 : RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1480 : {
1481 : GEN x, y, xd, yd, zd;
1482 : long a, lz, nx, ny;
1483 :
1484 9743626 : if (!signe(x0)) return y0;
1485 9588217 : ny = lgpol(y0);
1486 9588223 : nx = lgpol(x0);
1487 9588306 : zd = (GEN)avma;
1488 9588306 : x = x0 + 2; y = y0 + 2; a = ny-d;
1489 9588306 : if (a <= 0)
1490 : {
1491 1482369 : lz = nx+d+2;
1492 1482369 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1493 3406015 : while (xd > x) gel(--zd,0) = gel(--xd,0);
1494 1482368 : x = zd + a;
1495 1502907 : while (zd > x) gel(--zd,0) = gen_0;
1496 : }
1497 : else
1498 : {
1499 8105937 : xd = new_chunk(d); yd = y+d;
1500 8105938 : x = RgX_addspec_shallow(x,yd, nx,a);
1501 8105848 : lz = (a>nx)? ny+2: lg(x)+d;
1502 40009826 : x += 2; while (xd > x) *--zd = *--xd;
1503 : }
1504 23285805 : while (yd > y) *--zd = *--yd;
1505 9588216 : *--zd = x0[1];
1506 9588216 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1507 : }
1508 : GEN
1509 514933 : RgX_addmulXn(GEN x0, GEN y0, long d)
1510 : {
1511 : GEN x, y, xd, yd, zd;
1512 : long a, lz, nx, ny;
1513 :
1514 514933 : if (!signe(x0)) return RgX_copy(y0);
1515 514149 : nx = lgpol(x0);
1516 514149 : ny = lgpol(y0);
1517 514149 : zd = (GEN)avma;
1518 514149 : x = x0 + 2; y = y0 + 2; a = ny-d;
1519 514149 : if (a <= 0)
1520 : {
1521 291686 : lz = nx+d+2;
1522 291686 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1523 4293251 : while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1524 291686 : x = zd + a;
1525 757796 : while (zd > x) gel(--zd,0) = gen_0;
1526 : }
1527 : else
1528 : {
1529 222463 : xd = new_chunk(d); yd = y+d;
1530 222463 : x = RgX_addspec(x,yd, nx,a);
1531 222463 : lz = (a>nx)? ny+2: lg(x)+d;
1532 8418433 : x += 2; while (xd > x) *--zd = *--xd;
1533 : }
1534 2616093 : while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1535 514149 : *--zd = x0[1];
1536 514149 : *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1537 : }
1538 :
1539 : /* return x * y mod t^n */
1540 : static GEN
1541 7591743 : RgXn_mul_basecase(GEN x, GEN y, long n)
1542 : {
1543 7591743 : long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1544 : GEN z;
1545 7591743 : if (lx < 0) return pol_0(varn(x));
1546 7591743 : if (ly < 0) return pol_0(varn(x));
1547 7591743 : z = cgetg(lz, t_POL) + 2;
1548 7591743 : x+=2; if (lx > n) lx = n;
1549 7591743 : y+=2; if (ly > n) ly = n;
1550 7591743 : z[-1] = x[-1];
1551 7591743 : if (ly > lx) { swap(x,y); lswap(lx,ly); }
1552 7591743 : x = RgXspec_kill0(x, lx);
1553 7591743 : y = RgXspec_kill0(y, ly);
1554 : /* x:y:z [i] = term of degree i */
1555 29888193 : for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1556 13639717 : for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1557 7640692 : for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1558 7591743 : return normalizepol_lg(z - 2, lz);
1559 : }
1560 : /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1561 : static GEN
1562 10559969 : RgXn_mul2(GEN f, GEN g, long n)
1563 : {
1564 10559969 : pari_sp av = avma;
1565 : GEN fe,fo, ge,go, l,h,m;
1566 : long n0, n1;
1567 10559969 : if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1568 7627240 : if (n < 80) return RgXn_mul_basecase(f,g,n);
1569 35497 : n0 = n>>1; n1 = n-n0;
1570 35497 : RgX_even_odd(f, &fe, &fo);
1571 35497 : RgX_even_odd(g, &ge, &go);
1572 35497 : l = RgXn_mul2(fe,ge,n1);
1573 35497 : h = RgXn_mul2(fo,go,n0);
1574 35497 : m = RgX_sub(RgXn_mul2(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1575 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1576 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1577 35497 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1578 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1579 35497 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1580 35497 : m = RgX_inflate(m,2);
1581 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1582 35497 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1583 35497 : h = RgX_inflate(h,2);
1584 35497 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1585 35497 : return gc_upto(av, h);
1586 : }
1587 : /* (f*g) \/ x^n */
1588 : static GEN
1589 1801193 : RgX_mulhigh_i2(GEN f, GEN g, long n)
1590 : {
1591 1801193 : long d = degpol(f)+degpol(g) + 1 - n;
1592 : GEN h;
1593 1801193 : if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1594 29654 : h = RgX_recip_i(RgXn_mul2(RgX_recip_i(f),
1595 : RgX_recip_i(g), d));
1596 29654 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1597 : }
1598 :
1599 : /* (f*g) \/ x^n */
1600 : static GEN
1601 0 : RgX_sqrhigh_i2(GEN f, long n)
1602 : {
1603 0 : long d = 2*degpol(f)+ 1 - n;
1604 : GEN h;
1605 0 : if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1606 0 : h = RgX_recip_i(RgXn_sqr(RgX_recip_i(f), d));
1607 0 : return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1608 : }
1609 :
1610 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1611 : * b+2 were sent instead. na, nb = number of terms of a, b.
1612 : * Only c, c0, c1, c2 are genuine GEN.
1613 : */
1614 : GEN
1615 34251920 : RgX_mulspec(GEN a, GEN b, long na, long nb)
1616 : {
1617 : GEN a0, c, c0;
1618 34251920 : long n0, n0a, i, v = 0;
1619 : pari_sp av;
1620 :
1621 52326638 : while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1622 54510617 : while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1623 34251895 : if (na < nb) swapspec(a,b, na,nb);
1624 34251895 : if (!nb) return pol_0(0);
1625 :
1626 34090667 : if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1627 479390 : RgX_shift_inplace_init(v);
1628 479397 : i = (na>>1); n0 = na-i; na = i;
1629 479397 : av = avma; a0 = a+n0; n0a = n0;
1630 1334256 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1631 :
1632 479397 : if (nb > n0)
1633 : {
1634 : GEN b0,c1,c2;
1635 : long n0b;
1636 :
1637 478014 : nb -= n0; b0 = b+n0; n0b = n0;
1638 1445144 : while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1639 478014 : c = RgX_mulspec(a,b,n0a,n0b);
1640 478014 : c0 = RgX_mulspec(a0,b0, na,nb);
1641 :
1642 478014 : c2 = RgX_addspec_shallow(a0,a, na,n0a);
1643 478014 : c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1644 :
1645 478014 : c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1646 478014 : c2 = RgX_sub(c1, RgX_add(c0,c));
1647 478014 : c0 = RgX_addmulXn_shallow(c0, c2, n0);
1648 : }
1649 : else
1650 : {
1651 1383 : c = RgX_mulspec(a,b,n0a,nb);
1652 1383 : c0 = RgX_mulspec(a0,b,na,nb);
1653 : }
1654 479397 : c0 = RgX_addmulXn(c0,c,n0);
1655 479397 : return RgX_shift_inplace(gc_upto(av,c0), v);
1656 : }
1657 :
1658 : INLINE GEN
1659 94825 : RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1660 : {
1661 94825 : pari_sp av = avma;
1662 94825 : GEN s = NULL;
1663 94825 : long j, l = (i+1)>>1;
1664 197435 : for (j=a; j<l; j++)
1665 : {
1666 102610 : GEN xj = gel(x,j), xx = gel(x,i-j);
1667 102610 : if (xj && xx)
1668 : {
1669 93475 : GEN t = gmul(xj, xx);
1670 93475 : s = s? gadd(s, t): t;
1671 : }
1672 : }
1673 94825 : if (s) s = gshift(s,1);
1674 94825 : if ((i&1) == 0)
1675 : {
1676 64050 : GEN t = gel(x, i>>1);
1677 64050 : if (t) {
1678 60991 : t = gsqr(t);
1679 60991 : s = s? gadd(s, t): t;
1680 : }
1681 : }
1682 94825 : return RgX_coeff_simplify(av, s);
1683 : }
1684 : static GEN
1685 33275 : RgX_sqrspec_basecase(GEN x, long nx, long v)
1686 : {
1687 : long i, lz, nz;
1688 : GEN z;
1689 :
1690 33275 : if (!nx) return pol_0(0);
1691 33275 : x = RgXspec_kill0(x,nx);
1692 33275 : lz = (nx << 1) + 1, nz = lz-2;
1693 33275 : lz += v;
1694 33275 : z = cgetg(lz,t_POL) + 2;
1695 88911 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1696 97325 : for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1697 64050 : for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1698 33275 : z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1699 : }
1700 : /* return x^2 mod t^n */
1701 : static GEN
1702 0 : RgXn_sqr_basecase(GEN x, long n)
1703 : {
1704 0 : long i, lz = n+2, lx = lgpol(x);
1705 : GEN z;
1706 0 : if (lx < 0) return pol_0(varn(x));
1707 0 : z = cgetg(lz, t_POL);
1708 0 : z[1] = x[1];
1709 0 : x+=2; if (lx > n) lx = n;
1710 0 : x = RgXspec_kill0(x,lx);
1711 0 : z+=2;/* x:z [i] = term of degree i */
1712 0 : for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1713 0 : for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1714 0 : z -= 2; return normalizepol_lg(z, lz);
1715 : }
1716 : /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1717 : static GEN
1718 315 : RgXn_sqr2(GEN f, long n)
1719 : {
1720 315 : pari_sp av = avma;
1721 : GEN fe,fo, l,h,m;
1722 : long n0, n1;
1723 315 : if (2*degpol(f) < n) return RgX_sqr_i(f);
1724 0 : if (n < 80) return RgXn_sqr_basecase(f,n);
1725 0 : n0 = n>>1; n1 = n-n0;
1726 0 : RgX_even_odd(f, &fe, &fo);
1727 0 : l = RgXn_sqr(fe,n1);
1728 0 : h = RgXn_sqr(fo,n0);
1729 0 : m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1730 : /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1731 : * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1732 0 : l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1733 : /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1734 0 : if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1735 0 : m = RgX_inflate(m,2);
1736 : /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1737 0 : if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1738 0 : h = RgX_inflate(h,2);
1739 0 : h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1740 0 : return gc_upto(av, h);
1741 : }
1742 : GEN
1743 33314 : RgX_sqrspec(GEN a, long na)
1744 : {
1745 : GEN a0, c, c0, c1;
1746 33314 : long n0, n0a, i, v = 0;
1747 : pari_sp av;
1748 :
1749 61132 : while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1750 33314 : if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1751 39 : RgX_shift_inplace_init(v);
1752 39 : i = (na>>1); n0 = na-i; na = i;
1753 39 : av = avma; a0 = a+n0; n0a = n0;
1754 39 : while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1755 :
1756 39 : c = RgX_sqrspec(a,n0a);
1757 39 : c0 = RgX_sqrspec(a0,na);
1758 39 : c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1759 39 : c0 = RgX_addmulXn_shallow(c0,c1, n0);
1760 39 : c0 = RgX_addmulXn(c0,c,n0);
1761 39 : return RgX_shift_inplace(gc_upto(av,c0), v);
1762 : }
1763 :
1764 : /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1765 : GEN
1766 1821786 : RgX_mul_normalized(GEN A, long a, GEN B, long b)
1767 : {
1768 1821786 : GEN z = RgX_mul(A, B);
1769 1821775 : if (a < b)
1770 10330 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1771 1811445 : else if (a > b)
1772 1207069 : z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1773 : else
1774 604376 : z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1775 1821782 : return z;
1776 : }
1777 :
1778 : GEN
1779 32815118 : RgX_mul_i(GEN x, GEN y)
1780 : {
1781 32815118 : GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1782 32814792 : setvarn(z, varn(x)); return z;
1783 : }
1784 :
1785 : GEN
1786 33236 : RgX_sqr_i(GEN x)
1787 : {
1788 33236 : GEN z = RgX_sqrspec(x+2, lgpol(x));
1789 33236 : setvarn(z,varn(x)); return z;
1790 : }
1791 :
1792 : /*******************************************************************/
1793 : /* */
1794 : /* DIVISION */
1795 : /* */
1796 : /*******************************************************************/
1797 : GEN
1798 5238144 : RgX_Rg_divexact(GEN x, GEN y) {
1799 5238144 : long i, lx = lg(x);
1800 : GEN z;
1801 5238144 : if (lx == 2) return gcopy(x);
1802 5207860 : switch(typ(y))
1803 : {
1804 5097129 : case t_INT:
1805 5097129 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1806 5031753 : break;
1807 5257 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1808 : }
1809 5137227 : z = cgetg(lx, t_POL); z[1] = x[1];
1810 29086823 : for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1811 5137018 : return z;
1812 : }
1813 : /* cf RgX_coeff_simplify */
1814 : INLINE GEN
1815 36967866 : RgX_coeff_div(GEN a, GEN b)
1816 : {
1817 36967866 : pari_sp av = avma;
1818 36967866 : GEN s = gdiv(a, b);
1819 36967862 : if (typ(s) != t_POL) return s;
1820 19885097 : switch(lg(s))
1821 : {
1822 7265418 : case 2: return gc_const(av, gen_0);
1823 9559391 : case 3: return gc_GEN(av, gel(s,2));
1824 3060288 : default: return s;
1825 : }
1826 : }
1827 : GEN
1828 29637058 : RgX_Rg_div(GEN x, GEN y) {
1829 29637058 : if (lg(x) == 2) return gcopy(x);
1830 29559417 : switch(typ(y))
1831 : {
1832 20341874 : case t_INT:
1833 20341874 : if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1834 4219252 : break;
1835 5740 : case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1836 : }
1837 50398892 : pari_APPLY_pol(RgX_coeff_div(gel(x,i), y));
1838 : }
1839 : GEN
1840 39711 : RgX_normalize(GEN x)
1841 : {
1842 39711 : GEN z, d = NULL;
1843 39711 : long i, n = lg(x)-1;
1844 39711 : for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1845 39711 : if (i == 1) return pol_0(varn(x));
1846 39711 : if (i == n && isint1(d)) return x;
1847 18200 : n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1848 32571 : for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1849 18200 : gel(z,n) = Rg_get_1(d); return z;
1850 : }
1851 : GEN
1852 10948 : RgX_divs(GEN x, long y) { pari_APPLY_pol(gdivgs(gel(x,i),y)); }
1853 : GEN
1854 265216 : RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1855 : {
1856 265216 : long l = lg(a), i;
1857 : GEN a0, z0, z;
1858 :
1859 265216 : if (l <= 3)
1860 : {
1861 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1862 0 : return pol_0(varn(a));
1863 : }
1864 265216 : z = cgetg(l-1, t_POL);
1865 265221 : z[1] = a[1];
1866 265221 : a0 = a + l-1;
1867 265221 : z0 = z + l-2; *z0 = *a0--;
1868 3025195 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1869 : {
1870 2759985 : GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1871 2759974 : gel(z0,0) = t;
1872 : }
1873 265210 : if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1874 265210 : return z;
1875 : }
1876 : /* Polynomial division x / y:
1877 : * if pr = ONLY_REM return remainder, otherwise return quotient
1878 : * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1879 : * if pr != NULL set *pr to remainder, as the last object on stack */
1880 : /* assume, typ(x) = typ(y) = t_POL, same variable */
1881 : static GEN
1882 23448329 : RgX_divrem_i(GEN x, GEN y, GEN *pr)
1883 : {
1884 : pari_sp avy, av, av1;
1885 : long dx,dy,dz,i,j,sx,lr;
1886 : GEN z,p1,p2,rem,y_lead,mod,p;
1887 : GEN (*f)(GEN,GEN);
1888 :
1889 23448329 : if (!signe(y)) pari_err_INV("RgX_divrem",y);
1890 :
1891 23448329 : dy = degpol(y);
1892 23448322 : y_lead = gel(y,dy+2);
1893 23448322 : if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1894 : {
1895 0 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1896 0 : for (dy--; dy>=0; dy--)
1897 : {
1898 0 : y_lead = gel(y,dy+2);
1899 0 : if (!gequal0(y_lead)) break;
1900 : }
1901 : }
1902 23448319 : if (!dy) /* y is constant */
1903 : {
1904 6816 : if (pr == ONLY_REM) return pol_0(varn(x));
1905 6809 : z = RgX_Rg_div(x, y_lead);
1906 6809 : if (pr == ONLY_DIVIDES) return z;
1907 2595 : if (pr) *pr = pol_0(varn(x));
1908 2595 : return z;
1909 : }
1910 23441503 : dx = degpol(x);
1911 23441472 : if (dx < dy)
1912 : {
1913 3849872 : if (pr == ONLY_REM) return RgX_copy(x);
1914 346211 : if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1915 346190 : z = pol_0(varn(x));
1916 346190 : if (pr) *pr = RgX_copy(x);
1917 346190 : return z;
1918 : }
1919 :
1920 : /* x,y in R[X], y non constant */
1921 19591600 : av = avma;
1922 19591600 : p = NULL;
1923 19591600 : if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1924 : {
1925 127407 : z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1926 127407 : if (!z) return gc_NULL(av);
1927 127407 : z = FpX_to_mod(z, p);
1928 127407 : if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1929 68677 : return gc_upto(av, z);
1930 58730 : *pr = FpX_to_mod(*pr, p);
1931 58730 : return gc_all(av, 2, &z, pr);
1932 : }
1933 19464290 : switch(typ(y_lead))
1934 : {
1935 0 : case t_REAL:
1936 0 : y_lead = ginv(y_lead);
1937 0 : f = gmul; mod = NULL;
1938 0 : break;
1939 2795 : case t_INTMOD:
1940 2795 : case t_POLMOD: y_lead = ginv(y_lead);
1941 2795 : f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1942 2795 : break;
1943 19461495 : default: if (gequal1(y_lead)) y_lead = NULL;
1944 19461486 : f = gdiv; mod = NULL;
1945 : }
1946 :
1947 19464281 : if (y_lead == NULL)
1948 17503005 : p2 = gel(x,dx+2);
1949 : else {
1950 : for(;;) {
1951 1961276 : p2 = f(gel(x,dx+2),y_lead);
1952 1961277 : p2 = simplify_shallow(p2);
1953 1961277 : if (!isexactzero(p2) || (--dx < 0)) break;
1954 : }
1955 1961277 : if (dx < dy) /* leading coeff of x was in fact zero */
1956 : {
1957 0 : if (pr == ONLY_DIVIDES) {
1958 0 : set_avma(av);
1959 0 : return (dx < 0)? pol_0(varn(x)) : NULL;
1960 : }
1961 0 : if (pr == ONLY_REM)
1962 : {
1963 0 : if (dx < 0)
1964 0 : return gc_GEN(av, scalarpol(p2, varn(x)));
1965 : else
1966 : {
1967 : GEN t;
1968 0 : set_avma(av);
1969 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1970 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1971 0 : return t;
1972 : }
1973 : }
1974 0 : if (pr) /* cf ONLY_REM above */
1975 : {
1976 0 : if (dx < 0)
1977 : {
1978 0 : p2 = gclone(p2);
1979 0 : set_avma(av);
1980 0 : z = pol_0(varn(x));
1981 0 : x = scalarpol(p2, varn(x));
1982 0 : gunclone(p2);
1983 : }
1984 : else
1985 : {
1986 : GEN t;
1987 0 : set_avma(av);
1988 0 : z = pol_0(varn(x));
1989 0 : t = cgetg(dx + 3, t_POL); t[1] = x[1];
1990 0 : for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1991 0 : x = t;
1992 : }
1993 0 : *pr = x;
1994 : }
1995 : else
1996 : {
1997 0 : set_avma(av);
1998 0 : z = pol_0(varn(x));
1999 : }
2000 0 : return z;
2001 : }
2002 : }
2003 : /* dx >= dy */
2004 19464282 : avy = avma;
2005 19464282 : dz = dx-dy;
2006 19464282 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2007 19464252 : x += 2;
2008 19464252 : z += 2;
2009 19464252 : y += 2;
2010 19464252 : gel(z,dz) = gcopy(p2);
2011 :
2012 53752941 : for (i=dx-1; i>=dy; i--)
2013 : {
2014 34288825 : av1=avma; p1=gel(x,i);
2015 1139519874 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2016 34283990 : if (y_lead) p1 = simplify(f(p1,y_lead));
2017 :
2018 34283990 : if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
2019 : else
2020 24296515 : p1 = avma==av1? gcopy(p1): gc_upto(av1,p1);
2021 34288285 : gel(z,i-dy) = p1;
2022 : }
2023 19464116 : if (!pr) return gc_upto(av,z-2);
2024 :
2025 12539072 : rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
2026 12894835 : for (sx=0; ; i--)
2027 : {
2028 12894835 : p1 = gel(x,i);
2029 : /* we always enter this loop at least once */
2030 32177310 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2031 12893746 : if (mod && avma==av1) p1 = gmul(p1,mod);
2032 12893746 : if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
2033 3575980 : if (!isexactzero(p1)) break;
2034 3425623 : if (!i) break;
2035 355742 : set_avma(av1);
2036 : }
2037 12538455 : if (pr == ONLY_DIVIDES)
2038 : {
2039 8078 : if (sx) return gc_NULL(av);
2040 8057 : set_avma((pari_sp)rem); return gc_upto(av,z-2);
2041 : }
2042 12530377 : lr=i+3; rem -= lr;
2043 12530377 : if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
2044 12421081 : else p1 = gc_upto((pari_sp)rem,p1);
2045 12530967 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2046 12530967 : rem[1] = z[-1];
2047 12530967 : rem += 2;
2048 12530967 : gel(rem,i) = p1;
2049 17249296 : for (i--; i>=0; i--)
2050 : {
2051 4718369 : av1=avma; p1 = gel(x,i);
2052 13592933 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2053 4717868 : if (mod && avma==av1) p1 = gmul(p1,mod);
2054 4718090 : gel(rem,i) = avma==av1? gcopy(p1):gc_upto(av1,p1);
2055 : }
2056 12530927 : rem -= 2;
2057 12530927 : if (!sx) (void)normalizepol_lg(rem, lr);
2058 12531157 : if (pr == ONLY_REM) return gc_upto(av,rem);
2059 6768141 : z -= 2; *pr = rem; return gc_all_unsafe(av,avy,2,&z,pr);
2060 : }
2061 :
2062 : GEN
2063 14181721 : RgX_divrem(GEN x, GEN y, GEN *pr)
2064 : {
2065 14181721 : if (pr == ONLY_REM) return RgX_rem(x, y);
2066 14181721 : return RgX_divrem_i(x, y, pr);
2067 : }
2068 :
2069 : /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
2070 : GEN
2071 156752 : RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
2072 : {
2073 156752 : long vx = varn(x), dx = degpol(x), dy = degpol(y), dz, i, j, sx, lr;
2074 : pari_sp av0, av;
2075 : GEN z, p1, rem, lead;
2076 :
2077 156752 : if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
2078 156752 : if (dx < dy)
2079 : {
2080 41644 : if (pr)
2081 : {
2082 41616 : av0 = avma; x = RgXQX_red(x, T);
2083 41616 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
2084 41609 : if (pr == ONLY_REM) return x;
2085 0 : *pr = x;
2086 : }
2087 28 : return pol_0(vx);
2088 : }
2089 115108 : lead = leading_coeff(y);
2090 115108 : if (!dy) /* y is constant */
2091 : {
2092 602 : if (pr && pr != ONLY_DIVIDES)
2093 : {
2094 0 : if (pr == ONLY_REM) return pol_0(vx);
2095 0 : *pr = pol_0(vx);
2096 : }
2097 602 : if (gequal1(lead)) return RgX_copy(x);
2098 0 : av0 = avma; x = gmul(x, ginvmod(lead,T));
2099 0 : return gc_upto(av0, RgXQX_red(x,T));
2100 : }
2101 114506 : av0 = avma; dz = dx-dy;
2102 114506 : lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
2103 114506 : set_avma(av0);
2104 114506 : z = cgetg(dz+3,t_POL); z[1] = x[1];
2105 114506 : x += 2; y += 2; z += 2;
2106 :
2107 114506 : p1 = gel(x,dx); av = avma;
2108 114506 : gel(z,dz) = lead? gc_upto(av, grem(gmul(p1,lead), T)): gcopy(p1);
2109 581323 : for (i=dx-1; i>=dy; i--)
2110 : {
2111 466817 : av = avma; p1 = gel(x,i);
2112 2416804 : for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2113 466781 : if (lead) p1 = gmul(grem(p1, T), lead);
2114 466786 : gel(z,i-dy) = gc_upto(av, grem(p1, T));
2115 : }
2116 114506 : if (!pr) { guncloneNULL(lead); return z-2; }
2117 :
2118 113232 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
2119 191267 : for (sx=0; ; i--)
2120 : {
2121 191267 : p1 = gel(x,i);
2122 661291 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2123 191265 : p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
2124 105538 : if (!i) break;
2125 78035 : set_avma(av);
2126 : }
2127 113232 : if (pr == ONLY_DIVIDES)
2128 : {
2129 27038 : guncloneNULL(lead);
2130 27038 : if (sx) return gc_NULL(av0);
2131 24479 : return gc_const((pari_sp)rem, z-2);
2132 : }
2133 86194 : lr=i+3; rem -= lr; av = (pari_sp)rem;
2134 86194 : rem[0] = evaltyp(t_POL) | _evallg(lr);
2135 86194 : rem[1] = z[-1];
2136 86194 : rem += 2; gel(rem,i) = gc_upto(av, p1);
2137 187610 : for (i--; i>=0; i--)
2138 : {
2139 101416 : av = avma; p1 = gel(x,i);
2140 337537 : for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
2141 101416 : gel(rem,i) = gc_upto(av, grem(p1, T));
2142 : }
2143 86194 : rem -= 2;
2144 86194 : guncloneNULL(lead);
2145 86194 : if (!sx) (void)normalizepol_lg(rem, lr);
2146 86194 : if (pr == ONLY_REM) return gc_upto(av0,rem);
2147 168 : *pr = rem; return z-2;
2148 : }
2149 :
2150 : /*******************************************************************/
2151 : /* */
2152 : /* PSEUDO-DIVISION */
2153 : /* */
2154 : /*******************************************************************/
2155 : INLINE GEN
2156 1050958 : rem(GEN c, GEN T)
2157 : {
2158 1050958 : if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
2159 1050958 : return c;
2160 : }
2161 :
2162 : /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
2163 : int
2164 17720 : ZXQX_dvd(GEN x, GEN y, GEN T)
2165 : {
2166 : long dx, dy, i, T_ismonic;
2167 17720 : pari_sp av = avma, av2;
2168 : GEN y_lead;
2169 :
2170 17720 : if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
2171 17720 : dy = degpol(y); y_lead = gel(y,dy+2);
2172 17720 : if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
2173 : /* if monic, no point in using pseudo-division */
2174 17720 : if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
2175 14738 : T_ismonic = gequal1(leading_coeff(T));
2176 14738 : dx = degpol(x);
2177 14738 : if (dx < dy) return !signe(x);
2178 14738 : (void)new_chunk(2);
2179 14738 : x = RgX_recip_i(x)+2;
2180 14738 : y = RgX_recip_i(y)+2;
2181 : /* pay attention to sparse divisors */
2182 29977 : for (i = 1; i <= dy; i++)
2183 15239 : if (!signe(gel(y,i))) gel(y,i) = NULL;
2184 14738 : av2 = avma;
2185 : for (;;)
2186 72953 : {
2187 87691 : GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
2188 87691 : x0 = gneg(x0);
2189 87691 : m = gcdii(cx, y0);
2190 87691 : if (!equali1(m))
2191 : {
2192 85189 : x0 = gdiv(x0, m);
2193 85189 : y0 = diviiexact(y0, m);
2194 85189 : if (equali1(y0)) y0 = NULL;
2195 : }
2196 179155 : for (i=1; i<=dy; i++)
2197 : {
2198 91464 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2199 91464 : if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
2200 91464 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2201 91464 : gel(x,i) = c;
2202 : }
2203 687972 : for ( ; i<=dx; i++)
2204 : {
2205 600281 : GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2206 600281 : if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2207 600281 : gel(x,i) = c;
2208 : }
2209 102167 : do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2210 87691 : if (dx < dy) break;
2211 72953 : if (gc_needed(av2,1))
2212 : {
2213 28 : if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2214 28 : gc_slice(av2,x,dx+1);
2215 : }
2216 : }
2217 14738 : return gc_bool(av, dx < 0);
2218 : }
2219 :
2220 : /* T either NULL or a t_POL. */
2221 : GEN
2222 138802 : RgXQX_pseudorem(GEN x, GEN y, GEN T)
2223 : {
2224 138802 : long vx = varn(x), dx, dy, dz, i, lx, p;
2225 138802 : pari_sp av = avma, av2;
2226 : GEN y_lead;
2227 :
2228 138802 : if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2229 138802 : dy = degpol(y); y_lead = gel(y,dy+2);
2230 : /* if monic, no point in using pseudo-division */
2231 138802 : if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2232 125088 : dx = degpol(x);
2233 125088 : if (dx < dy) return RgX_copy(x);
2234 125081 : (void)new_chunk(2);
2235 125081 : x = RgX_recip_i(x)+2;
2236 125081 : y = RgX_recip_i(y)+2;
2237 : /* pay attention to sparse divisors */
2238 535552 : for (i = 1; i <= dy; i++)
2239 410471 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2240 125081 : dz = dx-dy; p = dz+1;
2241 125081 : av2 = avma;
2242 : for (;;)
2243 : {
2244 247146 : gel(x,0) = gneg(gel(x,0)); p--;
2245 1061830 : for (i=1; i<=dy; i++)
2246 : {
2247 814686 : GEN c = gmul(y_lead, gel(x,i));
2248 814676 : if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2249 814685 : gel(x,i) = rem(c, T);
2250 : }
2251 373887 : for ( ; i<=dx; i++)
2252 : {
2253 126741 : GEN c = gmul(y_lead, gel(x,i));
2254 126742 : gel(x,i) = rem(c, T);
2255 : }
2256 253346 : do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2257 247146 : if (dx < dy) break;
2258 122065 : if (gc_needed(av2,1))
2259 : {
2260 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2261 0 : gc_slice(av2,x,dx+1);
2262 : }
2263 : }
2264 125081 : if (dx < 0) return pol_0(vx);
2265 124920 : lx = dx+3; x -= 2;
2266 124920 : x[0] = evaltyp(t_POL) | _evallg(lx);
2267 124920 : x[1] = evalsigne(1) | evalvarn(vx);
2268 124920 : x = RgX_recip_i(x);
2269 124920 : if (p)
2270 : { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2271 3884 : GEN t = y_lead;
2272 3884 : if (T && typ(t) == t_POL && varn(t) == varn(T))
2273 0 : t = RgXQ_powu(t, p, T);
2274 : else
2275 3884 : t = gpowgs(t, p);
2276 11803 : for (i=2; i<lx; i++)
2277 : {
2278 7919 : GEN c = gmul(gel(x,i), t);
2279 7919 : gel(x,i) = rem(c,T);
2280 : }
2281 3884 : if (!T) return gc_upto(av, x);
2282 : }
2283 121036 : return gc_GEN(av, x);
2284 : }
2285 :
2286 : GEN
2287 138802 : RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2288 :
2289 : /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2290 : GEN
2291 12077 : RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2292 : {
2293 12077 : long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2294 12077 : pari_sp av = avma, av2;
2295 : GEN z, r, ypow, y_lead;
2296 :
2297 12077 : if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2298 12077 : dy = degpol(y); y_lead = gel(y,dy+2);
2299 12077 : if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2300 7561 : dx = degpol(x);
2301 7561 : if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2302 7561 : if (dx == dy)
2303 : {
2304 98 : GEN x_lead = gel(x,lg(x)-1);
2305 98 : x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2306 98 : y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2307 98 : r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2308 98 : *ptr = gc_upto(av, r); return scalarpol(x_lead, vx);
2309 : }
2310 7463 : (void)new_chunk(2);
2311 7463 : x = RgX_recip_i(x)+2;
2312 7463 : y = RgX_recip_i(y)+2;
2313 : /* pay attention to sparse divisors */
2314 39000 : for (i = 1; i <= dy; i++)
2315 31537 : if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2316 7463 : dz = dx-dy; p = dz+1;
2317 7463 : lz = dz+3;
2318 7463 : z = cgetg(lz, t_POL);
2319 7463 : z[1] = evalsigne(1) | evalvarn(vx);
2320 28640 : for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2321 7463 : ypow = new_chunk(dz+1);
2322 7463 : gel(ypow,0) = gen_1;
2323 7463 : gel(ypow,1) = y_lead;
2324 13714 : for (i=2; i<=dz; i++)
2325 : {
2326 6251 : GEN c = gmul(gel(ypow,i-1), y_lead);
2327 6251 : gel(ypow,i) = rem(c,T);
2328 : }
2329 7463 : av2 = avma;
2330 7463 : for (iz=2;;)
2331 : {
2332 15605 : p--;
2333 15605 : gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2334 69915 : for (i=1; i<=dy; i++)
2335 : {
2336 54310 : GEN c = gmul(y_lead, gel(x,i));
2337 54310 : if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2338 54310 : gel(x,i) = rem(c, T);
2339 : }
2340 41058 : for ( ; i<=dx; i++)
2341 : {
2342 25453 : GEN c = gmul(y_lead, gel(x,i));
2343 25453 : gel(x,i) = rem(c,T);
2344 : }
2345 15605 : x++; dx--;
2346 21177 : while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2347 15605 : if (dx < dy) break;
2348 8142 : if (gc_needed(av2,1))
2349 : {
2350 0 : GEN X = x-2;
2351 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2352 0 : X[0] = evaltyp(t_POL)|_evallg(dx+3); X[1] = z[1]; /* hack */
2353 0 : (void)gc_all(av2,2, &X, &z); x = X+2;
2354 : }
2355 : }
2356 14008 : while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2357 7463 : if (dx < 0)
2358 182 : x = pol_0(vx);
2359 : else
2360 : {
2361 7281 : lx = dx+3; x -= 2;
2362 7281 : x[0] = evaltyp(t_POL) | _evallg(lx);
2363 7281 : x[1] = evalsigne(1) | evalvarn(vx);
2364 7281 : x = RgX_recip_i(x);
2365 : }
2366 7463 : z = RgX_recip_i(z);
2367 7463 : r = x;
2368 7463 : if (p)
2369 : {
2370 3339 : GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2371 3339 : if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2372 : }
2373 7463 : *ptr = r; return gc_all(av, 2, &z, ptr);
2374 : }
2375 : GEN
2376 11832 : RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2377 11832 : { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2378 :
2379 : GEN
2380 0 : RgXQX_mul(GEN x, GEN y, GEN T)
2381 0 : { return RgXQX_red(RgX_mul(x,y), T); }
2382 :
2383 : /* cf RgX_coeff_simplify */
2384 : INLINE GEN
2385 728329183 : RgX_coeff_mul(GEN a, GEN b)
2386 : {
2387 728329183 : pari_sp av = avma;
2388 728329183 : GEN s = gmul(a, b);
2389 728277937 : if (typ(s) != t_POL) return s;
2390 129342050 : switch(lg(s))
2391 : {
2392 72205518 : case 2: return gc_const(av, gen_0);
2393 730568 : case 3: return gc_GEN(av, gel(s,2));
2394 56405964 : default: return s;
2395 : }
2396 : }
2397 : GEN
2398 965050855 : RgX_Rg_mul(GEN x, GEN y) { pari_APPLY_pol(RgX_coeff_mul(y, gel(x,i))); }
2399 : GEN
2400 405218 : RgX_mul2n(GEN x, long n) { pari_APPLY_pol(gmul2n(gel(x,i), n)); }
2401 : GEN
2402 26285 : RgX_muls(GEN x, long y) { pari_APPLY_pol(gmulsg(y, gel(x,i))); }
2403 : GEN
2404 35 : RgXQX_RgXQ_mul(GEN x, GEN y, GEN T) { return RgXQX_red(RgX_Rg_mul(x,y), T); }
2405 : GEN
2406 168 : RgXQV_RgXQ_mul(GEN v, GEN x, GEN T) { return RgXQV_red(RgV_Rg_mul(v,x), T); }
2407 :
2408 : GEN
2409 0 : RgXQX_sqr(GEN x, GEN T) { return RgXQX_red(RgX_sqr(x), T); }
2410 :
2411 : GEN
2412 0 : RgXQX_powers(GEN P, long n, GEN T)
2413 : {
2414 0 : GEN v = cgetg(n+2, t_VEC);
2415 : long i;
2416 0 : gel(v, 1) = pol_1(varn(T));
2417 0 : if (n==0) return v;
2418 0 : gel(v, 2) = gcopy(P);
2419 0 : for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2420 0 : return v;
2421 : }
2422 :
2423 : static GEN
2424 623235 : _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2425 : static GEN
2426 0 : _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2427 : static GEN
2428 74805 : _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2429 : static GEN
2430 88307 : _pow(void *data, GEN x, GEN n) { return RgXQ_pow(x, n, (GEN)data); }
2431 : static GEN
2432 280959 : _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2433 : static GEN
2434 1040369 : _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2435 : static GEN
2436 1019354 : _one(void *data) { return pol_1(varn((GEN)data)); }
2437 : static GEN
2438 1232 : _zero(void *data) { return pol_0(varn((GEN)data)); }
2439 : static GEN
2440 724098 : _red(void *data, GEN x) { (void)data; return gcopy(x); }
2441 :
2442 : static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2443 : _mul, _sqr, _one, _zero };
2444 :
2445 : GEN
2446 0 : RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2447 : {
2448 0 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2449 : }
2450 :
2451 : GEN
2452 417394 : RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2453 : {
2454 417394 : int use_sqr = 2*degpol(x) >= degpol(T);
2455 417393 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2456 : }
2457 :
2458 : /* mod X^n */
2459 : struct modXn {
2460 : long v; /* varn(X) */
2461 : long n;
2462 : } ;
2463 : static GEN
2464 11893 : _sqrXn(void *data, GEN x) {
2465 11893 : struct modXn *S = (struct modXn*)data;
2466 11893 : return RgXn_sqr(x, S->n);
2467 : }
2468 : static GEN
2469 4528 : _mulXn(void *data, GEN x, GEN y) {
2470 4528 : struct modXn *S = (struct modXn*)data;
2471 4528 : return RgXn_mul(x,y, S->n);
2472 : }
2473 : static GEN
2474 1939 : _oneXn(void *data) {
2475 1939 : struct modXn *S = (struct modXn*)data;
2476 1939 : return pol_1(S->v);
2477 : }
2478 : static GEN
2479 0 : _zeroXn(void *data) {
2480 0 : struct modXn *S = (struct modXn*)data;
2481 0 : return pol_0(S->v);
2482 : }
2483 : static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2484 : _oneXn, _zeroXn };
2485 :
2486 : GEN
2487 483 : RgXn_powers(GEN x, long m, long n)
2488 : {
2489 483 : long d = degpol(x);
2490 483 : int use_sqr = (d<<1) >= n;
2491 : struct modXn S;
2492 483 : S.v = varn(x); S.n = n;
2493 483 : return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2494 : }
2495 :
2496 : GEN
2497 2286 : RgXn_powu_i(GEN x, ulong m, long n)
2498 : {
2499 : struct modXn S;
2500 : long v;
2501 2286 : if (n == 0) return x;
2502 2286 : v = RgX_valrem(x, &x);
2503 2286 : if (v) { n -= m * v; if (n <= 0) return pol_0(varn(x)); }
2504 2265 : S.v = varn(x); S.n = n;
2505 2265 : x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2506 2265 : if (v) x = RgX_shift_shallow(x, m * v);
2507 2265 : return x;
2508 : }
2509 : GEN
2510 0 : RgXn_powu(GEN x, ulong m, long n)
2511 : {
2512 : pari_sp av;
2513 0 : if (n == 0) return gcopy(x);
2514 0 : av = avma; return gc_GEN(av, RgXn_powu_i(x, m, n));
2515 : }
2516 :
2517 : GEN
2518 966 : RgX_RgXnV_eval(GEN Q, GEN x, long n)
2519 : {
2520 : struct modXn S;
2521 966 : S.v = varn(gel(x,2)); S.n = n;
2522 966 : return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2523 : }
2524 :
2525 : GEN
2526 0 : RgX_RgXn_eval(GEN Q, GEN x, long n)
2527 : {
2528 0 : int use_sqr = 2*degpol(x) >= n;
2529 : struct modXn S;
2530 0 : S.v = varn(x); S.n = n;
2531 0 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2532 : }
2533 :
2534 : /* Q(x) mod t^n, x in R[t], n >= 1 */
2535 : GEN
2536 5593 : RgXn_eval(GEN Q, GEN x, long n)
2537 : {
2538 5593 : long d = degpol(x);
2539 : int use_sqr;
2540 : struct modXn S;
2541 5593 : if (d == 1 && isrationalzero(gel(x,2)))
2542 : {
2543 5586 : GEN y = RgX_unscale(Q, gel(x,3));
2544 5586 : setvarn(y, varn(x)); return y;
2545 : }
2546 7 : S.v = varn(x);
2547 7 : S.n = n;
2548 7 : use_sqr = (d<<1) >= n;
2549 7 : return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2550 : }
2551 :
2552 : /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2553 : static GEN
2554 2848403 : RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2555 : {
2556 2848403 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2557 2848403 : return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2558 : }
2559 :
2560 : /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2561 : static GEN
2562 14 : RgXn_sqrhigh(GEN f, long n2, long n)
2563 : {
2564 14 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2565 14 : return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2566 : }
2567 :
2568 : static GEN
2569 2371582 : RgXn_div_gen(GEN g, GEN f, long e)
2570 : {
2571 : pari_sp av;
2572 : ulong mask;
2573 : GEN W, a;
2574 2371582 : long v = varn(f), n = 1;
2575 :
2576 2371582 : if (!signe(f)) pari_err_INV("RgXn_inv",f);
2577 2371575 : a = ginv(gel(f,2));
2578 2371567 : if (e == 1 && !g) return scalarpol(a, v);
2579 2365918 : else if (e == 2 && !g)
2580 : {
2581 : GEN b;
2582 929121 : if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2583 333543 : b = gneg(b);
2584 333541 : if (!gequal1(a)) b = gmul(b, gsqr(a));
2585 333541 : return deg1pol(b, a, v);
2586 : }
2587 1436797 : av = avma;
2588 1436797 : W = scalarpol_shallow(a,v);
2589 1436796 : mask = quadratic_prec_mask(e);
2590 4279242 : while (mask > 1)
2591 : {
2592 : GEN u, fr;
2593 2842446 : long n2 = n;
2594 2842446 : n<<=1; if (mask & 1) n--;
2595 2842446 : mask >>= 1;
2596 2842446 : fr = RgXn_red_shallow(f, n);
2597 2842446 : if (mask>1 || !g)
2598 : {
2599 1476123 : u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2600 1476123 : W = RgX_sub(W, RgX_shift_shallow(u, n2));
2601 : }
2602 : else
2603 : {
2604 1366323 : GEN y = RgXn_mul(g, W, n), yt = RgXn_red_shallow(y, n-n2);
2605 1366323 : u = RgXn_mul(yt, RgXn_mulhigh(fr, W, n2, n), n-n2);
2606 1366323 : W = RgX_sub(y, RgX_shift_shallow(u, n2));
2607 : }
2608 2842446 : if (gc_needed(av,2))
2609 : {
2610 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2611 0 : W = gc_upto(av, W);
2612 : }
2613 : }
2614 1436796 : return W;
2615 : }
2616 :
2617 : static GEN
2618 98 : RgXn_div_FpX(GEN x, GEN y, long e, GEN p)
2619 : {
2620 : GEN r;
2621 98 : if (lgefint(p) == 3)
2622 : {
2623 98 : ulong pp = uel(p, 2);
2624 98 : if (pp == 2)
2625 7 : r = F2x_to_ZX(F2xn_div(RgX_to_F2x(x), RgX_to_F2x(y), e));
2626 : else
2627 91 : r = Flx_to_ZX_inplace(Flxn_div(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), e, pp));
2628 : }
2629 : else
2630 0 : r = FpXn_div(RgX_to_FpX(x, p), RgX_to_FpX(y, p), e, p);
2631 98 : return FpX_to_mod(r, p);
2632 : }
2633 :
2634 : static GEN
2635 0 : RgXn_div_FpXQX(GEN x, GEN y, long n, GEN pol, GEN p)
2636 : {
2637 0 : GEN r, T = RgX_to_FpX(pol, p);
2638 0 : if (signe(T) == 0) pari_err_OP("/", x, y);
2639 0 : r = FpXQXn_div(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), n, T, p);
2640 0 : return FpXQX_to_mod(r, T, p);
2641 : }
2642 :
2643 : static GEN
2644 91 : RgXn_inv_FpX(GEN x, long e, GEN p)
2645 : {
2646 : GEN r;
2647 91 : if (lgefint(p) == 3)
2648 : {
2649 91 : ulong pp = uel(p, 2);
2650 91 : if (pp == 2)
2651 28 : r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2652 : else
2653 63 : r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2654 : }
2655 : else
2656 0 : r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2657 91 : return FpX_to_mod(r, p);
2658 : }
2659 :
2660 : static GEN
2661 0 : RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2662 : {
2663 0 : GEN r, T = RgX_to_FpX(pol, p);
2664 0 : if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2665 0 : r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2666 0 : return FpXQX_to_mod(r, T, p);
2667 : }
2668 :
2669 : static GEN
2670 1005341 : RgXn_inv_fast(GEN x, long e)
2671 : {
2672 : GEN p, pol;
2673 : long pa;
2674 1005341 : long t = RgX_type(x,&p,&pol,&pa);
2675 1005350 : switch(t)
2676 : {
2677 91 : case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2678 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2679 0 : return RgXn_inv_FpXQX(x, e, pol, p);
2680 1005259 : default: return NULL;
2681 : }
2682 : }
2683 :
2684 : static GEN
2685 1366421 : RgXn_div_fast(GEN x, GEN y, long e)
2686 : {
2687 : GEN p, pol;
2688 : long pa;
2689 1366421 : long t = RgX_type2(x,y,&p,&pol,&pa);
2690 1366421 : switch(t)
2691 : {
2692 98 : case t_INTMOD: return RgXn_div_FpX(x, y, e, p);
2693 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2694 0 : return RgXn_div_FpXQX(x, y, e, pol, p);
2695 1366323 : default: return NULL;
2696 : }
2697 : }
2698 :
2699 : GEN
2700 1366421 : RgXn_div_i(GEN g, GEN f, long e)
2701 : {
2702 1366421 : GEN h = RgXn_div_fast(g, f, e);
2703 1366421 : if (h) return h;
2704 1366323 : return RgXn_div_gen(g, f, e);
2705 : }
2706 :
2707 : GEN
2708 562607 : RgXn_div(GEN g, GEN f, long e)
2709 : {
2710 562607 : pari_sp av = avma;
2711 562607 : return gc_upto(av, RgXn_div_i(g, f, e));
2712 : }
2713 :
2714 : GEN
2715 1005341 : RgXn_inv_i(GEN f, long e)
2716 : {
2717 1005341 : GEN h = RgXn_inv_fast(f, e);
2718 1005350 : if (h) return h;
2719 1005259 : return RgXn_div_gen(NULL, f, e);
2720 : }
2721 :
2722 : GEN
2723 808301 : RgXn_inv(GEN f, long e)
2724 : {
2725 808301 : pari_sp av = avma;
2726 808301 : return gc_upto(av, RgXn_inv_i(f, e));
2727 : }
2728 :
2729 : /* intformal(x^n*S) / x^(n+1) */
2730 : static GEN
2731 56376 : RgX_integXn(GEN x, long n)
2732 114355 : { pari_APPLY_pol_normalized(gdivgs(gel(x,i), n+i-1)); }
2733 :
2734 : GEN
2735 52920 : RgXn_expint(GEN h, long e)
2736 : {
2737 52920 : pari_sp av = avma, av2;
2738 52920 : long v = varn(h), n;
2739 52920 : GEN f = pol_1(v), g;
2740 : ulong mask;
2741 :
2742 52919 : if (!signe(h)) return f;
2743 50427 : g = pol_1(v);
2744 50427 : n = 1; mask = quadratic_prec_mask(e);
2745 50427 : av2 = avma;
2746 56377 : for (;mask>1;)
2747 : {
2748 : GEN u, w;
2749 56377 : long n2 = n;
2750 56377 : n<<=1; if (mask & 1) n--;
2751 56377 : mask >>= 1;
2752 56377 : u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2753 56375 : u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2754 56376 : w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2755 56377 : f = RgX_add(f, RgX_shift_shallow(w, n2));
2756 56376 : if (mask<=1) break;
2757 5950 : u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2758 5950 : g = RgX_sub(g, RgX_shift_shallow(u, n2));
2759 5950 : if (gc_needed(av2,2))
2760 : {
2761 0 : if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2762 0 : (void)gc_all(av2, 2, &f, &g);
2763 : }
2764 : }
2765 50426 : return gc_upto(av, f);
2766 : }
2767 :
2768 : GEN
2769 0 : RgXn_exp(GEN h, long e)
2770 : {
2771 0 : long d = degpol(h);
2772 0 : if (d < 0) return pol_1(varn(h));
2773 0 : if (!d || !gequal0(gel(h,2)))
2774 0 : pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2775 0 : return RgXn_expint(RgX_deriv(h), e);
2776 : }
2777 :
2778 : GEN
2779 154 : RgXn_reverse(GEN f, long e)
2780 : {
2781 154 : pari_sp av = avma, av2;
2782 : ulong mask;
2783 : GEN fi, a, df, W, an;
2784 154 : long v = varn(f), n=1;
2785 154 : if (degpol(f)<1 || !gequal0(gel(f,2)))
2786 0 : pari_err_INV("serreverse",f);
2787 154 : fi = ginv(gel(f,3));
2788 154 : a = deg1pol_shallow(fi,gen_0,v);
2789 154 : if (e <= 2) return gc_GEN(av, a);
2790 133 : W = scalarpol(fi,v);
2791 133 : df = RgX_deriv(f);
2792 133 : mask = quadratic_prec_mask(e);
2793 133 : av2 = avma;
2794 616 : for (;mask>1;)
2795 : {
2796 : GEN u, fa, fr;
2797 483 : long n2 = n, rt;
2798 483 : n<<=1; if (mask & 1) n--;
2799 483 : mask >>= 1;
2800 483 : fr = RgXn_red_shallow(f, n);
2801 483 : rt = brent_kung_optpow(degpol(fr), 4, 3);
2802 483 : an = RgXn_powers(a, rt, n);
2803 483 : if (n>1)
2804 : {
2805 483 : long n4 = (n2+1)>>1;
2806 483 : GEN dfr = RgXn_red_shallow(df, n2);
2807 483 : dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2808 483 : u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2809 483 : W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2810 : }
2811 483 : fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2812 483 : fa = RgX_shift(fa, -n2);
2813 483 : a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2814 483 : if (gc_needed(av2,2))
2815 : {
2816 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2817 0 : (void)gc_all(av2, 2, &a, &W);
2818 : }
2819 : }
2820 133 : return gc_upto(av, a);
2821 : }
2822 :
2823 : GEN
2824 7 : RgXn_sqrt(GEN h, long e)
2825 : {
2826 7 : pari_sp av = avma, av2;
2827 7 : long v = varn(h), n = 1;
2828 7 : GEN f = scalarpol(gen_1, v), df = f;
2829 7 : ulong mask = quadratic_prec_mask(e);
2830 7 : if (degpol(h)<0 || !gequal1(gel(h,2)))
2831 0 : pari_err_SQRTN("RgXn_sqrt",h);
2832 7 : av2 = avma;
2833 : while(1)
2834 7 : {
2835 14 : long n2 = n, m;
2836 : GEN g;
2837 14 : n<<=1; if (mask & 1) n--;
2838 14 : mask >>= 1;
2839 14 : m = n-n2;
2840 14 : g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2841 14 : f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2842 14 : if (mask==1) return gc_upto(av, f);
2843 7 : g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2844 7 : df = RgX_sub(df, RgX_shift_shallow(g, n2));
2845 7 : if (gc_needed(av2,2))
2846 : {
2847 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2848 0 : (void)gc_all(av2, 2, &f, &df);
2849 : }
2850 : }
2851 : }
2852 :
2853 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2854 : GEN
2855 120406 : RgXQ_powu(GEN x, ulong n, GEN T)
2856 : {
2857 120406 : pari_sp av = avma;
2858 :
2859 120406 : if (!n) return pol_1(varn(x));
2860 72159 : if (n == 1) return RgX_copy(x);
2861 23527 : x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2862 23527 : return gc_GEN(av, x);
2863 : }
2864 : /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2865 : GEN
2866 102160 : RgXQ_pow(GEN x, GEN n, GEN T)
2867 : {
2868 : pari_sp av;
2869 102160 : long s = signe(n);
2870 :
2871 102160 : if (!s) return pol_1(varn(x));
2872 102160 : if (is_pm1(n) == 1)
2873 88307 : return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2874 13853 : av = avma;
2875 13853 : if (s < 0) x = RgXQ_inv(x, T);
2876 13853 : x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2877 13853 : return gc_GEN(av, x);
2878 : }
2879 : static GEN
2880 200593 : _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2881 : static GEN
2882 111182 : _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2883 :
2884 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2885 : GEN
2886 12061 : ZXQ_powers(GEN x, long l, GEN T)
2887 : {
2888 12061 : int use_sqr = 2*degpol(x) >= degpol(T);
2889 12061 : return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2890 : }
2891 :
2892 : /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2893 : GEN
2894 168496 : ZXQ_powu(GEN x, ulong n, GEN T)
2895 : {
2896 168496 : pari_sp av = avma;
2897 :
2898 168496 : if (!n) return pol_1(varn(x));
2899 168496 : if (n == 1) return ZX_copy(x);
2900 113867 : x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2901 113875 : return gc_GEN(av, x);
2902 : }
2903 :
2904 : /* generates the list of powers of x of degree 0,1,2,...,l*/
2905 : GEN
2906 8603 : RgXQ_powers(GEN x, long l, GEN T)
2907 : {
2908 8603 : int use_sqr = 2*degpol(x) >= degpol(T);
2909 8603 : return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2910 : }
2911 :
2912 : GEN
2913 7118 : RgXQV_factorback(GEN L, GEN e, GEN T)
2914 : {
2915 7118 : return gen_factorback(L, e, (void*)T, &_mul, &_pow, &_one);
2916 : }
2917 :
2918 : /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2919 : GEN
2920 9506 : QXQ_powers(GEN a, long n, GEN T)
2921 : {
2922 : GEN den, v;
2923 9506 : if (!isint1(leading_coeff(T))) return RgXQ_powers(a, n, T);
2924 9478 : v = ZXQ_powers(Q_remove_denom(a, &den), n, T);
2925 : /* den*a integral; v[i+1] = (den*a)^i in K */
2926 9478 : if (den)
2927 : { /* restore denominators */
2928 5942 : GEN d = den;
2929 : long i;
2930 5942 : gel(v,2) = a;
2931 27679 : for (i=3; i<=n+1; i++) {
2932 21737 : d = mulii(d,den);
2933 21737 : gel(v,i) = RgX_Rg_div(gel(v,i), d);
2934 : }
2935 : }
2936 9478 : return v;
2937 : }
2938 :
2939 : static GEN
2940 3542 : do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2941 : {
2942 3542 : long l, i, m = 0;
2943 : GEN dz, z;
2944 3542 : GEN V = cgetg_copy(v, &l);
2945 12054 : for (i = imin; i < l; i++)
2946 : {
2947 8512 : GEN c = gel(v, i);
2948 8512 : if (typ(c) == t_POL) m = maxss(m, degpol(c));
2949 : }
2950 3542 : z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2951 3836 : for (i = 1; i < imin; i++) V[i] = v[i];
2952 12054 : for (i = imin; i < l; i++)
2953 : {
2954 8512 : GEN c = gel(v,i);
2955 8512 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2956 8512 : gel(V,i) = c;
2957 : }
2958 3542 : return V;
2959 : }
2960 : /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2961 : GEN
2962 3248 : QXV_QXQ_eval(GEN v, GEN a, GEN T)
2963 3248 : { return do_QXQ_eval(v, 1, a, T); }
2964 :
2965 : GEN
2966 294 : QXY_QXQ_evalx(GEN v, GEN a, GEN T)
2967 294 : { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2968 :
2969 : GEN
2970 2933 : RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2971 : {
2972 2933 : return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2973 : }
2974 :
2975 : GEN
2976 2659 : RgXQ_norm(GEN x, GEN T)
2977 : {
2978 : pari_sp av;
2979 2659 : long dx = degpol(x);
2980 : GEN L, y;
2981 2659 : if (degpol(T)==0) return gpowgs(x,0);
2982 2652 : av = avma; y = resultant(T, x);
2983 2652 : L = leading_coeff(T);
2984 2652 : if (gequal1(L) || !signe(x)) return y;
2985 0 : return gc_upto(av, gdiv(y, gpowgs(L, dx)));
2986 : }
2987 :
2988 : GEN
2989 476 : RgXQ_trace(GEN x, GEN T)
2990 : {
2991 476 : pari_sp av = avma;
2992 : GEN dT, z;
2993 : long n;
2994 476 : if (degpol(T)==0) return gmulgs(x,0);
2995 469 : dT = RgX_deriv(T); n = degpol(dT);
2996 469 : z = RgXQ_mul(x, dT, T);
2997 469 : if (degpol(z)<n) return gc_const(av, gen_0);
2998 420 : return gc_upto(av, gdiv(gel(z,2+n), gel(T,3+n)));
2999 : }
3000 :
3001 : GEN
3002 3021591 : RgX_blocks(GEN P, long n, long m)
3003 : {
3004 3021591 : GEN z = cgetg(m+1,t_VEC);
3005 3021591 : long i,j, k=2, l = lg(P);
3006 9258351 : for(i=1; i<=m; i++)
3007 : {
3008 6236760 : GEN zi = cgetg(n+2,t_POL);
3009 6236760 : zi[1] = P[1];
3010 6236760 : gel(z,i) = zi;
3011 18250978 : for(j=2; j<n+2; j++)
3012 12014218 : gel(zi, j) = k==l ? gen_0 : gel(P,k++);
3013 6236760 : zi = RgX_renormalize_lg(zi, n+2);
3014 : }
3015 3021591 : return z;
3016 : }
3017 :
3018 : /* write p(X) = e(X^2) + Xo(X^2), shallow function */
3019 : void
3020 49866199 : RgX_even_odd(GEN p, GEN *pe, GEN *po)
3021 : {
3022 49866199 : long n = degpol(p), v = varn(p), n0, n1, i;
3023 : GEN p0, p1;
3024 :
3025 49867179 : if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
3026 :
3027 49772635 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
3028 49772635 : p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
3029 49775368 : p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
3030 147084880 : for (i=0; i<n1; i++)
3031 : {
3032 97310186 : p0[2+i] = p[2+(i<<1)];
3033 97310186 : p1[2+i] = p[3+(i<<1)];
3034 : }
3035 49774694 : if (n1 != n0)
3036 18073213 : p0[2+i] = p[2+(i<<1)];
3037 49774694 : *pe = normalizepol(p0);
3038 49777442 : *po = normalizepol(p1);
3039 : }
3040 :
3041 : /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
3042 : GEN
3043 43631 : RgX_splitting(GEN p, long k)
3044 : {
3045 43631 : long n = degpol(p), v = varn(p), m, i, j, l;
3046 : GEN r;
3047 :
3048 43631 : m = n/k;
3049 43631 : r = cgetg(k+1,t_VEC);
3050 234059 : for(i=1; i<=k; i++)
3051 : {
3052 190428 : gel(r,i) = cgetg(m+3, t_POL);
3053 190428 : mael(r,i,1) = evalvarn(v)|evalsigne(1);
3054 : }
3055 571991 : for (j=1, i=0, l=2; i<=n; i++)
3056 : {
3057 528360 : gmael(r,j,l) = gel(p,2+i);
3058 528360 : if (j==k) { j=1; l++; } else j++;
3059 : }
3060 234059 : for(i=1; i<=k; i++)
3061 190428 : gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
3062 43631 : return r;
3063 : }
3064 :
3065 : /*******************************************************************/
3066 : /* */
3067 : /* Kronecker form */
3068 : /* */
3069 : /*******************************************************************/
3070 :
3071 : /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
3072 : * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
3073 : * normalized coefficients */
3074 : GEN
3075 180561 : Kronecker_to_mod(GEN z, GEN T)
3076 : {
3077 180561 : long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
3078 180561 : GEN x, t = cgetg(N,t_POL);
3079 180561 : t[1] = T[1];
3080 180561 : lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
3081 180561 : x[1] = z[1];
3082 180561 : T = RgX_copy(T);
3083 1375389 : for (i=2; i<lx+2; i++, z+= N-2)
3084 : {
3085 5704332 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
3086 1194828 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
3087 : }
3088 180561 : N = (l-2) % (N-2) + 2;
3089 463295 : for (j=2; j<N; j++) t[j] = z[j];
3090 180561 : gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
3091 180561 : return normalizepol_lg(x, i+1);
3092 : }
3093 :
3094 : /*******************************************************************/
3095 : /* */
3096 : /* Domain detection */
3097 : /* */
3098 : /*******************************************************************/
3099 :
3100 : static GEN
3101 582113 : RgX_mul_FpX(GEN x, GEN y, GEN p)
3102 : {
3103 582113 : pari_sp av = avma;
3104 : GEN r;
3105 582113 : if (lgefint(p) == 3)
3106 : {
3107 533872 : ulong pp = uel(p, 2);
3108 533872 : r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
3109 : RgX_to_Flx(y, pp), pp));
3110 : }
3111 : else
3112 48241 : r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3113 582113 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3114 536347 : return gc_upto(av, FpX_to_mod(r, p));
3115 : }
3116 :
3117 : static GEN
3118 1477 : RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3119 : {
3120 1477 : pari_sp av = avma;
3121 : long dT;
3122 : GEN kx, ky, r;
3123 1477 : GEN T = RgX_to_FpX(pol, p);
3124 1477 : if (signe(T)==0) pari_err_OP("*", x, y);
3125 1470 : dT = degpol(T);
3126 1470 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3127 1470 : ky = RgXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
3128 1470 : r = FpX_mul(kx, ky, p);
3129 1470 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3130 1386 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3131 : }
3132 :
3133 : static GEN
3134 444256 : RgX_liftred(GEN x, GEN T)
3135 444256 : { return RgXQX_red(liftpol_shallow(x), T); }
3136 :
3137 : static GEN
3138 165568 : RgX_mul_QXQX(GEN x, GEN y, GEN T)
3139 : {
3140 165568 : pari_sp av = avma;
3141 165568 : long dT = degpol(T);
3142 165568 : GEN r = QX_mul(RgXX_to_Kronecker(RgX_liftred(x, T), dT),
3143 : RgXX_to_Kronecker(RgX_liftred(y, T), dT));
3144 165568 : return gc_upto(av, Kronecker_to_mod(r, T));
3145 : }
3146 :
3147 : static GEN
3148 16975 : RgX_sqr_FpX(GEN x, GEN p)
3149 : {
3150 16975 : pari_sp av = avma;
3151 : GEN r;
3152 16975 : if (lgefint(p) == 3)
3153 : {
3154 16771 : ulong pp = uel(p, 2);
3155 16771 : r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
3156 : }
3157 : else
3158 204 : r = FpX_sqr(RgX_to_FpX(x, p), p);
3159 16975 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3160 16345 : return gc_upto(av, FpX_to_mod(r, p));
3161 : }
3162 :
3163 : static GEN
3164 196 : RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
3165 : {
3166 196 : pari_sp av = avma;
3167 : long dT;
3168 196 : GEN kx, r, T = RgX_to_FpX(pol, p);
3169 196 : if (signe(T)==0) pari_err_OP("*",x,x);
3170 189 : dT = degpol(T);
3171 189 : kx = RgXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
3172 189 : r = FpX_sqr(kx, p);
3173 189 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3174 189 : return gc_upto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
3175 : }
3176 :
3177 : static GEN
3178 13418 : RgX_sqr_QXQX(GEN x, GEN T)
3179 : {
3180 13418 : pari_sp av = avma;
3181 13418 : long dT = degpol(T);
3182 13418 : GEN r = QX_sqr(RgXX_to_Kronecker(RgX_liftred(x, T), dT));
3183 13418 : return gc_upto(av, Kronecker_to_mod(r, T));
3184 : }
3185 :
3186 : static GEN
3187 67781 : RgX_rem_FpX(GEN x, GEN y, GEN p)
3188 : {
3189 67781 : pari_sp av = avma;
3190 : GEN r;
3191 67781 : if (lgefint(p) == 3)
3192 : {
3193 51502 : ulong pp = uel(p, 2);
3194 51502 : r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
3195 : RgX_to_Flx(y, pp), pp));
3196 : }
3197 : else
3198 16279 : r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3199 67781 : if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
3200 30905 : return gc_upto(av, FpX_to_mod(r, p));
3201 : }
3202 :
3203 : static GEN
3204 49851 : RgX_rem_QXQX(GEN x, GEN y, GEN T)
3205 : {
3206 49851 : pari_sp av = avma;
3207 : GEN r;
3208 49851 : r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
3209 49851 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3210 : }
3211 : static GEN
3212 70 : RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3213 : {
3214 70 : pari_sp av = avma;
3215 : GEN r;
3216 70 : GEN T = RgX_to_FpX(pol, p);
3217 70 : if (signe(T) == 0) pari_err_OP("%", x, y);
3218 63 : if (lgefint(p) == 3)
3219 : {
3220 55 : ulong pp = uel(p, 2);
3221 55 : GEN Tp = ZX_to_Flx(T, pp);
3222 55 : r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
3223 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3224 : }
3225 : else
3226 8 : r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3227 63 : if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
3228 56 : return gc_upto(av, FpXQX_to_mod(r, T, p));
3229 : }
3230 :
3231 : static GEN
3232 84367981 : RgX_mul_fast(GEN x, GEN y)
3233 : {
3234 : GEN p, pol;
3235 : long pa;
3236 84367981 : long t = RgX_type2(x,y, &p,&pol,&pa);
3237 84368207 : switch(t)
3238 : {
3239 34902828 : case t_INT: return ZX_mul(x,y);
3240 3572131 : case t_FRAC: return QX_mul(x,y);
3241 103959 : case t_FFELT: return FFX_mul(x, y, pol);
3242 582113 : case t_INTMOD: return RgX_mul_FpX(x, y, p);
3243 165568 : case RgX_type_code(t_POLMOD, t_INT):
3244 : case RgX_type_code(t_POLMOD, t_FRAC):
3245 165568 : return RgX_mul_QXQX(x, y, pol);
3246 1455 : case RgX_type_code(t_POLMOD, t_INTMOD):
3247 1455 : return RgX_mul_FpXQX(x, y, pol, p);
3248 45040153 : default: return NULL;
3249 : }
3250 : }
3251 : static GEN
3252 1511270 : RgX_sqr_fast(GEN x)
3253 : {
3254 : GEN p, pol;
3255 : long pa;
3256 1511270 : long t = RgX_type(x,&p,&pol,&pa);
3257 1511274 : switch(t)
3258 : {
3259 1355814 : case t_INT: return ZX_sqr(x);
3260 89137 : case t_FRAC: return QX_sqr(x);
3261 2499 : case t_FFELT: return FFX_sqr(x, pol);
3262 16975 : case t_INTMOD: return RgX_sqr_FpX(x, p);
3263 13418 : case RgX_type_code(t_POLMOD, t_INT):
3264 : case RgX_type_code(t_POLMOD, t_FRAC):
3265 13418 : return RgX_sqr_QXQX(x, pol);
3266 195 : case RgX_type_code(t_POLMOD, t_INTMOD):
3267 195 : return RgX_sqr_FpXQX(x, pol, p);
3268 33236 : default: return NULL;
3269 : }
3270 : }
3271 :
3272 : static GEN
3273 14658262 : RgX_rem_fast(GEN x, GEN y)
3274 : {
3275 : GEN p, pol;
3276 : long pa;
3277 14658262 : long t = RgX_type2(x,y, &p,&pol,&pa);
3278 14658342 : switch(t)
3279 : {
3280 3453764 : case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3281 1836586 : case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3282 84 : case t_FFELT: return FFX_rem(x, y, pol);
3283 67781 : case t_INTMOD: return RgX_rem_FpX(x, y, p);
3284 49851 : case RgX_type_code(t_POLMOD, t_INT):
3285 : case RgX_type_code(t_POLMOD, t_FRAC):
3286 49851 : return RgX_rem_QXQX(x, y, pol);
3287 61 : case RgX_type_code(t_POLMOD, t_INTMOD):
3288 61 : return RgX_rem_FpXQX(x, y, pol, p);
3289 9250215 : default: return NULL;
3290 : }
3291 : }
3292 :
3293 : GEN
3294 69243250 : RgX_mul(GEN x, GEN y)
3295 : {
3296 69243250 : GEN z = RgX_mul_fast(x,y);
3297 69243364 : if (!z) z = RgX_mul_i(x,y);
3298 69243033 : return z;
3299 : }
3300 :
3301 : GEN
3302 1499118 : RgX_sqr(GEN x)
3303 : {
3304 1499118 : GEN z = RgX_sqr_fast(x);
3305 1499115 : if (!z) z = RgX_sqr_i(x);
3306 1499115 : return z;
3307 : }
3308 :
3309 : GEN
3310 14658278 : RgX_rem(GEN x, GEN y)
3311 : {
3312 14658278 : GEN z = RgX_rem_fast(x, y);
3313 14658332 : if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3314 14658302 : return z;
3315 : }
3316 :
3317 : static GEN
3318 0 : _RgX_mul(void* E, GEN x, GEN y)
3319 0 : { (void) E; return RgX_mul(x, y); }
3320 :
3321 : GEN
3322 0 : RgXV_prod(GEN V)
3323 0 : { return gen_product(V, NULL, &_RgX_mul); }
3324 :
3325 : GEN
3326 12219932 : RgXn_mul(GEN f, GEN g, long n)
3327 : {
3328 12219932 : pari_sp av = avma;
3329 12219932 : GEN h = RgX_mul_fast(f,g);
3330 12219933 : if (!h) return RgXn_mul2(f,g,n);
3331 1796109 : if (degpol(h) < n) return h;
3332 428238 : return gc_GEN(av, RgXn_red_shallow(h, n));
3333 : }
3334 :
3335 : GEN
3336 12152 : RgXn_sqr(GEN f, long n)
3337 : {
3338 12152 : pari_sp av = avma;
3339 12152 : GEN g = RgX_sqr_fast(f);
3340 12152 : if (!g) return RgXn_sqr2(f,n);
3341 11837 : if (degpol(g) < n) return g;
3342 10640 : return gc_GEN(av, RgXn_red_shallow(g, n));
3343 : }
3344 :
3345 : /* (f*g) \/ x^n */
3346 : GEN
3347 2904794 : RgX_mulhigh_i(GEN f, GEN g, long n)
3348 : {
3349 2904794 : GEN h = RgX_mul_fast(f,g);
3350 2904793 : return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3351 : }
3352 :
3353 : /* (f*g) \/ x^n */
3354 : GEN
3355 0 : RgX_sqrhigh_i(GEN f, long n)
3356 : {
3357 0 : GEN h = RgX_sqr_fast(f);
3358 0 : return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3359 : }
|
