"Fossies" - the Fresh Open Source Software Archive 
Member "pari-2.13.1/src/basemath/RgX.c" (14 Jan 2021, 77499 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "RgX.c" see the
Fossies "Dox" file reference documentation and the latest
Fossies "Diffs" side-by-side code changes report:
2.13.0_vs_2.13.1.
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 /*******************************************************************/
19 /* */
20 /* GENERIC */
21 /* */
22 /*******************************************************************/
23
24 /* Return optimal parameter l for the evaluation of n/m polynomials of degree d
25 Fractional values can be used if the evaluations are done with different
26 accuracies, and thus have different weights.
27 */
28 long
29 brent_kung_optpow(long d, long n, long m)
30 {
31 long p, r;
32 long pold=1, rold=n*(d-1);
33 for(p=2; p<=d; p++)
34 {
35 r = m*(p-1) + n*((d-1)/p);
36 if (r<rold) { pold=p; rold=r; }
37 }
38 return pold;
39 }
40
41 static GEN
42 gen_RgXQ_eval_powers(GEN P, GEN V, long a, long n, void *E, const struct bb_algebra *ff,
43 GEN cmul(void *E, GEN P, long a, GEN x))
44 {
45 pari_sp av = avma;
46 long i;
47 GEN z = cmul(E,P,a,ff->one(E));
48 if (!z) z = gen_0;
49 for (i=1; i<=n; i++)
50 {
51 GEN t = cmul(E,P,a+i,gel(V,i+1));
52 if (t) {
53 z = ff->add(E, z, t);
54 if (gc_needed(av,2)) z = gerepileupto(av, z);
55 }
56 }
57 return ff->red(E,z);
58 }
59
60 /* Brent & Kung
61 * (Fast algorithms for manipulating formal power series, JACM 25:581-595, 1978)
62 *
63 * V as output by FpXQ_powers(x,l,T,p). For optimal performance, l is as given
64 * by brent_kung_optpow */
65 GEN
66 gen_bkeval_powers(GEN P, long d, GEN V, void *E, const struct bb_algebra *ff,
67 GEN cmul(void *E, GEN P, long a, GEN x))
68 {
69 pari_sp av = avma;
70 long l = lg(V)-1;
71 GEN z, u;
72
73 if (d < 0) return ff->zero(E);
74 if (d < l) return gerepileupto(av, gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul));
75 if (l<2) pari_err_DOMAIN("gen_RgX_bkeval_powers", "#powers", "<",gen_2,V);
76 if (DEBUGLEVEL>=8)
77 {
78 long cnt = 1 + (d - l) / (l-1);
79 err_printf("RgX_RgXQV_eval(%ld/%ld): %ld RgXQ_mul\n", d, l-1, cnt);
80 }
81 d -= l;
82 z = gen_RgXQ_eval_powers(P,V,d+1,l-1,E,ff,cmul);
83 while (d >= l-1)
84 {
85 d -= l-1;
86 u = gen_RgXQ_eval_powers(P,V,d+1,l-2,E,ff,cmul);
87 z = ff->add(E,u, ff->mul(E,z,gel(V,l)));
88 if (gc_needed(av,2))
89 z = gerepileupto(av, z);
90 }
91 u = gen_RgXQ_eval_powers(P,V,0,d,E,ff,cmul);
92 z = ff->add(E,u, ff->mul(E,z,gel(V,d+2)));
93 return gerepileupto(av, ff->red(E,z));
94 }
95
96 GEN
97 gen_bkeval(GEN Q, long d, GEN x, int use_sqr, void *E, const struct bb_algebra *ff,
98 GEN cmul(void *E, GEN P, long a, GEN x))
99 {
100 pari_sp av = avma;
101 GEN z, V;
102 long rtd;
103 if (d < 0) return ff->zero(E);
104 rtd = (long) sqrt((double)d);
105 V = gen_powers(x,rtd,use_sqr,E,ff->sqr,ff->mul,ff->one);
106 z = gen_bkeval_powers(Q, d, V, E, ff, cmul);
107 return gerepileupto(av, z);
108 }
109
110 static GEN
111 _gen_nored(void *E, GEN x) { (void)E; return x; }
112 static GEN
113 _gen_add(void *E, GEN x, GEN y) { (void)E; return gadd(x, y); }
114 static GEN
115 _gen_sub(void *E, GEN x, GEN y) { (void)E; return gsub(x, y); }
116 static GEN
117 _gen_mul(void *E, GEN x, GEN y) { (void)E; return gmul(x, y); }
118 static GEN
119 _gen_sqr(void *E, GEN x) { (void)E; return gsqr(x); }
120 static GEN
121 _gen_one(void *E) { (void)E; return gen_1; }
122 static GEN
123 _gen_zero(void *E) { (void)E; return gen_0; }
124
125 static struct bb_algebra Rg_algebra = { _gen_nored, _gen_add, _gen_sub,
126 _gen_mul, _gen_sqr,_gen_one,_gen_zero };
127
128 static GEN
129 _gen_cmul(void *E, GEN P, long a, GEN x)
130 {(void)E; return gmul(gel(P,a+2), x);}
131
132 GEN
133 RgX_RgV_eval(GEN Q, GEN x)
134 {
135 return gen_bkeval_powers(Q, degpol(Q), x, NULL, &Rg_algebra, _gen_cmul);
136 }
137
138 GEN
139 RgX_Rg_eval_bk(GEN Q, GEN x)
140 {
141 return gen_bkeval(Q, degpol(Q), x, 1, NULL, &Rg_algebra, _gen_cmul);
142 }
143
144 GEN
145 RgXV_RgV_eval(GEN Q, GEN x)
146 {
147 long i, l = lg(Q), vQ = gvar(Q);
148 GEN v = cgetg(l, t_VEC);
149 for (i = 1; i < l; i++)
150 {
151 GEN Qi = gel(Q, i);
152 gel(v, i) = typ(Qi)==t_POL && varn(Qi)==vQ? RgX_RgV_eval(Qi, x): gcopy(Qi);
153 }
154 return v;
155 }
156
157 GEN
158 RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
159 {
160 pari_sp av = avma;
161 long d, i, v;
162 GEN s;
163 if (typ(P)!=t_POL)
164 return mkvec2(P, gen_1);
165 d = degpol(P); v = varn(A);
166 s = scalarpol_shallow(gel(P, d+2), v);
167 for (i = d-1; i >= 0; i--)
168 {
169 s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
170 if (gc_needed(av,1))
171 {
172 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
173 s = gerepileupto(av, s);
174 }
175 }
176 s = gerepileupto(av, s);
177 return mkvec2(s, gel(B,d+1));
178 }
179
180 GEN
181 QXQX_homogenous_evalpow(GEN P, GEN A, GEN B, GEN T)
182 {
183 pari_sp av = avma;
184 long i, d = degpol(P), v = varn(A);
185 GEN s;
186 if (signe(P)==0) return mkvec2(pol_0(v), pol_1(v));
187 s = scalarpol_shallow(gel(P, d+2), v);
188 for (i = d-1; i >= 0; i--)
189 {
190 GEN c = gel(P,i+2), b = gel(B,d+1-i);
191 s = RgX_add(QXQX_mul(s, A, T), typ(c)==t_POL ? QXQX_QXQ_mul(b, c, T): gmul(b, c));
192 if (gc_needed(av,1))
193 {
194 if(DEBUGMEM>1) pari_warn(warnmem,"QXQX_homogenous_eval(%ld)",i);
195 s = gerepileupto(av, s);
196 }
197 }
198 s = gerepileupto(av, s);
199 return mkvec2(s, gel(B,d+1));
200 }
201
202 const struct bb_algebra *
203 get_Rg_algebra(void)
204 {
205 return &Rg_algebra;
206 }
207
208 static struct bb_ring Rg_ring = { _gen_add, _gen_mul, _gen_sqr };
209
210 static GEN
211 _RgX_divrem(void *E, GEN x, GEN y, GEN *r)
212 {
213 (void) E;
214 return RgX_divrem(x, y, r);
215 }
216
217 GEN
218 RgX_digits(GEN x, GEN T)
219 {
220 pari_sp av = avma;
221 long d = degpol(T), n = (lgpol(x)+d-1)/d;
222 GEN z = gen_digits(x,T,n,NULL, &Rg_ring, _RgX_divrem);
223 return gerepileupto(av, z);
224 }
225
226 /*******************************************************************/
227 /* */
228 /* RgX */
229 /* */
230 /*******************************************************************/
231
232 long
233 RgX_equal(GEN x, GEN y)
234 {
235 long i = lg(x);
236
237 if (i != lg(y)) return 0;
238 for (i--; i > 1; i--)
239 if (!gequal(gel(x,i),gel(y,i))) return 0;
240 return 1;
241 }
242
243 /* Returns 1 in the base ring over which x is defined */
244 /* HACK: this also works for t_SER */
245 GEN
246 Rg_get_1(GEN x)
247 {
248 GEN p, T;
249 long i, lx, tx = Rg_type(x, &p, &T, &lx);
250 if (RgX_type_is_composite(tx))
251 RgX_type_decode(tx, &i /*junk*/, &tx);
252 switch(tx)
253 {
254 case t_INTMOD: retmkintmod(is_pm1(p)? gen_0: gen_1, icopy(p));
255 case t_PADIC: return cvtop(gen_1, p, lx);
256 case t_FFELT: return FF_1(T);
257 default: return gen_1;
258 }
259 }
260 /* Returns 0 in the base ring over which x is defined */
261 /* HACK: this also works for t_SER */
262 GEN
263 Rg_get_0(GEN x)
264 {
265 GEN p, T;
266 long i, lx, tx = Rg_type(x, &p, &T, &lx);
267 if (RgX_type_is_composite(tx))
268 RgX_type_decode(tx, &i /*junk*/, &tx);
269 switch(tx)
270 {
271 case t_INTMOD: retmkintmod(gen_0, icopy(p));
272 case t_PADIC: return zeropadic(p, lx);
273 case t_FFELT: return FF_zero(T);
274 default: return gen_0;
275 }
276 }
277
278 GEN
279 QX_ZXQV_eval(GEN P, GEN V, GEN dV)
280 {
281 long i, n = degpol(P);
282 GEN z, dz, dP;
283 if (n < 0) return gen_0;
284 P = Q_remove_denom(P, &dP);
285 z = gel(P,2); if (n == 0) return icopy(z);
286 if (dV) z = mulii(dV, z); /* V[1] = dV */
287 z = ZX_Z_add_shallow(ZX_Z_mul(gel(V,2),gel(P,3)), z);
288 for (i=2; i<=n; i++) z = ZX_add(ZX_Z_mul(gel(V,i+1),gel(P,2+i)), z);
289 dz = mul_denom(dP, dV);
290 return dz? RgX_Rg_div(z, dz): z;
291 }
292
293 /* Return P(h * x), not memory clean */
294 GEN
295 RgX_unscale(GEN P, GEN h)
296 {
297 long i, l = lg(P);
298 GEN hi = gen_1, Q = cgetg(l, t_POL);
299 Q[1] = P[1];
300 if (l == 2) return Q;
301 gel(Q,2) = gcopy(gel(P,2));
302 for (i=3; i<l; i++)
303 {
304 hi = gmul(hi,h);
305 gel(Q,i) = gmul(gel(P,i), hi);
306 }
307 return Q;
308 }
309 /* P a ZX, Return P(h * x), not memory clean; optimize for h = -1 */
310 GEN
311 ZX_z_unscale(GEN P, long h)
312 {
313 long i, l = lg(P);
314 GEN Q = cgetg(l, t_POL);
315 Q[1] = P[1];
316 if (l == 2) return Q;
317 gel(Q,2) = gel(P,2);
318 if (l == 3) return Q;
319 if (h == -1)
320 for (i = 3; i < l; i++)
321 {
322 gel(Q,i) = negi(gel(P,i));
323 if (++i == l) break;
324 gel(Q,i) = gel(P,i);
325 }
326 else
327 {
328 GEN hi;
329 gel(Q,3) = mulis(gel(P,3), h);
330 hi = sqrs(h);
331 for (i = 4; i < l; i++)
332 {
333 gel(Q,i) = mulii(gel(P,i), hi);
334 if (i != l-1) hi = mulis(hi,h);
335 }
336 }
337 return Q;
338 }
339 /* P a ZX, h a t_INT. Return P(h * x), not memory clean; optimize for h = -1 */
340 GEN
341 ZX_unscale(GEN P, GEN h)
342 {
343 long i, l;
344 GEN Q, hi;
345 i = itos_or_0(h); if (i) return ZX_z_unscale(P, i);
346 l = lg(P); Q = cgetg(l, t_POL);
347 Q[1] = P[1];
348 if (l == 2) return Q;
349 gel(Q,2) = gel(P,2);
350 if (l == 3) return Q;
351 hi = h;
352 gel(Q,3) = mulii(gel(P,3), hi);
353 for (i = 4; i < l; i++)
354 {
355 hi = mulii(hi,h);
356 gel(Q,i) = mulii(gel(P,i), hi);
357 }
358 return Q;
359 }
360 /* P a ZX. Return P(x << n), not memory clean */
361 GEN
362 ZX_unscale2n(GEN P, long n)
363 {
364 long i, ni = n, l = lg(P);
365 GEN Q = cgetg(l, t_POL);
366 Q[1] = P[1];
367 if (l == 2) return Q;
368 gel(Q,2) = gel(P,2);
369 if (l == 3) return Q;
370 gel(Q,3) = shifti(gel(P,3), ni);
371 for (i=4; i<l; i++)
372 {
373 ni += n;
374 gel(Q,i) = shifti(gel(P,i), ni);
375 }
376 return Q;
377 }
378 /* P(h*X) / h, assuming h | P(0), i.e. the result is a ZX */
379 GEN
380 ZX_unscale_div(GEN P, GEN h)
381 {
382 long i, l = lg(P);
383 GEN hi, Q = cgetg(l, t_POL);
384 Q[1] = P[1];
385 if (l == 2) return Q;
386 gel(Q,2) = diviiexact(gel(P,2), h);
387 if (l == 3) return Q;
388 gel(Q,3) = gel(P,3);
389 if (l == 4) return Q;
390 hi = h;
391 gel(Q,4) = mulii(gel(P,4), hi);
392 for (i=5; i<l; i++)
393 {
394 hi = mulii(hi,h);
395 gel(Q,i) = mulii(gel(P,i), hi);
396 }
397 return Q;
398 }
399
400 GEN
401 RgXV_unscale(GEN v, GEN h)
402 {
403 long i, l;
404 GEN w;
405 if (!h || isint1(h)) return v;
406 w = cgetg_copy(v, &l);
407 for (i=1; i<l; i++) gel(w,i) = RgX_unscale(gel(v,i), h);
408 return w;
409 }
410
411 /* Return h^degpol(P) P(x / h), not memory clean */
412 GEN
413 RgX_rescale(GEN P, GEN h)
414 {
415 long i, l = lg(P);
416 GEN Q = cgetg(l,t_POL), hi = h;
417 gel(Q,l-1) = gel(P,l-1);
418 for (i=l-2; i>=2; i--)
419 {
420 gel(Q,i) = gmul(gel(P,i), hi);
421 if (i == 2) break;
422 hi = gmul(hi,h);
423 }
424 Q[1] = P[1]; return Q;
425 }
426
427 /* A(X^d) --> A(X) */
428 GEN
429 RgX_deflate(GEN x0, long d)
430 {
431 GEN z, y, x;
432 long i,id, dy, dx = degpol(x0);
433 if (d == 1 || dx <= 0) return leafcopy(x0);
434 dy = dx/d;
435 y = cgetg(dy+3, t_POL); y[1] = x0[1];
436 z = y + 2;
437 x = x0+ 2;
438 for (i=id=0; i<=dy; i++,id+=d) gel(z,i) = gel(x,id);
439 return y;
440 }
441
442 /* F a t_RFRAC */
443 long
444 rfrac_deflate_order(GEN F)
445 {
446 GEN N = gel(F,1), D = gel(F,2);
447 long m = (degpol(D) <= 0)? 0: RgX_deflate_order(D);
448 if (m == 1) return 1;
449 if (typ(N) == t_POL && varn(N) == varn(D))
450 m = cgcd(m, RgX_deflate_order(N));
451 return m;
452 }
453 /* F a t_RFRAC */
454 GEN
455 rfrac_deflate_max(GEN F, long *m)
456 {
457 *m = rfrac_deflate_order(F);
458 return rfrac_deflate(F, *m);
459 }
460 /* F a t_RFRAC */
461 GEN
462 rfrac_deflate(GEN F, long m)
463 {
464 GEN N = gel(F,1), D = gel(F,2);
465 if (m == 1) return F;
466 if (typ(N) == t_POL && varn(N) == varn(D)) N = RgX_deflate(N, m);
467 D = RgX_deflate(D, m); return mkrfrac(N, D);
468 }
469
470 /* return x0(X^d) */
471 GEN
472 RgX_inflate(GEN x0, long d)
473 {
474 long i, id, dy, dx = degpol(x0);
475 GEN x = x0 + 2, z, y;
476 if (dx <= 0) return leafcopy(x0);
477 dy = dx*d;
478 y = cgetg(dy+3, t_POL); y[1] = x0[1];
479 z = y + 2;
480 for (i=0; i<=dy; i++) gel(z,i) = gen_0;
481 for (i=id=0; i<=dx; i++,id+=d) gel(z,id) = gel(x,i);
482 return y;
483 }
484
485 /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
486 static GEN
487 RgX_translate_basecase(GEN P, GEN c)
488 {
489 pari_sp av = avma;
490 GEN Q, R;
491 long i, k, n;
492
493 if (!signe(P) || gequal0(c)) return RgX_copy(P);
494 Q = leafcopy(P);
495 R = Q+2; n = degpol(P);
496 if (isint1(c))
497 {
498 for (i=1; i<=n; i++)
499 {
500 for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gel(R,k+1));
501 if (gc_needed(av,2))
502 {
503 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(1), i = %ld/%ld", i,n);
504 Q = gerepilecopy(av, Q); R = Q+2;
505 }
506 }
507 }
508 else if (isintm1(c))
509 {
510 for (i=1; i<=n; i++)
511 {
512 for (k=n-i; k<n; k++) gel(R,k) = gsub(gel(R,k), gel(R,k+1));
513 if (gc_needed(av,2))
514 {
515 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate(-1), i = %ld/%ld", i,n);
516 Q = gerepilecopy(av, Q); R = Q+2;
517 }
518 }
519 }
520 else
521 {
522 for (i=1; i<=n; i++)
523 {
524 for (k=n-i; k<n; k++) gel(R,k) = gadd(gel(R,k), gmul(c, gel(R,k+1)));
525 if (gc_needed(av,2))
526 {
527 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_translate, i = %ld/%ld", i,n);
528 Q = gerepilecopy(av, Q); R = Q+2;
529 }
530 }
531 }
532 return gerepilecopy(av, Q);
533 }
534 GEN
535 RgX_translate(GEN P, GEN c)
536 {
537 pari_sp av = avma;
538 long n = degpol(P);
539 if (n < 40)
540 return RgX_translate_basecase(P, c);
541 else
542 {
543 long d = n >> 1;
544 GEN Q = RgX_translate(RgX_shift_shallow(P, -d), c);
545 GEN R = RgX_translate(RgXn_red_shallow(P, d), c);
546 GEN S = gpowgs(deg1pol_shallow(gen_1, c, varn(P)), d);
547 return gerepileupto(av, RgX_add(RgX_mul(Q, S), R));
548 }
549 }
550
551 /* return lift( P(X + c) ) using Horner, c in R[y]/(T) */
552 GEN
553 RgXQX_translate(GEN P, GEN c, GEN T)
554 {
555 pari_sp av = avma;
556 GEN Q, R;
557 long i, k, n;
558
559 if (!signe(P) || gequal0(c)) return RgX_copy(P);
560 Q = leafcopy(P);
561 R = Q+2; n = degpol(P);
562 for (i=1; i<=n; i++)
563 {
564 for (k=n-i; k<n; k++)
565 {
566 pari_sp av2 = avma;
567 gel(R,k) = gerepileupto(av2,
568 RgX_rem(gadd(gel(R,k), gmul(c, gel(R,k+1))), T));
569 }
570 if (gc_needed(av,2))
571 {
572 if(DEBUGMEM>1) pari_warn(warnmem,"RgXQX_translate, i = %ld/%ld", i,n);
573 Q = gerepilecopy(av, Q); R = Q+2;
574 }
575 }
576 return gerepilecopy(av, Q);
577 }
578
579 /********************************************************************/
580 /** **/
581 /** CONVERSIONS **/
582 /** (not memory clean) **/
583 /** **/
584 /********************************************************************/
585 /* to INT / FRAC / (POLMOD mod T), not memory clean because T not copied,
586 * but everything else is */
587 static GEN
588 QXQ_to_mod(GEN x, GEN T)
589 {
590 long d;
591 switch(typ(x))
592 {
593 case t_INT: return icopy(x);
594 case t_FRAC: return gcopy(x);
595 case t_POL:
596 d = degpol(x);
597 if (d < 0) return gen_0;
598 if (d == 0) return gcopy(gel(x,2));
599 return mkpolmod(RgX_copy(x), T);
600 default: pari_err_TYPE("QXQ_to_mod",x);
601 return NULL;/* LCOV_EXCL_LINE */
602 }
603 }
604 /* pure shallow version */
605 GEN
606 QXQ_to_mod_shallow(GEN x, GEN T)
607 {
608 long d;
609 switch(typ(x))
610 {
611 case t_INT:
612 case t_FRAC: return x;
613 case t_POL:
614 d = degpol(x);
615 if (d < 0) return gen_0;
616 if (d == 0) return gel(x,2);
617 return mkpolmod(x, T);
618 default: pari_err_TYPE("QXQ_to_mod",x);
619 return NULL;/* LCOV_EXCL_LINE */
620 }
621 }
622 /* T a ZX, z lifted from (Q[Y]/(T(Y)))[X], apply QXQ_to_mod to all coeffs.
623 * Not memory clean because T not copied, but everything else is */
624 static GEN
625 QXQX_to_mod(GEN z, GEN T)
626 {
627 long i,l = lg(z);
628 GEN x = cgetg(l,t_POL);
629 for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod(gel(z,i), T);
630 x[1] = z[1]; return normalizepol_lg(x,l);
631 }
632 /* pure shallow version */
633 GEN
634 QXQX_to_mod_shallow(GEN z, GEN T)
635 {
636 long i,l = lg(z);
637 GEN x = cgetg(l,t_POL);
638 for (i=2; i<l; i++) gel(x,i) = QXQ_to_mod_shallow(gel(z,i), T);
639 x[1] = z[1]; return normalizepol_lg(x,l);
640 }
641 /* Apply QXQX_to_mod to all entries. Memory-clean ! */
642 GEN
643 QXQXV_to_mod(GEN V, GEN T)
644 {
645 long i, l = lg(V);
646 GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
647 for (i=1;i<l; i++) gel(z,i) = QXQX_to_mod(gel(V,i), T);
648 return z;
649 }
650 /* Apply QXQ_to_mod to all entries. Memory-clean ! */
651 GEN
652 QXQV_to_mod(GEN V, GEN T)
653 {
654 long i, l = lg(V);
655 GEN z = cgetg(l, t_VEC); T = ZX_copy(T);
656 for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod(gel(V,i), T);
657 return z;
658 }
659
660 /* Apply QXQ_to_mod to all entries. Memory-clean ! */
661 GEN
662 QXQC_to_mod_shallow(GEN V, GEN T)
663 {
664 long i, l = lg(V);
665 GEN z = cgetg(l, t_COL);
666 for (i=1;i<l; i++) gel(z,i) = QXQ_to_mod_shallow(gel(V,i), T);
667 return z;
668 }
669
670 GEN
671 QXQM_to_mod_shallow(GEN V, GEN T)
672 {
673 long i, l = lg(V);
674 GEN z = cgetg(l, t_MAT);
675 for (i=1; i<l; i++) gel(z,i) = QXQC_to_mod_shallow(gel(V,i), T);
676 return z;
677 }
678
679 GEN
680 RgX_renormalize_lg(GEN x, long lx)
681 {
682 long i;
683 for (i = lx-1; i>1; i--)
684 if (! gequal0(gel(x,i))) break; /* _not_ isexactzero */
685 stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
686 setlg(x, i+1); setsigne(x, i != 1); return x;
687 }
688
689 GEN
690 RgV_to_RgX(GEN x, long v)
691 {
692 long i, k = lg(x);
693 GEN p;
694
695 while (--k && gequal0(gel(x,k)));
696 if (!k) return pol_0(v);
697 i = k+2; p = cgetg(i,t_POL);
698 p[1] = evalsigne(1) | evalvarn(v);
699 x--; for (k=2; k<i; k++) gel(p,k) = gel(x,k);
700 return p;
701 }
702 GEN
703 RgV_to_RgX_reverse(GEN x, long v)
704 {
705 long j, k, l = lg(x);
706 GEN p;
707
708 for (k = 1; k < l; k++)
709 if (!gequal0(gel(x,k))) break;
710 if (k == l) return pol_0(v);
711 k -= 1;
712 l -= k;
713 x += k;
714 p = cgetg(l+1,t_POL);
715 p[1] = evalsigne(1) | evalvarn(v);
716 for (j=2, k=l; j<=l; j++) gel(p,j) = gel(x,--k);
717 return p;
718 }
719
720 /* return the (N-dimensional) vector of coeffs of p */
721 GEN
722 RgX_to_RgC(GEN x, long N)
723 {
724 long i, l;
725 GEN z;
726 l = lg(x)-1; x++;
727 if (l > N+1) l = N+1; /* truncate higher degree terms */
728 z = cgetg(N+1,t_COL);
729 for (i=1; i<l ; i++) gel(z,i) = gel(x,i);
730 for ( ; i<=N; i++) gel(z,i) = gen_0;
731 return z;
732 }
733 GEN
734 Rg_to_RgC(GEN x, long N)
735 {
736 return (typ(x) == t_POL)? RgX_to_RgC(x,N): scalarcol_shallow(x, N);
737 }
738
739 /* vector of polynomials (in v) whose coeffs are given by the columns of x */
740 GEN
741 RgM_to_RgXV(GEN x, long v)
742 { pari_APPLY_type(t_VEC, RgV_to_RgX(gel(x,i), v)) }
743
744 /* matrix whose entries are given by the coeffs of the polynomials in
745 * vector v (considered as degree n-1 polynomials) */
746 GEN
747 RgV_to_RgM(GEN x, long n)
748 { pari_APPLY_type(t_MAT, Rg_to_RgC(gel(x,i), n)) }
749
750 GEN
751 RgXV_to_RgM(GEN x, long n)
752 { pari_APPLY_type(t_MAT, RgX_to_RgC(gel(x,i), n)) }
753
754 /* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */
755 GEN
756 RgM_to_RgXX(GEN x, long v,long w)
757 {
758 long j, lx = lg(x);
759 GEN y = cgetg(lx+1, t_POL);
760 y[1] = evalsigne(1) | evalvarn(v);
761 y++;
762 for (j=1; j<lx; j++) gel(y,j) = RgV_to_RgX(gel(x,j), w);
763 return normalizepol_lg(--y, lx+1);
764 }
765
766 /* matrix whose entries are given by the coeffs of the polynomial v in
767 * two variables (considered as degree n-1 polynomials) */
768 GEN
769 RgXX_to_RgM(GEN v, long n)
770 {
771 long j, N = lg(v)-1;
772 GEN y = cgetg(N, t_MAT);
773 for (j=1; j<N; j++) gel(y,j) = Rg_to_RgC(gel(v,j+1), n);
774 return y;
775 }
776
777 /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
778 GEN
779 RgXY_swapspec(GEN x, long n, long w, long nx)
780 {
781 long j, ly = n+3;
782 GEN y = cgetg(ly, t_POL);
783 y[1] = evalsigne(1);
784 for (j=2; j<ly; j++)
785 {
786 long k;
787 GEN a = cgetg(nx+2,t_POL);
788 a[1] = evalsigne(1) | evalvarn(w);
789 for (k=0; k<nx; k++)
790 {
791 GEN xk = gel(x,k);
792 if (typ(xk)==t_POL)
793 gel(a,k+2) = j<lg(xk)? gel(xk,j): gen_0;
794 else
795 gel(a,k+2) = j==2 ? xk: gen_0;
796 }
797 gel(y,j) = normalizepol_lg(a, nx+2);
798 }
799 return normalizepol_lg(y,ly);
800 }
801
802 /* P(X,Y) --> P(Y,X), n is an upper bound for deg_Y(P) */
803 GEN
804 RgXY_swap(GEN x, long n, long w)
805 {
806 GEN z = RgXY_swapspec(x+2, n, w, lgpol(x));
807 setvarn(z, varn(x)); return z;
808 }
809
810 long
811 RgXY_degreex(GEN b)
812 {
813 long deg = 0, i;
814 if (!signe(b)) return -1;
815 for (i = 2; i < lg(b); ++i)
816 {
817 GEN bi = gel(b, i);
818 if (typ(bi) == t_POL)
819 deg = maxss(deg, degpol(bi));
820 }
821 return deg;
822 }
823
824 /* return (x % X^n). Shallow */
825 GEN
826 RgXn_red_shallow(GEN a, long n)
827 {
828 long i, L = n+2, l = lg(a);
829 GEN b;
830 if (L >= l) return a; /* deg(x) < n */
831 b = cgetg(L, t_POL); b[1] = a[1];
832 for (i=2; i<L; i++) gel(b,i) = gel(a,i);
833 return normalizepol_lg(b,L);
834 }
835
836 GEN
837 RgXnV_red_shallow(GEN x, long n)
838 { pari_APPLY_type(t_VEC, RgXn_red_shallow(gel(x,i), n)) }
839
840 /* return (x * X^n). Shallow */
841 GEN
842 RgX_shift_shallow(GEN a, long n)
843 {
844 long i, l = lg(a);
845 GEN b;
846 if (l == 2 || !n) return a;
847 l += n;
848 if (n < 0)
849 {
850 if (l <= 2) return pol_0(varn(a));
851 b = cgetg(l, t_POL); b[1] = a[1];
852 a -= n;
853 for (i=2; i<l; i++) gel(b,i) = gel(a,i);
854 } else {
855 b = cgetg(l, t_POL); b[1] = a[1];
856 a -= n; n += 2;
857 for (i=2; i<n; i++) gel(b,i) = gen_0;
858 for ( ; i<l; i++) gel(b,i) = gel(a,i);
859 }
860 return b;
861 }
862 /* return (x * X^n). */
863 GEN
864 RgX_shift(GEN a, long n)
865 {
866 long i, l = lg(a);
867 GEN b;
868 if (l == 2 || !n) return RgX_copy(a);
869 l += n;
870 if (n < 0)
871 {
872 if (l <= 2) return pol_0(varn(a));
873 b = cgetg(l, t_POL); b[1] = a[1];
874 a -= n;
875 for (i=2; i<l; i++) gel(b,i) = gcopy(gel(a,i));
876 } else {
877 b = cgetg(l, t_POL); b[1] = a[1];
878 a -= n; n += 2;
879 for (i=2; i<n; i++) gel(b,i) = gen_0;
880 for ( ; i<l; i++) gel(b,i) = gcopy(gel(a,i));
881 }
882 return b;
883 }
884
885 GEN
886 RgX_rotate_shallow(GEN P, long k, long p)
887 {
888 long i, l = lgpol(P);
889 GEN r;
890 if (signe(P)==0)
891 return pol_0(varn(P));
892 r = cgetg(p+2,t_POL); r[1] = P[1];
893 for(i=0; i<p; i++)
894 {
895 long s = 2+(i+k)%p;
896 gel(r,s) = i<l? gel(P,2+i): gen_0;
897 }
898 return RgX_renormalize(r);
899 }
900
901 GEN
902 RgX_mulXn(GEN x, long d)
903 {
904 pari_sp av;
905 GEN z;
906 long v;
907 if (d >= 0) return RgX_shift(x, d);
908 d = -d;
909 v = RgX_val(x);
910 if (v >= d) return RgX_shift(x, -d);
911 av = avma;
912 z = gred_rfrac_simple(RgX_shift_shallow(x, -v), pol_xn(d - v, varn(x)));
913 return gerepileupto(av, z);
914 }
915
916 long
917 RgXV_maxdegree(GEN x)
918 {
919 long d = -1, i, l = lg(x);
920 for (i = 1; i < l; i++)
921 d = maxss(d, degpol(gel(x,i)));
922 return d;
923 }
924
925 long
926 RgX_val(GEN x)
927 {
928 long i, lx = lg(x);
929 if (lx == 2) return LONG_MAX;
930 for (i = 2; i < lx; i++)
931 if (!isexactzero(gel(x,i))) break;
932 if (i == lx) return LONG_MAX;/* possible with nonrational zeros */
933 return i - 2;
934 }
935 long
936 RgX_valrem(GEN x, GEN *Z)
937 {
938 long v, i, lx = lg(x);
939 if (lx == 2) { *Z = pol_0(varn(x)); return LONG_MAX; }
940 for (i = 2; i < lx; i++)
941 if (!isexactzero(gel(x,i))) break;
942 /* possible with nonrational zeros */
943 if (i == lx) { *Z = pol_0(varn(x)); return LONG_MAX; }
944 v = i - 2;
945 *Z = RgX_shift_shallow(x, -v);
946 return v;
947 }
948 long
949 RgX_valrem_inexact(GEN x, GEN *Z)
950 {
951 long v;
952 if (!signe(x)) { if (Z) *Z = pol_0(varn(x)); return LONG_MAX; }
953 for (v = 0;; v++)
954 if (!gequal0(gel(x,2+v))) break;
955 if (Z) *Z = RgX_shift_shallow(x, -v);
956 return v;
957 }
958
959 GEN
960 RgXQC_red(GEN x, GEN T)
961 { pari_APPLY_type(t_COL, grem(gel(x,i), T)) }
962
963 GEN
964 RgXQV_red(GEN x, GEN T)
965 { pari_APPLY_type(t_VEC, grem(gel(x,i), T)) }
966
967 GEN
968 RgXQM_red(GEN x, GEN T)
969 { pari_APPLY_same(RgXQC_red(gel(x,i), T)) }
970
971 GEN
972 RgXQM_mul(GEN P, GEN Q, GEN T)
973 {
974 return RgXQM_red(RgM_mul(P, Q), T);
975 }
976
977 GEN
978 RgXQX_red(GEN P, GEN T)
979 {
980 long i, l = lg(P);
981 GEN Q = cgetg(l, t_POL);
982 Q[1] = P[1];
983 for (i=2; i<l; i++) gel(Q,i) = grem(gel(P,i), T);
984 return normalizepol_lg(Q, l);
985 }
986
987 GEN
988 RgX_deriv(GEN x)
989 {
990 long i,lx = lg(x)-1;
991 GEN y;
992
993 if (lx<3) return pol_0(varn(x));
994 y = cgetg(lx,t_POL); gel(y,2) = gcopy(gel(x,3));
995 for (i=3; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
996 y[1] = x[1]; return normalizepol_lg(y,i);
997 }
998
999 GEN
1000 RgX_recipspec_shallow(GEN x, long l, long n)
1001 {
1002 long i;
1003 GEN z = cgetg(n+2,t_POL);
1004 z[1] = 0; z += 2;
1005 for(i=0; i<l; i++)
1006 gel(z,n-i-1) = gel(x,i);
1007 for( ; i<n; i++)
1008 gel(z, n-i-1) = gen_0;
1009 return normalizepol_lg(z-2,n+2);
1010 }
1011
1012 GEN
1013 RgXn_recip_shallow(GEN P, long n)
1014 {
1015 GEN Q = RgX_recipspec_shallow(P+2, lgpol(P), n);
1016 setvarn(Q, varn(P));
1017 return Q;
1018 }
1019
1020 /* return coefficients s.t x = x_0 X^n + ... + x_n */
1021 GEN
1022 RgX_recip(GEN x)
1023 {
1024 long lx, i, j;
1025 GEN y = cgetg_copy(x, &lx);
1026 y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gcopy(gel(x,j));
1027 return normalizepol_lg(y,lx);
1028 }
1029 /* shallow version */
1030 GEN
1031 RgX_recip_shallow(GEN x)
1032 {
1033 long lx, i, j;
1034 GEN y = cgetg_copy(x, &lx);
1035 y[1] = x[1]; for (i=2,j=lx-1; i<lx; i++,j--) gel(y,i) = gel(x,j);
1036 return y;
1037 }
1038 /*******************************************************************/
1039 /* */
1040 /* ADDITION / SUBTRACTION */
1041 /* */
1042 /*******************************************************************/
1043 /* same variable */
1044 GEN
1045 RgX_add(GEN x, GEN y)
1046 {
1047 long i, lx = lg(x), ly = lg(y);
1048 GEN z;
1049 if (ly <= lx) {
1050 z = cgetg(lx,t_POL); z[1] = x[1];
1051 for (i=2; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1052 for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1053 z = normalizepol_lg(z, lx);
1054 } else {
1055 z = cgetg(ly,t_POL); z[1] = y[1];
1056 for (i=2; i < lx; i++) gel(z,i) = gadd(gel(x,i),gel(y,i));
1057 for ( ; i < ly; i++) gel(z,i) = gcopy(gel(y,i));
1058 z = normalizepol_lg(z, ly);
1059 }
1060 return z;
1061 }
1062 GEN
1063 RgX_sub(GEN x, GEN y)
1064 {
1065 long i, lx = lg(x), ly = lg(y);
1066 GEN z;
1067 if (ly <= lx) {
1068 z = cgetg(lx,t_POL); z[1] = x[1];
1069 for (i=2; i < ly; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1070 for ( ; i < lx; i++) gel(z,i) = gcopy(gel(x,i));
1071 z = normalizepol_lg(z, lx);
1072 } else {
1073 z = cgetg(ly,t_POL); z[1] = y[1];
1074 for (i=2; i < lx; i++) gel(z,i) = gsub(gel(x,i),gel(y,i));
1075 for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
1076 z = normalizepol_lg(z, ly);
1077 }
1078 return z;
1079 }
1080 GEN
1081 RgX_neg(GEN x)
1082 {
1083 long i, lx = lg(x);
1084 GEN y = cgetg(lx, t_POL); y[1] = x[1];
1085 for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
1086 return y;
1087 }
1088
1089 GEN
1090 RgX_Rg_add(GEN y, GEN x)
1091 {
1092 GEN z;
1093 long lz = lg(y), i;
1094 if (lz == 2) return scalarpol(x,varn(y));
1095 z = cgetg(lz,t_POL); z[1] = y[1];
1096 gel(z,2) = gadd(gel(y,2),x);
1097 for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1098 /* probably useless unless lz = 3, but cannot be skipped if y is
1099 * an inexact 0 */
1100 return normalizepol_lg(z,lz);
1101 }
1102 GEN
1103 RgX_Rg_add_shallow(GEN y, GEN x)
1104 {
1105 GEN z;
1106 long lz = lg(y), i;
1107 if (lz == 2) return scalarpol(x,varn(y));
1108 z = cgetg(lz,t_POL); z[1] = y[1];
1109 gel(z,2) = gadd(gel(y,2),x);
1110 for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
1111 return z = normalizepol_lg(z,lz);
1112 }
1113 GEN
1114 RgX_Rg_sub(GEN y, GEN x)
1115 {
1116 GEN z;
1117 long lz = lg(y), i;
1118 if (lz == 2)
1119 { /* scalarpol(gneg(x),varn(y)) optimized */
1120 long v = varn(y);
1121 if (isrationalzero(x)) return pol_0(v);
1122 z = cgetg(3,t_POL);
1123 z[1] = gequal0(x)? evalvarn(v)
1124 : evalvarn(v) | evalsigne(1);
1125 gel(z,2) = gneg(x); return z;
1126 }
1127 z = cgetg(lz,t_POL); z[1] = y[1];
1128 gel(z,2) = gsub(gel(y,2),x);
1129 for(i=3; i<lz; i++) gel(z,i) = gcopy(gel(y,i));
1130 return z = normalizepol_lg(z,lz);
1131 }
1132 GEN
1133 Rg_RgX_sub(GEN x, GEN y)
1134 {
1135 GEN z;
1136 long lz = lg(y), i;
1137 if (lz == 2) return scalarpol(x,varn(y));
1138 z = cgetg(lz,t_POL); z[1] = y[1];
1139 gel(z,2) = gsub(x, gel(y,2));
1140 for(i=3; i<lz; i++) gel(z,i) = gneg(gel(y,i));
1141 return z = normalizepol_lg(z,lz);
1142 }
1143 /*******************************************************************/
1144 /* */
1145 /* KARATSUBA MULTIPLICATION */
1146 /* */
1147 /*******************************************************************/
1148 #if 0
1149 /* to debug Karatsuba-like routines */
1150 GEN
1151 zx_debug_spec(GEN x, long nx)
1152 {
1153 GEN z = cgetg(nx+2,t_POL);
1154 long i;
1155 for (i=0; i<nx; i++) gel(z,i+2) = stoi(x[i]);
1156 z[1] = evalsigne(1); return z;
1157 }
1158
1159 GEN
1160 RgX_debug_spec(GEN x, long nx)
1161 {
1162 GEN z = cgetg(nx+2,t_POL);
1163 long i;
1164 for (i=0; i<nx; i++) z[i+2] = x[i];
1165 z[1] = evalsigne(1); return z;
1166 }
1167 #endif
1168
1169 /* generic multiplication */
1170 GEN
1171 RgX_addspec_shallow(GEN x, GEN y, long nx, long ny)
1172 {
1173 GEN z, t;
1174 long i;
1175 if (nx == ny) {
1176 z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1177 for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1178 return normalizepol_lg(z, nx+2);
1179 }
1180 if (ny < nx) {
1181 z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1182 for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1183 for ( ; i < nx; i++) gel(t,i) = gel(x,i);
1184 return normalizepol_lg(z, nx+2);
1185 } else {
1186 z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1187 for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1188 for ( ; i < ny; i++) gel(t,i) = gel(y,i);
1189 return normalizepol_lg(z, ny+2);
1190 }
1191 }
1192 GEN
1193 RgX_addspec(GEN x, GEN y, long nx, long ny)
1194 {
1195 GEN z, t;
1196 long i;
1197 if (nx == ny) {
1198 z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1199 for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1200 return normalizepol_lg(z, nx+2);
1201 }
1202 if (ny < nx) {
1203 z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1204 for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1205 for ( ; i < nx; i++) gel(t,i) = gcopy(gel(x,i));
1206 return normalizepol_lg(z, nx+2);
1207 } else {
1208 z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
1209 for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
1210 for ( ; i < ny; i++) gel(t,i) = gcopy(gel(y,i));
1211 return normalizepol_lg(z, ny+2);
1212 }
1213 }
1214
1215 /* Return the vector of coefficients of x, where we replace rational 0s by NULL
1216 * [ to speed up basic operation s += x[i]*y[j] ]. We create a proper
1217 * t_VECSMALL, to hold this, which can be left on stack: gerepile
1218 * will not crash on it. The returned vector itself is not a proper GEN,
1219 * we access the coefficients as x[i], i = 0..deg(x) */
1220 static GEN
1221 RgXspec_kill0(GEN x, long lx)
1222 {
1223 GEN z = cgetg(lx+1, t_VECSMALL) + 1; /* inhibit gerepile-wise */
1224 long i;
1225 for (i=0; i <lx; i++)
1226 {
1227 GEN c = gel(x,i);
1228 z[i] = (long)(isrationalzero(c)? NULL: c);
1229 }
1230 return z;
1231 }
1232
1233 INLINE GEN
1234 RgX_mulspec_basecase_limb(GEN x, GEN y, long a, long b)
1235 {
1236 pari_sp av = avma;
1237 GEN s = NULL;
1238 long i;
1239
1240 for (i=a; i<b; i++)
1241 if (gel(y,i) && gel(x,-i))
1242 {
1243 GEN t = gmul(gel(y,i), gel(x,-i));
1244 s = s? gadd(s, t): t;
1245 }
1246 return s? gerepileupto(av, s): gen_0;
1247 }
1248
1249 /* assume nx >= ny > 0, return x * y * t^v */
1250 static GEN
1251 RgX_mulspec_basecase(GEN x, GEN y, long nx, long ny, long v)
1252 {
1253 long i, lz, nz;
1254 GEN z;
1255
1256 x = RgXspec_kill0(x,nx);
1257 y = RgXspec_kill0(y,ny);
1258 lz = nx + ny + 1; nz = lz-2;
1259 lz += v;
1260 z = cgetg(lz, t_POL) + 2; /* x:y:z [i] = term of degree i */
1261 for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1262 for (i=0; i<ny; i++)gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0, i+1);
1263 for ( ; i<nx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ny);
1264 for ( ; i<nz; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-nx+1,ny);
1265 z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1266 }
1267
1268 /* return (x * X^d) + y. Assume d > 0 */
1269 GEN
1270 RgX_addmulXn_shallow(GEN x0, GEN y0, long d)
1271 {
1272 GEN x, y, xd, yd, zd;
1273 long a, lz, nx, ny;
1274
1275 if (!signe(x0)) return y0;
1276 ny = lgpol(y0);
1277 nx = lgpol(x0);
1278 zd = (GEN)avma;
1279 x = x0 + 2; y = y0 + 2; a = ny-d;
1280 if (a <= 0)
1281 {
1282 lz = nx+d+2;
1283 (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1284 while (xd > x) gel(--zd,0) = gel(--xd,0);
1285 x = zd + a;
1286 while (zd > x) gel(--zd,0) = gen_0;
1287 }
1288 else
1289 {
1290 xd = new_chunk(d); yd = y+d;
1291 x = RgX_addspec_shallow(x,yd, nx,a);
1292 lz = (a>nx)? ny+2: lg(x)+d;
1293 x += 2; while (xd > x) *--zd = *--xd;
1294 }
1295 while (yd > y) *--zd = *--yd;
1296 *--zd = x0[1];
1297 *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1298 }
1299 GEN
1300 RgX_addmulXn(GEN x0, GEN y0, long d)
1301 {
1302 GEN x, y, xd, yd, zd;
1303 long a, lz, nx, ny;
1304
1305 if (!signe(x0)) return RgX_copy(y0);
1306 nx = lgpol(x0);
1307 ny = lgpol(y0);
1308 zd = (GEN)avma;
1309 x = x0 + 2; y = y0 + 2; a = ny-d;
1310 if (a <= 0)
1311 {
1312 lz = nx+d+2;
1313 (void)new_chunk(lz); xd = x+nx; yd = y+ny;
1314 while (xd > x) gel(--zd,0) = gcopy(gel(--xd,0));
1315 x = zd + a;
1316 while (zd > x) gel(--zd,0) = gen_0;
1317 }
1318 else
1319 {
1320 xd = new_chunk(d); yd = y+d;
1321 x = RgX_addspec(x,yd, nx,a);
1322 lz = (a>nx)? ny+2: lg(x)+d;
1323 x += 2; while (xd > x) *--zd = *--xd;
1324 }
1325 while (yd > y) gel(--zd,0) = gcopy(gel(--yd,0));
1326 *--zd = x0[1];
1327 *--zd = evaltyp(t_POL) | evallg(lz); return zd;
1328 }
1329
1330 /* return x * y mod t^n */
1331 static GEN
1332 RgXn_mul_basecase(GEN x, GEN y, long n)
1333 {
1334 long i, lz = n+2, lx = lgpol(x), ly = lgpol(y);
1335 GEN z;
1336 if (lx < 0) return pol_0(varn(x));
1337 if (ly < 0) return pol_0(varn(x));
1338 z = cgetg(lz, t_POL) + 2;
1339 x+=2; if (lx > n) lx = n;
1340 y+=2; if (ly > n) ly = n;
1341 z[-1] = x[-1];
1342 if (ly > lx) { swap(x,y); lswap(lx,ly); }
1343 x = RgXspec_kill0(x, lx);
1344 y = RgXspec_kill0(y, ly);
1345 /* x:y:z [i] = term of degree i */
1346 for (i=0;i<ly; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,i+1);
1347 for ( ; i<lx; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, 0,ly);
1348 for ( ; i<n; i++) gel(z,i) = RgX_mulspec_basecase_limb(x+i,y, i-lx+1,ly);
1349 return normalizepol_lg(z - 2, lz);
1350 }
1351 /* Mulders / Karatsuba product f*g mod t^n (Hanrot-Zimmermann variant) */
1352 static GEN
1353 RgXn_mul2(GEN f, GEN g, long n)
1354 {
1355 pari_sp av = avma;
1356 GEN fe,fo, ge,go, l,h,m;
1357 long n0, n1;
1358 if (degpol(f) + degpol(g) < n) return RgX_mul(f,g);
1359 if (n < 80) return RgXn_mul_basecase(f,g,n);
1360 n0 = n>>1; n1 = n-n0;
1361 RgX_even_odd(f, &fe, &fo);
1362 RgX_even_odd(g, &ge, &go);
1363 l = RgXn_mul(fe,ge,n1);
1364 h = RgXn_mul(fo,go,n0);
1365 m = RgX_sub(RgXn_mul(RgX_add(fe,fo),RgX_add(ge,go),n0), RgX_add(l,h));
1366 /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1367 * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1368 l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1369 /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1370 if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1371 m = RgX_inflate(m,2);
1372 /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1373 if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1374 h = RgX_inflate(h,2);
1375 h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1376 return gerepileupto(av, h);
1377 }
1378 /* (f*g) \/ x^n */
1379 static GEN
1380 RgX_mulhigh_i2(GEN f, GEN g, long n)
1381 {
1382 long d = degpol(f)+degpol(g) + 1 - n;
1383 GEN h;
1384 if (d <= 2) return RgX_shift_shallow(RgX_mul(f,g), -n);
1385 h = RgX_recip_shallow(RgXn_mul(RgX_recip_shallow(f),
1386 RgX_recip_shallow(g), d));
1387 return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1388 }
1389
1390 /* (f*g) \/ x^n */
1391 static GEN
1392 RgX_sqrhigh_i2(GEN f, long n)
1393 {
1394 long d = 2*degpol(f)+ 1 - n;
1395 GEN h;
1396 if (d <= 2) return RgX_shift_shallow(RgX_sqr(f), -n);
1397 h = RgX_recip_shallow(RgXn_sqr(RgX_recip_shallow(f), d));
1398 return RgX_shift_shallow(h, d-1-degpol(h)); /* possibly (fg)(0) = 0 */
1399 }
1400
1401 /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
1402 * b+2 were sent instead. na, nb = number of terms of a, b.
1403 * Only c, c0, c1, c2 are genuine GEN.
1404 */
1405 GEN
1406 RgX_mulspec(GEN a, GEN b, long na, long nb)
1407 {
1408 GEN a0, c, c0;
1409 long n0, n0a, i, v = 0;
1410 pari_sp av;
1411
1412 while (na && isrationalzero(gel(a,0))) { a++; na--; v++; }
1413 while (nb && isrationalzero(gel(b,0))) { b++; nb--; v++; }
1414 if (na < nb) swapspec(a,b, na,nb);
1415 if (!nb) return pol_0(0);
1416
1417 if (nb < RgX_MUL_LIMIT) return RgX_mulspec_basecase(a,b,na,nb, v);
1418 RgX_shift_inplace_init(v);
1419 i = (na>>1); n0 = na-i; na = i;
1420 av = avma; a0 = a+n0; n0a = n0;
1421 while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1422
1423 if (nb > n0)
1424 {
1425 GEN b0,c1,c2;
1426 long n0b;
1427
1428 nb -= n0; b0 = b+n0; n0b = n0;
1429 while (n0b && isrationalzero(gel(b,n0b-1))) n0b--;
1430 c = RgX_mulspec(a,b,n0a,n0b);
1431 c0 = RgX_mulspec(a0,b0, na,nb);
1432
1433 c2 = RgX_addspec_shallow(a0,a, na,n0a);
1434 c1 = RgX_addspec_shallow(b0,b, nb,n0b);
1435
1436 c1 = RgX_mulspec(c1+2,c2+2, lgpol(c1),lgpol(c2));
1437 c2 = RgX_sub(c1, RgX_add(c0,c));
1438 c0 = RgX_addmulXn_shallow(c0, c2, n0);
1439 }
1440 else
1441 {
1442 c = RgX_mulspec(a,b,n0a,nb);
1443 c0 = RgX_mulspec(a0,b,na,nb);
1444 }
1445 c0 = RgX_addmulXn(c0,c,n0);
1446 return RgX_shift_inplace(gerepileupto(av,c0), v);
1447 }
1448
1449 INLINE GEN
1450 RgX_sqrspec_basecase_limb(GEN x, long a, long i)
1451 {
1452 pari_sp av = avma;
1453 GEN s = NULL;
1454 long j, l = (i+1)>>1;
1455 for (j=a; j<l; j++)
1456 {
1457 GEN xj = gel(x,j), xx = gel(x,i-j);
1458 if (xj && xx)
1459 {
1460 GEN t = gmul(xj, xx);
1461 s = s? gadd(s, t): t;
1462 }
1463 }
1464 if (s) s = gshift(s,1);
1465 if ((i&1) == 0)
1466 {
1467 GEN t = gel(x, i>>1);
1468 if (t) {
1469 t = gsqr(t);
1470 s = s? gadd(s, t): t;
1471 }
1472 }
1473 return s? gerepileupto(av,s): gen_0;
1474 }
1475 static GEN
1476 RgX_sqrspec_basecase(GEN x, long nx, long v)
1477 {
1478 long i, lz, nz;
1479 GEN z;
1480
1481 if (!nx) return pol_0(0);
1482 x = RgXspec_kill0(x,nx);
1483 lz = (nx << 1) + 1, nz = lz-2;
1484 lz += v;
1485 z = cgetg(lz,t_POL) + 2;
1486 for (i=0; i<v; i++) gel(z++, 0) = gen_0;
1487 for (i=0; i<nx; i++)gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1488 for ( ; i<nz; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-nx+1, i);
1489 z -= v+2; z[1] = 0; return normalizepol_lg(z, lz);
1490 }
1491 /* return x^2 mod t^n */
1492 static GEN
1493 RgXn_sqr_basecase(GEN x, long n)
1494 {
1495 long i, lz = n+2, lx = lgpol(x);
1496 GEN z;
1497 if (lx < 0) return pol_0(varn(x));
1498 z = cgetg(lz, t_POL);
1499 z[1] = x[1];
1500 x+=2; if (lx > n) lx = n;
1501 x = RgXspec_kill0(x,lx);
1502 z+=2;/* x:z [i] = term of degree i */
1503 for (i=0;i<lx; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, 0, i);
1504 for ( ; i<n; i++) gel(z,i) = RgX_sqrspec_basecase_limb(x, i-lx+1, i);
1505 z -= 2; return normalizepol_lg(z, lz);
1506 }
1507 /* Mulders / Karatsuba product f^2 mod t^n (Hanrot-Zimmermann variant) */
1508 static GEN
1509 RgXn_sqr2(GEN f, long n)
1510 {
1511 pari_sp av = avma;
1512 GEN fe,fo, l,h,m;
1513 long n0, n1;
1514 if (2*degpol(f) < n) return RgX_sqr_i(f);
1515 if (n < 80) return RgXn_sqr_basecase(f,n);
1516 n0 = n>>1; n1 = n-n0;
1517 RgX_even_odd(f, &fe, &fo);
1518 l = RgXn_sqr(fe,n1);
1519 h = RgXn_sqr(fo,n0);
1520 m = RgX_sub(RgXn_sqr(RgX_add(fe,fo),n0), RgX_add(l,h));
1521 /* n1-1 <= n0 <= n1, deg l,m <= n1-1, deg h <= n0-1
1522 * result is t^2 h(t^2) + t m(t^2) + l(t^2) */
1523 l = RgX_inflate(l,2); /* deg l <= 2n1 - 2 <= n-1 */
1524 /* deg(t m(t^2)) <= 2n1 - 1 <= n, truncate to < n */
1525 if (2*degpol(m)+1 == n) m = normalizepol_lg(m, lg(m)-1);
1526 m = RgX_inflate(m,2);
1527 /* deg(t^2 h(t^2)) <= 2n0 <= n, truncate to < n */
1528 if (2*degpol(h)+2 == n) h = normalizepol_lg(h, lg(h)-1);
1529 h = RgX_inflate(h,2);
1530 h = RgX_addmulXn(RgX_addmulXn_shallow(h,m,1), l,1);
1531 return gerepileupto(av, h);
1532 }
1533 GEN
1534 RgX_sqrspec(GEN a, long na)
1535 {
1536 GEN a0, c, c0, c1;
1537 long n0, n0a, i, v = 0;
1538 pari_sp av;
1539
1540 while (na && isrationalzero(gel(a,0))) { a++; na--; v += 2; }
1541 if (na<RgX_SQR_LIMIT) return RgX_sqrspec_basecase(a, na, v);
1542 RgX_shift_inplace_init(v);
1543 i = (na>>1); n0 = na-i; na = i;
1544 av = avma; a0 = a+n0; n0a = n0;
1545 while (n0a && isrationalzero(gel(a,n0a-1))) n0a--;
1546
1547 c = RgX_sqrspec(a,n0a);
1548 c0 = RgX_sqrspec(a0,na);
1549 c1 = gmul2n(RgX_mulspec(a0,a, na,n0a), 1);
1550 c0 = RgX_addmulXn_shallow(c0,c1, n0);
1551 c0 = RgX_addmulXn(c0,c,n0);
1552 return RgX_shift_inplace(gerepileupto(av,c0), v);
1553 }
1554
1555 /* (X^a + A)(X^b + B) - X^(a+b), where deg A < a, deg B < b */
1556 GEN
1557 RgX_mul_normalized(GEN A, long a, GEN B, long b)
1558 {
1559 GEN z = RgX_mul(A, B);
1560 if (a < b)
1561 z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(A, B, b-a), z, a);
1562 else if (a > b)
1563 z = RgX_addmulXn_shallow(RgX_addmulXn_shallow(B, A, a-b), z, b);
1564 else
1565 z = RgX_addmulXn_shallow(RgX_add(A, B), z, a);
1566 return z;
1567 }
1568
1569 GEN
1570 RgX_mul_i(GEN x, GEN y)
1571 {
1572 GEN z = RgX_mulspec(x+2, y+2, lgpol(x), lgpol(y));
1573 setvarn(z, varn(x)); return z;
1574 }
1575
1576 GEN
1577 RgX_sqr_i(GEN x)
1578 {
1579 GEN z = RgX_sqrspec(x+2, lgpol(x));
1580 setvarn(z,varn(x)); return z;
1581 }
1582
1583 /*******************************************************************/
1584 /* */
1585 /* DIVISION */
1586 /* */
1587 /*******************************************************************/
1588 GEN
1589 RgX_Rg_divexact(GEN x, GEN y) {
1590 long i, lx = lg(x);
1591 GEN z;
1592 if (lx == 2) return gcopy(x);
1593 switch(typ(y))
1594 {
1595 case t_INT:
1596 if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1597 break;
1598 case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1599 }
1600 z = cgetg(lx, t_POL); z[1] = x[1];
1601 for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1602 return z;
1603 }
1604 GEN
1605 RgX_Rg_div(GEN x, GEN y) {
1606 long i, lx = lg(x);
1607 GEN z;
1608 if (lx == 2) return gcopy(x);
1609 switch(typ(y))
1610 {
1611 case t_INT:
1612 if (is_pm1(y)) return signe(y) < 0 ? RgX_neg(x): RgX_copy(x);
1613 break;
1614 case t_INTMOD: case t_POLMOD: return RgX_Rg_mul(x, ginv(y));
1615 }
1616 z = cgetg(lx, t_POL); z[1] = x[1];
1617 for (i=2; i<lx; i++) gel(z,i) = gdiv(gel(x,i),y);
1618 return normalizepol_lg(z, lx);
1619 }
1620 GEN
1621 RgX_normalize(GEN x)
1622 {
1623 GEN z, d = NULL;
1624 long i, n = lg(x)-1;
1625 for (i = n; i > 1; i--) { d = gel(x,i); if (!gequal0(d)) break; }
1626 if (i == 1) return pol_0(varn(x));
1627 if (i == n && isint1(d)) return x;
1628 n = i; z = cgetg(n+1, t_POL); z[1] = x[1];
1629 for (i=2; i<n; i++) gel(z,i) = gdiv(gel(x,i),d);
1630 gel(z,n) = Rg_get_1(d); return z;
1631 }
1632 GEN
1633 RgX_divs(GEN x, long y) {
1634 long i, lx;
1635 GEN z = cgetg_copy(x, &lx); z[1] = x[1];
1636 for (i=2; i<lx; i++) gel(z,i) = gdivgs(gel(x,i),y);
1637 return normalizepol_lg(z, lx);
1638 }
1639 GEN
1640 RgX_div_by_X_x(GEN a, GEN x, GEN *r)
1641 {
1642 long l = lg(a), i;
1643 GEN a0, z0, z;
1644
1645 if (l <= 3)
1646 {
1647 if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1648 return pol_0(0);
1649 }
1650 z = cgetg(l-1, t_POL);
1651 z[1] = a[1];
1652 a0 = a + l-1;
1653 z0 = z + l-2; *z0 = *a0--;
1654 for (i=l-3; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1655 {
1656 GEN t = gadd(gel(a0--,0), gmul(x, gel(z0--,0)));
1657 gel(z0,0) = t;
1658 }
1659 if (r) *r = gadd(gel(a0,0), gmul(x, gel(z0,0)));
1660 return z;
1661 }
1662 /* Polynomial division x / y:
1663 * if pr = ONLY_REM return remainder, otherwise return quotient
1664 * if pr = ONLY_DIVIDES return quotient if division is exact, else NULL
1665 * if pr != NULL set *pr to remainder, as the last object on stack */
1666 /* assume, typ(x) = typ(y) = t_POL, same variable */
1667 static GEN
1668 RgX_divrem_i(GEN x, GEN y, GEN *pr)
1669 {
1670 pari_sp avy, av, av1;
1671 long dx,dy,dz,i,j,sx,lr;
1672 GEN z,p1,p2,rem,y_lead,mod,p;
1673 GEN (*f)(GEN,GEN);
1674
1675 if (!signe(y)) pari_err_INV("RgX_divrem",y);
1676
1677 dy = degpol(y);
1678 y_lead = gel(y,dy+2);
1679 if (gequal0(y_lead)) /* normalize denominator if leading term is 0 */
1680 {
1681 pari_warn(warner,"normalizing a polynomial with 0 leading term");
1682 for (dy--; dy>=0; dy--)
1683 {
1684 y_lead = gel(y,dy+2);
1685 if (!gequal0(y_lead)) break;
1686 }
1687 }
1688 if (!dy) /* y is constant */
1689 {
1690 if (pr == ONLY_REM) return pol_0(varn(x));
1691 z = RgX_Rg_div(x, y_lead);
1692 if (pr == ONLY_DIVIDES) return z;
1693 if (pr) *pr = pol_0(varn(x));
1694 return z;
1695 }
1696 dx = degpol(x);
1697 if (dx < dy)
1698 {
1699 if (pr == ONLY_REM) return RgX_copy(x);
1700 if (pr == ONLY_DIVIDES) return signe(x)? NULL: pol_0(varn(x));
1701 z = pol_0(varn(x));
1702 if (pr) *pr = RgX_copy(x);
1703 return z;
1704 }
1705
1706 /* x,y in R[X], y non constant */
1707 av = avma;
1708 p = NULL;
1709 if (RgX_is_FpX(x, &p) && RgX_is_FpX(y, &p) && p)
1710 {
1711 z = FpX_divrem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, pr);
1712 if (!z) return gc_NULL(av);
1713 z = FpX_to_mod(z, p);
1714 if (!pr || pr == ONLY_REM || pr == ONLY_DIVIDES)
1715 return gerepileupto(av, z);
1716 *pr = FpX_to_mod(*pr, p);
1717 gerepileall(av, 2, pr, &z);
1718 return z;
1719 }
1720 switch(typ(y_lead))
1721 {
1722 case t_REAL:
1723 y_lead = ginv(y_lead);
1724 f = gmul; mod = NULL;
1725 break;
1726 case t_INTMOD:
1727 case t_POLMOD: y_lead = ginv(y_lead);
1728 f = gmul; mod = gmodulo(gen_1, gel(y_lead,1));
1729 break;
1730 default: if (gequal1(y_lead)) y_lead = NULL;
1731 f = gdiv; mod = NULL;
1732 }
1733
1734 if (y_lead == NULL)
1735 p2 = gel(x,dx+2);
1736 else {
1737 for(;;) {
1738 p2 = f(gel(x,dx+2),y_lead);
1739 p2 = simplify_shallow(p2);
1740 if (!isexactzero(p2) || (--dx < 0)) break;
1741 }
1742 if (dx < dy) /* leading coeff of x was in fact zero */
1743 {
1744 if (pr == ONLY_DIVIDES) {
1745 set_avma(av);
1746 return (dx < 0)? pol_0(varn(x)) : NULL;
1747 }
1748 if (pr == ONLY_REM)
1749 {
1750 if (dx < 0)
1751 return gerepilecopy(av, scalarpol(p2, varn(x)));
1752 else
1753 {
1754 GEN t;
1755 set_avma(av);
1756 t = cgetg(dx + 3, t_POL); t[1] = x[1];
1757 for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1758 return t;
1759 }
1760 }
1761 if (pr) /* cf ONLY_REM above */
1762 {
1763 if (dx < 0)
1764 {
1765 p2 = gclone(p2);
1766 set_avma(av);
1767 z = pol_0(varn(x));
1768 x = scalarpol(p2, varn(x));
1769 gunclone(p2);
1770 }
1771 else
1772 {
1773 GEN t;
1774 set_avma(av);
1775 z = pol_0(varn(x));
1776 t = cgetg(dx + 3, t_POL); t[1] = x[1];
1777 for (i = 2; i < dx + 3; i++) gel(t,i) = gcopy(gel(x,i));
1778 x = t;
1779 }
1780 *pr = x;
1781 }
1782 else
1783 {
1784 set_avma(av);
1785 z = pol_0(varn(x));
1786 }
1787 return z;
1788 }
1789 }
1790 /* dx >= dy */
1791 avy = avma;
1792 dz = dx-dy;
1793 z = cgetg(dz+3,t_POL); z[1] = x[1];
1794 x += 2;
1795 z += 2;
1796 y += 2;
1797 gel(z,dz) = gcopy(p2);
1798
1799 for (i=dx-1; i>=dy; i--)
1800 {
1801 av1=avma; p1=gel(x,i);
1802 for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1803 if (y_lead) p1 = simplify(f(p1,y_lead));
1804
1805 if (isrationalzero(p1)) { set_avma(av1); p1 = gen_0; }
1806 else
1807 p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);
1808 gel(z,i-dy) = p1;
1809 }
1810 if (!pr) return gerepileupto(av,z-2);
1811
1812 rem = (GEN)avma; av1 = (pari_sp)new_chunk(dx+3);
1813 for (sx=0; ; i--)
1814 {
1815 p1 = gel(x,i);
1816 /* we always enter this loop at least once */
1817 for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1818 if (mod && avma==av1) p1 = gmul(p1,mod);
1819 if (!gequal0(p1)) { sx = 1; break; } /* remainder is nonzero */
1820 if (!isexactzero(p1)) break;
1821 if (!i) break;
1822 set_avma(av1);
1823 }
1824 if (pr == ONLY_DIVIDES)
1825 {
1826 if (sx) return gc_NULL(av);
1827 set_avma((pari_sp)rem); return gerepileupto(av,z-2);
1828 }
1829 lr=i+3; rem -= lr;
1830 if (avma==av1) { set_avma((pari_sp)rem); p1 = gcopy(p1); }
1831 else p1 = gerepileupto((pari_sp)rem,p1);
1832 rem[0] = evaltyp(t_POL) | evallg(lr);
1833 rem[1] = z[-1];
1834 rem += 2;
1835 gel(rem,i) = p1;
1836 for (i--; i>=0; i--)
1837 {
1838 av1=avma; p1 = gel(x,i);
1839 for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1840 if (mod && avma==av1) p1 = gmul(p1,mod);
1841 gel(rem,i) = avma==av1? gcopy(p1):gerepileupto(av1,p1);
1842 }
1843 rem -= 2;
1844 if (!sx) (void)normalizepol_lg(rem, lr);
1845 if (pr == ONLY_REM) return gerepileupto(av,rem);
1846 z -= 2;
1847 {
1848 GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;
1849 gerepilemanysp(av,avy,gptr,2); *pr = rem; return z;
1850 }
1851 }
1852
1853 GEN
1854 RgX_divrem(GEN x, GEN y, GEN *pr)
1855 {
1856 if (pr == ONLY_REM) return RgX_rem(x, y);
1857 return RgX_divrem_i(x, y, pr);
1858 }
1859
1860 /* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */
1861 GEN
1862 RgXQX_divrem(GEN x, GEN y, GEN T, GEN *pr)
1863 {
1864 long vx, dx, dy, dz, i, j, sx, lr;
1865 pari_sp av0, av, tetpil;
1866 GEN z,p1,rem,lead;
1867
1868 if (!signe(y)) pari_err_INV("RgXQX_divrem",y);
1869 vx = varn(x);
1870 dx = degpol(x);
1871 dy = degpol(y);
1872 if (dx < dy)
1873 {
1874 if (pr)
1875 {
1876 av0 = avma; x = RgXQX_red(x, T);
1877 if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: gen_0; }
1878 if (pr == ONLY_REM) return x;
1879 *pr = x;
1880 }
1881 return pol_0(vx);
1882 }
1883 lead = leading_coeff(y);
1884 if (!dy) /* y is constant */
1885 {
1886 if (pr && pr != ONLY_DIVIDES)
1887 {
1888 if (pr == ONLY_REM) return pol_0(vx);
1889 *pr = pol_0(vx);
1890 }
1891 if (gequal1(lead)) return RgX_copy(x);
1892 av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma;
1893 return gerepile(av0,tetpil,RgXQX_red(x,T));
1894 }
1895 av0 = avma; dz = dx-dy;
1896 lead = gequal1(lead)? NULL: gclone(ginvmod(lead,T));
1897 set_avma(av0);
1898 z = cgetg(dz+3,t_POL); z[1] = x[1];
1899 x += 2; y += 2; z += 2;
1900
1901 p1 = gel(x,dx); av = avma;
1902 gel(z,dz) = lead? gerepileupto(av, grem(gmul(p1,lead), T)): gcopy(p1);
1903 for (i=dx-1; i>=dy; i--)
1904 {
1905 av=avma; p1=gel(x,i);
1906 for (j=i-dy+1; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1907 if (lead) p1 = gmul(grem(p1, T), lead);
1908 tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil, grem(p1, T));
1909 }
1910 if (!pr) { guncloneNULL(lead); return z-2; }
1911
1912 rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
1913 for (sx=0; ; i--)
1914 {
1915 p1 = gel(x,i);
1916 for (j=0; j<=i && j<=dz; j++) p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1917 tetpil=avma; p1 = grem(p1, T); if (!gequal0(p1)) { sx = 1; break; }
1918 if (!i) break;
1919 set_avma(av);
1920 }
1921 if (pr == ONLY_DIVIDES)
1922 {
1923 guncloneNULL(lead);
1924 if (sx) return gc_NULL(av0);
1925 return gc_const((pari_sp)rem, z-2);
1926 }
1927 lr=i+3; rem -= lr;
1928 rem[0] = evaltyp(t_POL) | evallg(lr);
1929 rem[1] = z[-1];
1930 p1 = gerepile((pari_sp)rem,tetpil,p1);
1931 rem += 2; gel(rem,i) = p1;
1932 for (i--; i>=0; i--)
1933 {
1934 av=avma; p1 = gel(x,i);
1935 for (j=0; j<=i && j<=dz; j++)
1936 p1 = gsub(p1, gmul(gel(z,j),gel(y,i-j)));
1937 tetpil=avma; gel(rem,i) = gerepile(av,tetpil, grem(p1, T));
1938 }
1939 rem -= 2;
1940 guncloneNULL(lead);
1941 if (!sx) (void)normalizepol_lg(rem, lr);
1942 if (pr == ONLY_REM) return gerepileupto(av0,rem);
1943 *pr = rem; return z-2;
1944 }
1945
1946 /*******************************************************************/
1947 /* */
1948 /* PSEUDO-DIVISION */
1949 /* */
1950 /*******************************************************************/
1951 INLINE GEN
1952 rem(GEN c, GEN T)
1953 {
1954 if (T && typ(c) == t_POL && varn(c) == varn(T)) c = RgX_rem(c, T);
1955 return c;
1956 }
1957
1958 /* x, y, are ZYX, lc(y) is an integer, T is a ZY */
1959 int
1960 ZXQX_dvd(GEN x, GEN y, GEN T)
1961 {
1962 long dx, dy, dz, i, p, T_ismonic;
1963 pari_sp av = avma, av2;
1964 GEN y_lead;
1965
1966 if (!signe(y)) pari_err_INV("ZXQX_dvd",y);
1967 dy = degpol(y); y_lead = gel(y,dy+2);
1968 if (typ(y_lead) == t_POL) y_lead = gel(y_lead, 2); /* t_INT */
1969 /* if monic, no point in using pseudo-division */
1970 if (gequal1(y_lead)) return signe(RgXQX_rem(x, y, T)) == 0;
1971 T_ismonic = gequal1(leading_coeff(T));
1972 dx = degpol(x);
1973 if (dx < dy) return !signe(x);
1974 (void)new_chunk(2);
1975 x = RgX_recip_shallow(x)+2;
1976 y = RgX_recip_shallow(y)+2;
1977 /* pay attention to sparse divisors */
1978 for (i = 1; i <= dy; i++)
1979 if (!signe(gel(y,i))) gel(y,i) = NULL;
1980 dz = dx-dy; p = dz+1;
1981 av2 = avma;
1982 for (;;)
1983 {
1984 GEN m, x0 = gel(x,0), y0 = y_lead, cx = content(x0);
1985 x0 = gneg(x0); p--;
1986 m = gcdii(cx, y0);
1987 if (!equali1(m))
1988 {
1989 x0 = gdiv(x0, m);
1990 y0 = diviiexact(y0, m);
1991 if (equali1(y0)) y0 = NULL;
1992 }
1993 for (i=1; i<=dy; i++)
1994 {
1995 GEN c = gel(x,i); if (y0) c = gmul(y0, c);
1996 if (gel(y,i)) c = gadd(c, gmul(x0,gel(y,i)));
1997 if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
1998 gel(x,i) = c;
1999 }
2000 for ( ; i<=dx; i++)
2001 {
2002 GEN c = gel(x,i); if (y0) c = gmul(y0, c);
2003 if (typ(c) == t_POL) c = T_ismonic ? ZX_rem(c, T): RgX_rem(c, T);
2004 gel(x,i) = c;
2005 }
2006 do { x++; dx--; } while (dx >= 0 && !signe(gel(x,0)));
2007 if (dx < dy) break;
2008 if (gc_needed(av2,1))
2009 {
2010 if(DEBUGMEM>1) pari_warn(warnmem,"ZXQX_dvd dx = %ld >= %ld",dx,dy);
2011 gerepilecoeffs(av2,x,dx+1);
2012 }
2013 }
2014 return gc_bool(av, dx < 0);
2015 }
2016
2017 /* T either NULL or a t_POL. */
2018 GEN
2019 RgXQX_pseudorem(GEN x, GEN y, GEN T)
2020 {
2021 long vx = varn(x), dx, dy, dz, i, lx, p;
2022 pari_sp av = avma, av2;
2023 GEN y_lead;
2024
2025 if (!signe(y)) pari_err_INV("RgXQX_pseudorem",y);
2026 dy = degpol(y); y_lead = gel(y,dy+2);
2027 /* if monic, no point in using pseudo-division */
2028 if (gequal1(y_lead)) return T? RgXQX_rem(x, y, T): RgX_rem(x, y);
2029 dx = degpol(x);
2030 if (dx < dy) return RgX_copy(x);
2031 (void)new_chunk(2);
2032 x = RgX_recip_shallow(x)+2;
2033 y = RgX_recip_shallow(y)+2;
2034 /* pay attention to sparse divisors */
2035 for (i = 1; i <= dy; i++)
2036 if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2037 dz = dx-dy; p = dz+1;
2038 av2 = avma;
2039 for (;;)
2040 {
2041 gel(x,0) = gneg(gel(x,0)); p--;
2042 for (i=1; i<=dy; i++)
2043 {
2044 GEN c = gmul(y_lead, gel(x,i));
2045 if (gel(y,i)) c = gadd(c, gmul(gel(x,0),gel(y,i)));
2046 gel(x,i) = rem(c, T);
2047 }
2048 for ( ; i<=dx; i++)
2049 {
2050 GEN c = gmul(y_lead, gel(x,i));
2051 gel(x,i) = rem(c, T);
2052 }
2053 do { x++; dx--; } while (dx >= 0 && gequal0(gel(x,0)));
2054 if (dx < dy) break;
2055 if (gc_needed(av2,1))
2056 {
2057 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudorem dx = %ld >= %ld",dx,dy);
2058 gerepilecoeffs(av2,x,dx+1);
2059 }
2060 }
2061 if (dx < 0) return pol_0(vx);
2062 lx = dx+3; x -= 2;
2063 x[0] = evaltyp(t_POL) | evallg(lx);
2064 x[1] = evalsigne(1) | evalvarn(vx);
2065 x = RgX_recip_shallow(x);
2066 if (p)
2067 { /* multiply by y[0]^p [beware dummy vars from FpX_FpXY_resultant] */
2068 GEN t = y_lead;
2069 if (T && typ(t) == t_POL && varn(t) == varn(T))
2070 t = RgXQ_powu(t, p, T);
2071 else
2072 t = gpowgs(t, p);
2073 for (i=2; i<lx; i++)
2074 {
2075 GEN c = gmul(gel(x,i), t);
2076 gel(x,i) = rem(c,T);
2077 }
2078 if (!T) return gerepileupto(av, x);
2079 }
2080 return gerepilecopy(av, x);
2081 }
2082
2083 GEN
2084 RgX_pseudorem(GEN x, GEN y) { return RgXQX_pseudorem(x,y, NULL); }
2085
2086 /* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
2087 GEN
2088 RgXQX_pseudodivrem(GEN x, GEN y, GEN T, GEN *ptr)
2089 {
2090 long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
2091 pari_sp av = avma, av2;
2092 GEN z, r, ypow, y_lead;
2093
2094 if (!signe(y)) pari_err_INV("RgXQX_pseudodivrem",y);
2095 dy = degpol(y); y_lead = gel(y,dy+2);
2096 if (gequal1(y_lead)) return T? RgXQX_divrem(x,y, T, ptr): RgX_divrem(x,y, ptr);
2097 dx = degpol(x);
2098 if (dx < dy) { *ptr = RgX_copy(x); return pol_0(vx); }
2099 if (dx == dy)
2100 {
2101 GEN x_lead = gel(x,lg(x)-1);
2102 x = RgX_renormalize_lg(leafcopy(x), lg(x)-1);
2103 y = RgX_renormalize_lg(leafcopy(y), lg(y)-1);
2104 r = RgX_sub(RgX_Rg_mul(x, y_lead), RgX_Rg_mul(y, x_lead));
2105 *ptr = gerepileupto(av, r); return scalarpol(x_lead, vx);
2106 }
2107 (void)new_chunk(2);
2108 x = RgX_recip_shallow(x)+2;
2109 y = RgX_recip_shallow(y)+2;
2110 /* pay attention to sparse divisors */
2111 for (i = 1; i <= dy; i++)
2112 if (isexactzero(gel(y,i))) gel(y,i) = NULL;
2113 dz = dx-dy; p = dz+1;
2114 lz = dz+3;
2115 z = cgetg(lz, t_POL);
2116 z[1] = evalsigne(1) | evalvarn(vx);
2117 for (i = 2; i < lz; i++) gel(z,i) = gen_0;
2118 ypow = new_chunk(dz+1);
2119 gel(ypow,0) = gen_1;
2120 gel(ypow,1) = y_lead;
2121 for (i=2; i<=dz; i++)
2122 {
2123 GEN c = gmul(gel(ypow,i-1), y_lead);
2124 gel(ypow,i) = rem(c,T);
2125 }
2126 av2 = avma;
2127 for (iz=2;;)
2128 {
2129 p--;
2130 gel(z,iz++) = rem(gmul(gel(x,0), gel(ypow,p)), T);
2131 for (i=1; i<=dy; i++)
2132 {
2133 GEN c = gmul(y_lead, gel(x,i));
2134 if (gel(y,i)) c = gsub(c, gmul(gel(x,0),gel(y,i)));
2135 gel(x,i) = rem(c, T);
2136 }
2137 for ( ; i<=dx; i++)
2138 {
2139 GEN c = gmul(y_lead, gel(x,i));
2140 gel(x,i) = rem(c,T);
2141 }
2142 x++; dx--;
2143 while (dx >= dy && gequal0(gel(x,0))) { x++; dx--; iz++; }
2144 if (dx < dy) break;
2145 if (gc_needed(av2,1))
2146 {
2147 GEN X = x-2;
2148 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_pseudodivrem dx=%ld >= %ld",dx,dy);
2149 X[0] = evaltyp(t_POL)|evallg(dx+3); X[1] = z[1]; /* hack */
2150 gerepileall(av2,2, &X, &z); x = X+2;
2151 }
2152 }
2153 while (dx >= 0 && gequal0(gel(x,0))) { x++; dx--; }
2154 if (dx < 0)
2155 x = pol_0(vx);
2156 else
2157 {
2158 lx = dx+3; x -= 2;
2159 x[0] = evaltyp(t_POL) | evallg(lx);
2160 x[1] = evalsigne(1) | evalvarn(vx);
2161 x = RgX_recip_shallow(x);
2162 }
2163 z = RgX_recip_shallow(z);
2164 r = x;
2165 if (p)
2166 {
2167 GEN c = gel(ypow,p); r = RgX_Rg_mul(r, c);
2168 if (T && typ(c) == t_POL && varn(c) == varn(T)) r = RgXQX_red(r, T);
2169 }
2170 gerepileall(av, 2, &z, &r);
2171 *ptr = r; return z;
2172 }
2173 GEN
2174 RgX_pseudodivrem(GEN x, GEN y, GEN *ptr)
2175 { return RgXQX_pseudodivrem(x,y,NULL,ptr); }
2176
2177 GEN
2178 RgXQX_mul(GEN x, GEN y, GEN T)
2179 {
2180 return RgXQX_red(RgX_mul(x,y), T);
2181 }
2182 GEN
2183 RgX_Rg_mul(GEN y, GEN x) {
2184 long i, ly;
2185 GEN z = cgetg_copy(y, &ly); z[1] = y[1];
2186 if (ly == 2) return z;
2187 for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
2188 return normalizepol_lg(z,ly);
2189 }
2190 GEN
2191 RgX_muls(GEN y, long x) {
2192 long i, ly;
2193 GEN z = cgetg_copy(y, &ly); z[1] = y[1];
2194 if (ly == 2) return z;
2195 for (i = 2; i < ly; i++) gel(z,i) = gmulsg(x,gel(y,i));
2196 return normalizepol_lg(z,ly);
2197 }
2198 GEN
2199 RgXQX_RgXQ_mul(GEN x, GEN y, GEN T)
2200 {
2201 return RgXQX_red(RgX_Rg_mul(x,y), T);
2202 }
2203 GEN
2204 RgXQV_RgXQ_mul(GEN v, GEN x, GEN T)
2205 {
2206 return RgXQV_red(RgV_Rg_mul(v,x), T);
2207 }
2208
2209 GEN
2210 RgXQX_sqr(GEN x, GEN T)
2211 {
2212 return RgXQX_red(RgX_sqr(x), T);
2213 }
2214
2215 GEN
2216 RgXQX_powers(GEN P, long n, GEN T)
2217 {
2218 GEN v = cgetg(n+2, t_VEC);
2219 long i;
2220 gel(v, 1) = pol_1(varn(T));
2221 if (n==0) return v;
2222 gel(v, 2) = gcopy(P);
2223 for (i = 2; i <= n; i++) gel(v,i+1) = RgXQX_mul(P, gel(v,i), T);
2224 return v;
2225 }
2226
2227 static GEN
2228 _add(void *data, GEN x, GEN y) { (void)data; return RgX_add(x, y); }
2229 static GEN
2230 _sub(void *data, GEN x, GEN y) { (void)data; return RgX_sub(x, y); }
2231 static GEN
2232 _sqr(void *data, GEN x) { return RgXQ_sqr(x, (GEN)data); }
2233 static GEN
2234 _mul(void *data, GEN x, GEN y) { return RgXQ_mul(x,y, (GEN)data); }
2235 static GEN
2236 _cmul(void *data, GEN P, long a, GEN x) { (void)data; return RgX_Rg_mul(x,gel(P,a+2)); }
2237 static GEN
2238 _one(void *data) { return pol_1(varn((GEN)data)); }
2239 static GEN
2240 _zero(void *data) { return pol_0(varn((GEN)data)); }
2241 static GEN
2242 _red(void *data, GEN x) { (void)data; return gcopy(x); }
2243
2244 static struct bb_algebra RgXQ_algebra = { _red, _add, _sub,
2245 _mul, _sqr, _one, _zero };
2246
2247 GEN
2248 RgX_RgXQV_eval(GEN Q, GEN x, GEN T)
2249 {
2250 return gen_bkeval_powers(Q,degpol(Q),x,(void*)T,&RgXQ_algebra,_cmul);
2251 }
2252
2253 GEN
2254 RgX_RgXQ_eval(GEN Q, GEN x, GEN T)
2255 {
2256 int use_sqr = 2*degpol(x) >= degpol(T);
2257 return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)T,&RgXQ_algebra,_cmul);
2258 }
2259
2260 /* mod X^n */
2261 struct modXn {
2262 long v; /* varn(X) */
2263 long n;
2264 } ;
2265 static GEN
2266 _sqrXn(void *data, GEN x) {
2267 struct modXn *S = (struct modXn*)data;
2268 return RgXn_sqr(x, S->n);
2269 }
2270 static GEN
2271 _mulXn(void *data, GEN x, GEN y) {
2272 struct modXn *S = (struct modXn*)data;
2273 return RgXn_mul(x,y, S->n);
2274 }
2275 static GEN
2276 _oneXn(void *data) {
2277 struct modXn *S = (struct modXn*)data;
2278 return pol_1(S->v);
2279 }
2280 static GEN
2281 _zeroXn(void *data) {
2282 struct modXn *S = (struct modXn*)data;
2283 return pol_0(S->v);
2284 }
2285 static struct bb_algebra RgXn_algebra = { _red, _add, _sub, _mulXn, _sqrXn,
2286 _oneXn, _zeroXn };
2287
2288 GEN
2289 RgXn_powers(GEN x, long m, long n)
2290 {
2291 long d = degpol(x);
2292 int use_sqr = (d<<1) >= n;
2293 struct modXn S;
2294 S.v = varn(x); S.n = n;
2295 return gen_powers(x,m,use_sqr,(void*)&S,_sqrXn,_mulXn,_oneXn);
2296 }
2297
2298 GEN
2299 RgXn_powu_i(GEN x, ulong m, long n)
2300 {
2301 struct modXn S;
2302 long v;
2303 if (n == 1) return x;
2304 v = RgX_valrem(x, &x);
2305 if (v) n -= m * v;
2306 S.v = varn(x); S.n = n;
2307 x = gen_powu_i(x, m, (void*)&S,_sqrXn,_mulXn);
2308 if (v) x = RgX_shift_shallow(x, m * v);
2309 return x;
2310 }
2311 GEN
2312 RgXn_powu(GEN x, ulong m, long n)
2313 {
2314 pari_sp av;
2315 if (n == 1) return gcopy(x);
2316 av = avma; return gerepilecopy(av, RgXn_powu_i(x, m, n));
2317 }
2318
2319 GEN
2320 RgX_RgXnV_eval(GEN Q, GEN x, long n)
2321 {
2322 struct modXn S;
2323 S.v = varn(gel(x,2)); S.n = n;
2324 return gen_bkeval_powers(Q,degpol(Q),x,(void*)&S,&RgXn_algebra,_cmul);
2325 }
2326
2327 GEN
2328 RgX_RgXn_eval(GEN Q, GEN x, long n)
2329 {
2330 int use_sqr = 2*degpol(x) >= n;
2331 struct modXn S;
2332 S.v = varn(x); S.n = n;
2333 return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2334 }
2335
2336 /* Q(x) mod t^n, x in R[t], n >= 1 */
2337 GEN
2338 RgXn_eval(GEN Q, GEN x, long n)
2339 {
2340 long d = degpol(x);
2341 int use_sqr;
2342 struct modXn S;
2343 if (d == 1 && isrationalzero(gel(x,2)))
2344 {
2345 GEN y = RgX_unscale(Q, gel(x,3));
2346 setvarn(y, varn(x)); return y;
2347 }
2348 S.v = varn(x);
2349 S.n = n;
2350 use_sqr = (d<<1) >= n;
2351 return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&S,&RgXn_algebra,_cmul);
2352 }
2353
2354 /* (f*g mod t^n) \ t^n2, assuming 2*n2 >= n */
2355 static GEN
2356 RgXn_mulhigh(GEN f, GEN g, long n2, long n)
2357 {
2358 GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2359 return RgX_add(RgX_mulhigh_i(fl, g, n2), RgXn_mul(fh, g, n - n2));
2360 }
2361
2362 /* (f^2 mod t^n) \ t^n2, assuming 2*n2 >= n */
2363 static GEN
2364 RgXn_sqrhigh(GEN f, long n2, long n)
2365 {
2366 GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2367 return RgX_add(RgX_mulhigh_i(fl, f, n2), RgXn_mul(fh, f, n - n2));
2368 }
2369
2370 GEN
2371 RgXn_inv_i(GEN f, long e)
2372 {
2373 pari_sp av;
2374 ulong mask;
2375 GEN W, a;
2376 long v = varn(f), n = 1;
2377
2378 if (!signe(f)) pari_err_INV("RgXn_inv",f);
2379 a = ginv(gel(f,2));
2380 if (e == 1) return scalarpol(a, v);
2381 else if (e == 2)
2382 {
2383 GEN b;
2384 if (degpol(f) <= 0 || gequal0(b = gel(f,3))) return scalarpol(a, v);
2385 b = gneg(b);
2386 if (!gequal1(a)) b = gmul(b, gsqr(a));
2387 return deg1pol(b, a, v);
2388 }
2389 av = avma;
2390 W = scalarpol_shallow(a,v);
2391 mask = quadratic_prec_mask(e);
2392 while (mask > 1)
2393 {
2394 GEN u, fr;
2395 long n2 = n;
2396 n<<=1; if (mask & 1) n--;
2397 mask >>= 1;
2398 fr = RgXn_red_shallow(f, n);
2399 u = RgXn_mul(W, RgXn_mulhigh(fr, W, n2, n), n-n2);
2400 W = RgX_sub(W, RgX_shift_shallow(u, n2));
2401 if (gc_needed(av,2))
2402 {
2403 if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_inv, e = %ld", n);
2404 W = gerepileupto(av, W);
2405 }
2406 }
2407 return W;
2408 }
2409
2410 static GEN
2411 RgXn_inv_FpX(GEN x, long e, GEN p)
2412 {
2413 pari_sp av = avma;
2414 GEN r;
2415 if (lgefint(p) == 3)
2416 {
2417 ulong pp = uel(p, 2);
2418 if (pp == 2)
2419 r = F2x_to_ZX(F2xn_inv(RgX_to_F2x(x), e));
2420 else
2421 r = Flx_to_ZX_inplace(Flxn_inv(RgX_to_Flx(x, pp), e, pp));
2422 }
2423 else
2424 r = FpXn_inv(RgX_to_FpX(x, p), e, p);
2425 return gerepileupto(av, FpX_to_mod(r, p));
2426 }
2427
2428 static GEN
2429 RgXn_inv_FpXQX(GEN x, long n, GEN pol, GEN p)
2430 {
2431 pari_sp av = avma;
2432 GEN r, T = RgX_to_FpX(pol, p);
2433 if (signe(T) == 0) pari_err_OP("/", gen_1, x);
2434 r = FpXQXn_inv(RgX_to_FpXQX(x, T, p), n, T, p);
2435 return gerepileupto(av, FpXQX_to_mod(r, T, p));
2436 }
2437
2438 #define code(t1,t2) ((t1 << 6) | t2)
2439
2440 static GEN
2441 RgXn_inv_fast(GEN x, long e)
2442 {
2443 GEN p, pol;
2444 long pa;
2445 long t = RgX_type(x,&p,&pol,&pa);
2446 switch(t)
2447 {
2448 case t_INTMOD: return RgXn_inv_FpX(x, e, p);
2449 case code(t_POLMOD, t_INTMOD):
2450 return RgXn_inv_FpXQX(x, e, pol, p);
2451 default: return NULL;
2452 }
2453 }
2454 #undef code
2455
2456 GEN
2457 RgXn_inv(GEN f, long e)
2458 {
2459 pari_sp av = avma;
2460 GEN h = RgXn_inv_fast(f, e);
2461 if (h) return h;
2462 return gerepileupto(av, RgXn_inv_i(f, e));
2463 }
2464
2465 /* Compute intformal(x^n*S)/x^(n+1) */
2466 static GEN
2467 RgX_integXn(GEN x, long n)
2468 {
2469 long i, lx = lg(x);
2470 GEN y;
2471 if (lx == 2) return RgX_copy(x);
2472 y = cgetg(lx, t_POL); y[1] = x[1];
2473 for (i=2; i<lx; i++)
2474 gel(y,i) = gdivgs(gel(x,i), n+i-1);
2475 return RgX_renormalize_lg(y, lx);;
2476 }
2477
2478 GEN
2479 RgXn_expint(GEN h, long e)
2480 {
2481 pari_sp av = avma, av2;
2482 long v = varn(h), n;
2483 GEN f = pol_1(v), g;
2484 ulong mask;
2485
2486 if (!signe(h)) return f;
2487 g = pol_1(v);
2488 n = 1; mask = quadratic_prec_mask(e);
2489 av2 = avma;
2490 for (;mask>1;)
2491 {
2492 GEN u, w;
2493 long n2 = n;
2494 n<<=1; if (mask & 1) n--;
2495 mask >>= 1;
2496 u = RgXn_mul(g, RgX_mulhigh_i(f, RgXn_red_shallow(h, n2-1), n2-1), n-n2);
2497 u = RgX_add(u, RgX_shift_shallow(RgXn_red_shallow(h, n-1), 1-n2));
2498 w = RgXn_mul(f, RgX_integXn(u, n2-1), n-n2);
2499 f = RgX_add(f, RgX_shift_shallow(w, n2));
2500 if (mask<=1) break;
2501 u = RgXn_mul(g, RgXn_mulhigh(f, g, n2, n), n-n2);
2502 g = RgX_sub(g, RgX_shift_shallow(u, n2));
2503 if (gc_needed(av2,2))
2504 {
2505 if (DEBUGMEM>1) pari_warn(warnmem,"RgXn_expint, e = %ld", n);
2506 gerepileall(av2, 2, &f, &g);
2507 }
2508 }
2509 return gerepileupto(av, f);
2510 }
2511
2512 GEN
2513 RgXn_exp(GEN h, long e)
2514 {
2515 long d = degpol(h);
2516 if (d < 0) return pol_1(varn(h));
2517 if (!d || !gequal0(gel(h,2)))
2518 pari_err_DOMAIN("RgXn_exp","valuation", "<", gen_1, h);
2519 return RgXn_expint(RgX_deriv(h), e);
2520 }
2521
2522 GEN
2523 RgXn_reverse(GEN f, long e)
2524 {
2525 pari_sp av = avma, av2;
2526 ulong mask;
2527 GEN fi, a, df, W, an;
2528 long v = varn(f), n=1;
2529 if (degpol(f)<1 || !gequal0(gel(f,2)))
2530 pari_err_INV("serreverse",f);
2531 fi = ginv(gel(f,3));
2532 a = deg1pol_shallow(fi,gen_0,v);
2533 if (e <= 2) return gerepilecopy(av, a);
2534 W = scalarpol(fi,v);
2535 df = RgX_deriv(f);
2536 mask = quadratic_prec_mask(e);
2537 av2 = avma;
2538 for (;mask>1;)
2539 {
2540 GEN u, fa, fr;
2541 long n2 = n, rt;
2542 n<<=1; if (mask & 1) n--;
2543 mask >>= 1;
2544 fr = RgXn_red_shallow(f, n);
2545 rt = brent_kung_optpow(degpol(fr), 4, 3);
2546 an = RgXn_powers(a, rt, n);
2547 if (n>1)
2548 {
2549 long n4 = (n2+1)>>1;
2550 GEN dfr = RgXn_red_shallow(df, n2);
2551 dfr = RgX_RgXnV_eval(dfr, RgXnV_red_shallow(an, n2), n2);
2552 u = RgX_shift(RgX_Rg_sub(RgXn_mul(W, dfr, n2), gen_1), -n4);
2553 W = RgX_sub(W, RgX_shift(RgXn_mul(u, W, n2-n4), n4));
2554 }
2555 fa = RgX_sub(RgX_RgXnV_eval(fr, an, n), pol_x(v));
2556 fa = RgX_shift(fa, -n2);
2557 a = RgX_sub(a, RgX_shift(RgXn_mul(W, fa, n-n2), n2));
2558 if (gc_needed(av2,2))
2559 {
2560 if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_reverse, e = %ld", n);
2561 gerepileall(av2, 2, &a, &W);
2562 }
2563 }
2564 return gerepileupto(av, a);
2565 }
2566
2567 GEN
2568 RgXn_sqrt(GEN h, long e)
2569 {
2570 pari_sp av = avma, av2;
2571 long v = varn(h), n = 1;
2572 GEN f = scalarpol(gen_1, v), df = f;
2573 ulong mask = quadratic_prec_mask(e);
2574 if (degpol(h)<0 || !gequal1(gel(h,2)))
2575 pari_err_SQRTN("RgXn_sqrt",h);
2576 av2 = avma;
2577 while(1)
2578 {
2579 long n2 = n, m;
2580 GEN g;
2581 n<<=1; if (mask & 1) n--;
2582 mask >>= 1;
2583 m = n-n2;
2584 g = RgX_sub(RgXn_sqrhigh(f, n2, n), RgX_shift_shallow(RgXn_red_shallow(h, n),-n2));
2585 f = RgX_sub(f, RgX_shift_shallow(RgXn_mul(gmul2n(df, -1), g, m), n2));
2586 if (mask==1) return gerepileupto(av, f);
2587 g = RgXn_mul(df, RgXn_mulhigh(df, f, n2, n), m);
2588 df = RgX_sub(df, RgX_shift_shallow(g, n2));
2589 if (gc_needed(av2,2))
2590 {
2591 if(DEBUGMEM>1) pari_warn(warnmem,"RgXn_sqrt, e = %ld", n);
2592 gerepileall(av2, 2, &f, &df);
2593 }
2594 }
2595 }
2596
2597 /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2598 GEN
2599 RgXQ_powu(GEN x, ulong n, GEN T)
2600 {
2601 pari_sp av = avma;
2602
2603 if (!n) return pol_1(varn(x));
2604 if (n == 1) return RgX_copy(x);
2605 x = gen_powu_i(x, n, (void*)T, &_sqr, &_mul);
2606 return gerepilecopy(av, x);
2607 }
2608 /* x,T in Rg[X], n in N, compute lift(x^n mod T)) */
2609 GEN
2610 RgXQ_pow(GEN x, GEN n, GEN T)
2611 {
2612 pari_sp av;
2613 long s = signe(n);
2614
2615 if (!s) return pol_1(varn(x));
2616 if (is_pm1(n) == 1)
2617 return (s < 0)? RgXQ_inv(x, T): RgX_copy(x);
2618 av = avma;
2619 if (s < 0) x = RgXQ_inv(x, T);
2620 x = gen_pow_i(x, n, (void*)T, &_sqr, &_mul);
2621 return gerepilecopy(av, x);
2622 }
2623 static GEN
2624 _ZXQsqr(void *data, GEN x) { return ZXQ_sqr(x, (GEN)data); }
2625 static GEN
2626 _ZXQmul(void *data, GEN x, GEN y) { return ZXQ_mul(x,y, (GEN)data); }
2627
2628 /* generates the list of powers of x of degree 0,1,2,...,l*/
2629 GEN
2630 ZXQ_powers(GEN x, long l, GEN T)
2631 {
2632 int use_sqr = 2*degpol(x) >= degpol(T);
2633 return gen_powers(x, l, use_sqr, (void *)T,_ZXQsqr,_ZXQmul,_one);
2634 }
2635
2636 /* x,T in Z[X], n in N, compute lift(x^n mod T)) */
2637 GEN
2638 ZXQ_powu(GEN x, ulong n, GEN T)
2639 {
2640 pari_sp av = avma;
2641
2642 if (!n) return pol_1(varn(x));
2643 if (n == 1) return ZX_copy(x);
2644 x = gen_powu_i(x, n, (void*)T, &_ZXQsqr, &_ZXQmul);
2645 return gerepilecopy(av, x);
2646 }
2647
2648 /* generates the list of powers of x of degree 0,1,2,...,l*/
2649 GEN
2650 RgXQ_powers(GEN x, long l, GEN T)
2651 {
2652 int use_sqr = 2*degpol(x) >= degpol(T);
2653 return gen_powers(x, l, use_sqr, (void *)T,_sqr,_mul,_one);
2654 }
2655
2656 /* a in K = Q[X]/(T), returns [a^0, ..., a^n] */
2657 GEN
2658 QXQ_powers(GEN a, long n, GEN T)
2659 {
2660 GEN den, v = RgXQ_powers(Q_remove_denom(a, &den), n, T);
2661 /* den*a integral; v[i+1] = (den*a)^i in K */
2662 if (den)
2663 { /* restore denominators */
2664 GEN d = den;
2665 long i;
2666 gel(v,2) = a;
2667 for (i=3; i<=n+1; i++) {
2668 d = mulii(d,den);
2669 gel(v,i) = RgX_Rg_div(gel(v,i), d);
2670 }
2671 }
2672 return v;
2673 }
2674
2675 static GEN
2676 do_QXQ_eval(GEN v, long imin, GEN a, GEN T)
2677 {
2678 long l, i, m = 0;
2679 GEN dz, z;
2680 GEN V = cgetg_copy(v, &l);
2681 for (i = imin; i < l; i++)
2682 {
2683 GEN c = gel(v, i);
2684 if (typ(c) == t_POL) m = maxss(m, degpol(c));
2685 }
2686 z = Q_remove_denom(QXQ_powers(a, m, T), &dz);
2687 for (i = 1; i < imin; i++) V[i] = v[i];
2688 for (i = imin; i < l; i++)
2689 {
2690 GEN c = gel(v,i);
2691 if (typ(c) == t_POL) c = QX_ZXQV_eval(c, z, dz);
2692 gel(V,i) = c;
2693 }
2694 return V;
2695 }
2696 /* [ s(a mod T) | s <- lift(v) ], a,T are QX, v a QXV */
2697 GEN
2698 QXV_QXQ_eval(GEN v, GEN a, GEN T)
2699 { return do_QXQ_eval(v, 1, a, T); }
2700 GEN
2701 QXX_QXQ_eval(GEN v, GEN a, GEN T)
2702 { return normalizepol(do_QXQ_eval(v, 2, a, T)); }
2703
2704 GEN
2705 RgXQ_matrix_pow(GEN y, long n, long m, GEN P)
2706 {
2707 return RgXV_to_RgM(RgXQ_powers(y,m-1,P),n);
2708 }
2709
2710 GEN
2711 RgXQ_norm(GEN x, GEN T)
2712 {
2713 pari_sp av;
2714 long dx = degpol(x);
2715 GEN L, y;
2716
2717 av = avma; y = resultant(T, x);
2718 L = leading_coeff(T);
2719 if (gequal1(L) || !signe(x)) return y;
2720 return gerepileupto(av, gdiv(y, gpowgs(L, dx)));
2721 }
2722
2723 GEN
2724 RgX_blocks(GEN P, long n, long m)
2725 {
2726 GEN z = cgetg(m+1,t_VEC);
2727 long i,j, k=2, l = lg(P);
2728 for(i=1; i<=m; i++)
2729 {
2730 GEN zi = cgetg(n+2,t_POL);
2731 zi[1] = P[1];
2732 gel(z,i) = zi;
2733 for(j=2; j<n+2; j++)
2734 gel(zi, j) = k==l ? gen_0 : gel(P,k++);
2735 zi = RgX_renormalize_lg(zi, n+2);
2736 }
2737 return z;
2738 }
2739
2740 /* write p(X) = e(X^2) + Xo(X^2), shallow function */
2741 void
2742 RgX_even_odd(GEN p, GEN *pe, GEN *po)
2743 {
2744 long n = degpol(p), v = varn(p), n0, n1, i;
2745 GEN p0, p1;
2746
2747 if (n <= 0) { *pe = RgX_copy(p); *po = zeropol(v); return; }
2748
2749 n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
2750 p0 = cgetg(n0+2, t_POL); p0[1] = evalvarn(v)|evalsigne(1);
2751 p1 = cgetg(n1+2, t_POL); p1[1] = evalvarn(v)|evalsigne(1);
2752 for (i=0; i<n1; i++)
2753 {
2754 p0[2+i] = p[2+(i<<1)];
2755 p1[2+i] = p[3+(i<<1)];
2756 }
2757 if (n1 != n0)
2758 p0[2+i] = p[2+(i<<1)];
2759 *pe = normalizepol(p0);
2760 *po = normalizepol(p1);
2761 }
2762
2763 /* write p(X) = a_0(X^k) + Xa_1(X^k) + ... + X^(k-1)a_{k-1}(X^k), shallow function */
2764 GEN
2765 RgX_splitting(GEN p, long k)
2766 {
2767 long n = degpol(p), v = varn(p), m, i, j, l;
2768 GEN r;
2769
2770 m = n/k;
2771 r = cgetg(k+1,t_VEC);
2772 for(i=1; i<=k; i++)
2773 {
2774 gel(r,i) = cgetg(m+3, t_POL);
2775 mael(r,i,1) = evalvarn(v)|evalsigne(1);
2776 }
2777 for (j=1, i=0, l=2; i<=n; i++)
2778 {
2779 gmael(r,j,l) = gel(p,2+i);
2780 if (j==k) { j=1; l++; } else j++;
2781 }
2782 for(i=1; i<=k; i++)
2783 gel(r,i) = normalizepol_lg(gel(r,i),i<j?l+1:l);
2784 return r;
2785 }
2786
2787 /*******************************************************************/
2788 /* */
2789 /* Kronecker form */
2790 /* */
2791 /*******************************************************************/
2792
2793 /* z in R[Y] representing an elt in R[X,Y] mod T(Y) in Kronecker form,
2794 * i.e subst(lift(z), x, y^(2deg(z)-1)). Recover the "real" z, with
2795 * normalized coefficients */
2796 GEN
2797 Kronecker_to_mod(GEN z, GEN T)
2798 {
2799 long i,j,lx,l = lg(z), N = (degpol(T)<<1) + 1;
2800 GEN x, t = cgetg(N,t_POL);
2801 t[1] = T[1];
2802 lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL);
2803 x[1] = z[1];
2804 T = RgX_copy(T);
2805 for (i=2; i<lx+2; i++, z+= N-2)
2806 {
2807 for (j=2; j<N; j++) gel(t,j) = gel(z,j);
2808 gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2809 }
2810 N = (l-2) % (N-2) + 2;
2811 for (j=2; j<N; j++) t[j] = z[j];
2812 gel(x,i) = mkpolmod(RgX_rem(normalizepol_lg(t,N), T), T);
2813 return normalizepol_lg(x, i+1);
2814 }
2815
2816 /*******************************************************************/
2817 /* */
2818 /* Domain detection */
2819 /* */
2820 /*******************************************************************/
2821
2822 static GEN
2823 zero_FpX_mod(GEN p, long v)
2824 {
2825 GEN r = cgetg(3,t_POL);
2826 r[1] = evalvarn(v);
2827 gel(r,2) = mkintmod(gen_0, icopy(p));
2828 return r;
2829 }
2830
2831 static GEN
2832 RgX_mul_FpX(GEN x, GEN y, GEN p)
2833 {
2834 pari_sp av = avma;
2835 GEN r;
2836 if (lgefint(p) == 3)
2837 {
2838 ulong pp = uel(p, 2);
2839 r = Flx_to_ZX_inplace(Flx_mul(RgX_to_Flx(x, pp),
2840 RgX_to_Flx(y, pp), pp));
2841 }
2842 else
2843 r = FpX_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2844 if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2845 return gerepileupto(av, FpX_to_mod(r, p));
2846 }
2847
2848 static GEN
2849 zero_FpXQX_mod(GEN pol, GEN p, long v)
2850 {
2851 GEN r = cgetg(3,t_POL);
2852 r[1] = evalvarn(v);
2853 gel(r,2) = mkpolmod(mkintmod(gen_0, icopy(p)), gcopy(pol));
2854 return r;
2855 }
2856
2857 static GEN
2858 RgX_mul_FpXQX(GEN x, GEN y, GEN pol, GEN p)
2859 {
2860 pari_sp av = avma;
2861 long dT;
2862 GEN kx, ky, r;
2863 GEN T = RgX_to_FpX(pol, p);
2864 if (signe(T)==0) pari_err_OP("*", x, y);
2865 dT = degpol(T);
2866 kx = ZXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
2867 ky = ZXX_to_Kronecker(RgX_to_FpXQX(y, T, p), dT);
2868 r = FpX_mul(kx, ky, p);
2869 if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
2870 return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
2871 }
2872
2873 static GEN
2874 RgX_liftred(GEN x, GEN T)
2875 { return RgXQX_red(liftpol_shallow(x), T); }
2876
2877 static GEN
2878 RgX_mul_QXQX(GEN x, GEN y, GEN T)
2879 {
2880 pari_sp av = avma;
2881 long dT = degpol(T);
2882 GEN r = QX_mul(ZXX_to_Kronecker(RgX_liftred(x, T), dT),
2883 ZXX_to_Kronecker(RgX_liftred(y, T), dT));
2884 return gerepileupto(av, Kronecker_to_mod(r, T));
2885 }
2886
2887 static GEN
2888 RgX_sqr_FpX(GEN x, GEN p)
2889 {
2890 pari_sp av = avma;
2891 GEN r;
2892 if (lgefint(p) == 3)
2893 {
2894 ulong pp = uel(p, 2);
2895 r = Flx_to_ZX_inplace(Flx_sqr(RgX_to_Flx(x, pp), pp));
2896 }
2897 else
2898 r = FpX_sqr(RgX_to_FpX(x, p), p);
2899 if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2900 return gerepileupto(av, FpX_to_mod(r, p));
2901 }
2902
2903 static GEN
2904 RgX_sqr_FpXQX(GEN x, GEN pol, GEN p)
2905 {
2906 pari_sp av = avma;
2907 long dT;
2908 GEN kx, r, T = RgX_to_FpX(pol, p);
2909 if (signe(T)==0) pari_err_OP("*",x,x);
2910 dT = degpol(T);
2911 kx = ZXX_to_Kronecker(RgX_to_FpXQX(x, T, p), dT);
2912 r = FpX_sqr(kx, p);
2913 if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
2914 return gerepileupto(av, Kronecker_to_mod(FpX_to_mod(r, p), pol));
2915 }
2916
2917 static GEN
2918 RgX_sqr_QXQX(GEN x, GEN T)
2919 {
2920 pari_sp av = avma;
2921 long dT = degpol(T);
2922 GEN r = QX_sqr(ZXX_to_Kronecker(RgX_liftred(x, T), dT));
2923 return gerepileupto(av, Kronecker_to_mod(r, T));
2924 }
2925
2926 static GEN
2927 RgX_rem_FpX(GEN x, GEN y, GEN p)
2928 {
2929 pari_sp av = avma;
2930 GEN r;
2931 if (lgefint(p) == 3)
2932 {
2933 ulong pp = uel(p, 2);
2934 r = Flx_to_ZX_inplace(Flx_rem(RgX_to_Flx(x, pp),
2935 RgX_to_Flx(y, pp), pp));
2936 }
2937 else
2938 r = FpX_rem(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
2939 if (signe(r)==0) { set_avma(av); return zero_FpX_mod(p, varn(x)); }
2940 return gerepileupto(av, FpX_to_mod(r, p));
2941 }
2942
2943 static GEN
2944 RgX_rem_QXQX(GEN x, GEN y, GEN T)
2945 {
2946 pari_sp av = avma;
2947 GEN r;
2948 r = RgXQX_rem(RgX_liftred(x, T), RgX_liftred(y, T), T);
2949 return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
2950 }
2951 static GEN
2952 RgX_rem_FpXQX(GEN x, GEN y, GEN pol, GEN p)
2953 {
2954 pari_sp av = avma;
2955 GEN r;
2956 GEN T = RgX_to_FpX(pol, p);
2957 if (signe(T) == 0) pari_err_OP("%", x, y);
2958 if (lgefint(p) == 3)
2959 {
2960 ulong pp = uel(p, 2);
2961 GEN Tp = ZX_to_Flx(T, pp);
2962 r = FlxX_to_ZXX(FlxqX_rem(RgX_to_FlxqX(x, Tp, pp),
2963 RgX_to_FlxqX(y, Tp, pp), Tp, pp));
2964 }
2965 else
2966 r = FpXQX_rem(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
2967 if (signe(r)==0) { set_avma(av); return zero_FpXQX_mod(pol, p, varn(x)); }
2968 return gerepileupto(av, FpXQX_to_mod(r, T, p));
2969 }
2970
2971 #define code(t1,t2) ((t1 << 6) | t2)
2972 static GEN
2973 RgX_mul_fast(GEN x, GEN y)
2974 {
2975 GEN p, pol;
2976 long pa;
2977 long t = RgX_type2(x,y, &p,&pol,&pa);
2978 switch(t)
2979 {
2980 case t_INT: return ZX_mul(x,y);
2981 case t_FRAC: return QX_mul(x,y);
2982 case t_FFELT: return FFX_mul(x, y, pol);
2983 case t_INTMOD: return RgX_mul_FpX(x, y, p);
2984 case code(t_POLMOD, t_INT):
2985 case code(t_POLMOD, t_FRAC):
2986 return RgX_mul_QXQX(x, y, pol);
2987 case code(t_POLMOD, t_INTMOD):
2988 return RgX_mul_FpXQX(x, y, pol, p);
2989 default: return NULL;
2990 }
2991 }
2992 static GEN
2993 RgX_sqr_fast(GEN x)
2994 {
2995 GEN p, pol;
2996 long pa;
2997 long t = RgX_type(x,&p,&pol,&pa);
2998 switch(t)
2999 {
3000 case t_INT: return ZX_sqr(x);
3001 case t_FRAC: return QX_sqr(x);
3002 case t_FFELT: return FFX_sqr(x, pol);
3003 case t_INTMOD: return RgX_sqr_FpX(x, p);
3004 case code(t_POLMOD, t_INT):
3005 case code(t_POLMOD, t_FRAC):
3006 return RgX_sqr_QXQX(x, pol);
3007 case code(t_POLMOD, t_INTMOD):
3008 return RgX_sqr_FpXQX(x, pol, p);
3009 default: return NULL;
3010 }
3011 }
3012
3013 static GEN
3014 RgX_rem_fast(GEN x, GEN y)
3015 {
3016 GEN p, pol;
3017 long pa;
3018 long t = RgX_type2(x,y, &p,&pol,&pa);
3019 switch(t)
3020 {
3021 case t_INT: return ZX_is_monic(y) ? ZX_rem(x,y): NULL;
3022 case t_FRAC: return RgX_is_ZX(y) && ZX_is_monic(y) ? QX_ZX_rem(x,y): NULL;
3023 case t_FFELT: return FFX_rem(x, y, pol);
3024 case t_INTMOD: return RgX_rem_FpX(x, y, p);
3025 case code(t_POLMOD, t_INT):
3026 case code(t_POLMOD, t_FRAC):
3027 return RgX_rem_QXQX(x, y, pol);
3028 case code(t_POLMOD, t_INTMOD):
3029 return RgX_rem_FpXQX(x, y, pol, p);
3030 default: return NULL;
3031 }
3032 }
3033
3034 #undef code
3035
3036 GEN
3037 RgX_mul(GEN x, GEN y)
3038 {
3039 GEN z = RgX_mul_fast(x,y);
3040 if (!z) z = RgX_mul_i(x,y);
3041 return z;
3042 }
3043
3044 GEN
3045 RgX_sqr(GEN x)
3046 {
3047 GEN z = RgX_sqr_fast(x);
3048 if (!z) z = RgX_sqr_i(x);
3049 return z;
3050 }
3051
3052 GEN
3053 RgX_rem(GEN x, GEN y)
3054 {
3055 GEN z = RgX_rem_fast(x, y);
3056 if (!z) z = RgX_divrem_i(x, y, ONLY_REM);
3057 return z;
3058 }
3059
3060 GEN
3061 RgXn_mul(GEN f, GEN g, long n)
3062 {
3063 pari_sp av = avma;
3064 GEN h = RgX_mul_fast(f,g);
3065 if (!h) return RgXn_mul2(f,g,n);
3066 if (degpol(h) < n) return h;
3067 return gerepilecopy(av, RgXn_red_shallow(h, n));
3068 }
3069
3070 GEN
3071 RgXn_sqr(GEN f, long n)
3072 {
3073 pari_sp av = avma;
3074 GEN g = RgX_sqr_fast(f);
3075 if (!g) return RgXn_sqr2(f,n);
3076 if (degpol(g) < n) return g;
3077 return gerepilecopy(av, RgXn_red_shallow(g, n));
3078 }
3079
3080 /* (f*g) \/ x^n */
3081 GEN
3082 RgX_mulhigh_i(GEN f, GEN g, long n)
3083 {
3084 GEN h = RgX_mul_fast(f,g);
3085 return h? RgX_shift_shallow(h, -n): RgX_mulhigh_i2(f,g,n);
3086 }
3087
3088 /* (f*g) \/ x^n */
3089 GEN
3090 RgX_sqrhigh_i(GEN f, long n)
3091 {
3092 GEN h = RgX_sqr_fast(f);
3093 return h? RgX_shift_shallow(h, -n): RgX_sqrhigh_i2(f,n);
3094 }