Line data Source code
1 : /* Copyright (C) 2007 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 : /* ZX */
21 : /* */
22 : /*******************************************************************/
23 : void
24 369818 : RgX_check_QX(GEN x, const char *s)
25 369818 : { if (!RgX_is_QX(x)) pari_err_TYPE(stack_strcat(s," [not in Q[X]]"), x); }
26 : void
27 2979291 : RgX_check_ZX(GEN x, const char *s)
28 2979291 : { if (!RgX_is_ZX(x)) pari_err_TYPE(stack_strcat(s," [not in Z[X]]"), x); }
29 : long
30 109375 : ZX_max_lg(GEN x)
31 : {
32 109375 : long i, l = 0, lx = lg(x);
33 657614 : for (i = 2; i < lx; i++) l = maxss(l, lgefint(gel(x,i)));
34 109375 : return l;
35 : }
36 :
37 : GEN
38 39052817 : ZX_add(GEN x, GEN y)
39 : {
40 : long lx,ly,i;
41 : GEN z;
42 39052817 : lx = lg(x); ly = lg(y); if (lx < ly) swapspec(x,y, lx,ly);
43 39052817 : z = cgetg(lx,t_POL); z[1] = x[1];
44 232189403 : for (i=2; i<ly; i++) gel(z,i) = addii(gel(x,i),gel(y,i));
45 62240275 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
46 38955266 : if (lx == ly) z = ZX_renormalize(z, lx);
47 38954048 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
48 37740167 : return z;
49 : }
50 :
51 : GEN
52 29546750 : ZX_sub(GEN x,GEN y)
53 : {
54 29546750 : long i, lx = lg(x), ly = lg(y);
55 : GEN z;
56 29546750 : if (lx >= ly)
57 : {
58 29097030 : z = cgetg(lx,t_POL); z[1] = x[1];
59 118216576 : for (i=2; i<ly; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
60 29051949 : if (lx == ly)
61 : {
62 23410743 : z = ZX_renormalize(z, lx);
63 23409908 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); z = pol_0(varn(x)); }
64 : }
65 : else
66 17411479 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
67 : }
68 : else
69 : {
70 449720 : z = cgetg(ly,t_POL); z[1] = y[1];
71 1178121 : for (i=2; i<lx; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
72 1396466 : for ( ; i<ly; i++) gel(z,i) = negi(gel(y,i));
73 : }
74 29542326 : return z;
75 : }
76 :
77 : GEN
78 599389 : ZX_neg(GEN x)
79 : {
80 599389 : long i, l = lg(x);
81 599389 : GEN y = cgetg(l,t_POL);
82 2470814 : y[1] = x[1]; for(i=2; i<l; i++) gel(y,i) = negi(gel(x,i));
83 599432 : return y;
84 : }
85 : GEN
86 5940126 : ZX_copy(GEN x)
87 : {
88 5940126 : long i, l = lg(x);
89 5940126 : GEN y = cgetg(l, t_POL);
90 5939705 : y[1] = x[1];
91 23643652 : for (i=2; i<l; i++)
92 : {
93 17703550 : GEN c = gel(x,i);
94 17703550 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
95 : }
96 5940102 : return y;
97 : }
98 :
99 : GEN
100 55393 : scalar_ZX(GEN x, long v)
101 : {
102 : GEN z;
103 55393 : if (!signe(x)) return pol_0(v);
104 4488 : z = cgetg(3, t_POL);
105 4488 : z[1] = evalsigne(1) | evalvarn(v);
106 4488 : gel(z,2) = icopy(x); return z;
107 : }
108 :
109 : GEN
110 22897 : scalar_ZX_shallow(GEN x, long v)
111 : {
112 : GEN z;
113 22897 : if (!signe(x)) return pol_0(v);
114 22585 : z = cgetg(3, t_POL);
115 22585 : z[1] = evalsigne(1) | evalvarn(v);
116 22585 : gel(z,2) = x; return z;
117 : }
118 :
119 : GEN
120 87705 : ZX_Z_add(GEN y, GEN x)
121 : {
122 : long lz, i;
123 87705 : GEN z = cgetg_copy(y, &lz);
124 87707 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
125 76969 : z[1] = y[1];
126 76969 : gel(z,2) = addii(gel(y,2),x);
127 1831754 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
128 77086 : if (lz==3) z = ZX_renormalize(z,lz);
129 76969 : return z;
130 : }
131 : GEN
132 95030 : ZX_Z_add_shallow(GEN y, GEN x)
133 : {
134 : long lz, i;
135 95030 : GEN z = cgetg_copy(y, &lz);
136 95030 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX_shallow(x,varn(y)); }
137 94743 : z[1] = y[1];
138 94743 : gel(z,2) = addii(gel(y,2),x);
139 556891 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
140 94735 : if (lz==3) z = ZX_renormalize(z,lz);
141 94736 : return z;
142 : }
143 :
144 : GEN
145 161615 : ZX_Z_sub(GEN y, GEN x)
146 : {
147 : long lz, i;
148 161615 : GEN z = cgetg_copy(y, &lz);
149 161615 : if (lz == 2)
150 : { /* scalarpol(negi(x), v) */
151 441 : long v = varn(y);
152 441 : set_avma((pari_sp)(z + 2));
153 441 : if (!signe(x)) return pol_0(v);
154 441 : z = cgetg(3,t_POL);
155 441 : z[1] = evalvarn(v) | evalsigne(1);
156 441 : gel(z,2) = negi(x); return z;
157 : }
158 161174 : z[1] = y[1];
159 161174 : gel(z,2) = subii(gel(y,2),x);
160 787330 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
161 161174 : if (lz==3) z = ZX_renormalize(z,lz);
162 161174 : return z;
163 : }
164 :
165 : GEN
166 2739286 : Z_ZX_sub(GEN x, GEN y)
167 : {
168 : long lz, i;
169 2739286 : GEN z = cgetg_copy(y, &lz);
170 2739298 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
171 2739298 : z[1] = y[1];
172 2739298 : gel(z,2) = subii(x, gel(y,2));
173 13806141 : for(i=3; i<lz; i++) gel(z,i) = negi(gel(y,i));
174 2739638 : if (lz==3) z = ZX_renormalize(z,lz);
175 2739380 : return z;
176 : }
177 :
178 : GEN
179 5144811 : ZX_Z_divexact(GEN y,GEN x)
180 : {
181 5144811 : long i, l = lg(y);
182 5144811 : GEN z = cgetg(l,t_POL); z[1] = y[1];
183 32900959 : for(i=2; i<l; i++) gel(z,i) = diviiexact(gel(y,i),x);
184 5142674 : return z;
185 : }
186 :
187 : GEN
188 977809 : ZX_divuexact(GEN y, ulong x)
189 : {
190 977809 : long i, l = lg(y);
191 977809 : GEN z = cgetg(l,t_POL); z[1] = y[1];
192 9684745 : for(i=2; i<l; i++) gel(z,i) = diviuexact(gel(y,i),x);
193 977809 : return z;
194 : }
195 :
196 : GEN
197 112 : zx_z_divexact(GEN y, long x)
198 : {
199 112 : long i, l = lg(y);
200 112 : GEN z = cgetg(l,t_VECSMALL); z[1] = y[1];
201 595 : for (i=2; i<l; i++) z[i] = y[i]/x;
202 112 : return z;
203 : }
204 :
205 : GEN
206 15027438 : ZX_Z_mul(GEN y,GEN x)
207 : {
208 : GEN z;
209 : long i, l;
210 15027438 : if (!signe(x)) return pol_0(varn(y));
211 12394216 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
212 73306034 : for(i=2; i<l; i++) gel(z,i) = mulii(gel(y,i),x);
213 12391280 : return z;
214 : }
215 :
216 : GEN
217 17045 : ZX_mulu(GEN y, ulong x)
218 : {
219 : GEN z;
220 : long i, l;
221 17045 : if (!x) return pol_0(varn(y));
222 17045 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
223 146797 : for(i=2; i<l; i++) gel(z,i) = mului(x,gel(y,i));
224 17045 : return z;
225 : }
226 :
227 : GEN
228 6524801 : ZX_shifti(GEN y, long n)
229 : {
230 : GEN z;
231 : long i, l;
232 6524801 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
233 31338365 : for(i=2; i<l; i++) gel(z,i) = shifti(gel(y,i),n);
234 6520622 : return ZX_renormalize(z,l);
235 : }
236 :
237 : GEN
238 106764 : ZX_remi2n(GEN y, long n)
239 : {
240 : GEN z;
241 : long i, l;
242 106764 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
243 4904895 : for(i=2; i<l; i++) gel(z,i) = remi2n(gel(y,i),n);
244 104906 : return ZX_renormalize(z,l);
245 : }
246 :
247 : GEN
248 46052 : ZXT_remi2n(GEN z, long n)
249 : {
250 46052 : if (typ(z) == t_POL)
251 39775 : return ZX_remi2n(z, n);
252 : else
253 : {
254 6277 : long i,l = lg(z);
255 6277 : GEN x = cgetg(l, t_VEC);
256 18855 : for (i=1; i<l; i++) gel(x,i) = ZXT_remi2n(gel(z,i), n);
257 6281 : return x;
258 : }
259 : }
260 :
261 : GEN
262 20923 : zx_to_ZX(GEN z)
263 : {
264 20923 : long i, l = lg(z);
265 20923 : GEN x = cgetg(l,t_POL);
266 218512 : for (i=2; i<l; i++) gel(x,i) = stoi(z[i]);
267 20923 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
268 : }
269 :
270 : GEN
271 7944733 : ZX_deriv(GEN x)
272 : {
273 7944733 : long i,lx = lg(x)-1;
274 : GEN y;
275 :
276 7944733 : if (lx<3) return pol_0(varn(x));
277 7855758 : y = cgetg(lx,t_POL);
278 40599962 : for (i=2; i<lx ; i++) gel(y,i) = mului(i-1,gel(x,i+1));
279 7854607 : y[1] = x[1]; return y;
280 : }
281 :
282 : int
283 2582270 : ZX_equal(GEN V, GEN W)
284 : {
285 2582270 : long i, l = lg(V);
286 2582270 : if (lg(W) != l) return 0;
287 6511546 : for (i = 2; i < l; i++)
288 5182122 : if (!equalii(gel(V,i), gel(W,i))) return 0;
289 1329424 : return 1;
290 : }
291 :
292 : static long
293 336861824 : ZX_valspec(GEN x, long nx)
294 : {
295 : long vx;
296 504436392 : for (vx = 0; vx<nx ; vx++)
297 504457811 : if (signe(gel(x,vx))) break;
298 336861824 : return vx;
299 : }
300 :
301 : long
302 346893 : ZX_val(GEN x)
303 : {
304 : long vx;
305 346893 : if (!signe(x)) return LONG_MAX;
306 273997 : for (vx = 0;; vx++)
307 273997 : if (signe(gel(x,2+vx))) break;
308 222321 : return vx;
309 : }
310 : long
311 31525307 : ZX_valrem(GEN x, GEN *Z)
312 : {
313 : long vx;
314 31525307 : if (!signe(x)) { *Z = pol_0(varn(x)); return LONG_MAX; }
315 52249287 : for (vx = 0;; vx++)
316 52249287 : if (signe(gel(x,2+vx))) break;
317 31525307 : *Z = RgX_shift_shallow(x, -vx);
318 31525302 : return vx;
319 : }
320 :
321 : GEN
322 21 : ZX_div_by_X_1(GEN a, GEN *r)
323 : {
324 21 : long l = lg(a), i;
325 21 : GEN a0, z0, z = cgetg(l-1, t_POL);
326 21 : z[1] = a[1];
327 21 : a0 = a + l-1;
328 21 : z0 = z + l-2; *z0 = *a0--;
329 4732 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + z[i+1] */
330 : {
331 4711 : GEN t = addii(gel(a0--,0), gel(z0--,0));
332 4711 : gel(z0,0) = t;
333 : }
334 21 : if (r) *r = addii(gel(a0,0), gel(z0,0));
335 21 : return z;
336 : }
337 :
338 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
339 : static GEN
340 3906557 : ZX_Z_translate_basecase(GEN P, GEN c)
341 : {
342 3906557 : pari_sp av = avma;
343 : GEN Q, R;
344 : long i, k, n;
345 :
346 3906557 : if (!signe(P) || !signe(c)) return ZX_copy(P);
347 3860709 : Q = leafcopy(P);
348 3860749 : R = Q+2; n = degpol(P);
349 3860748 : if (equali1(c))
350 : {
351 18433761 : for (i=1; i<=n; i++)
352 : {
353 101457753 : for (k=n-i; k<n; k++) gel(R,k) = addii(gel(R,k), gel(R,k+1));
354 15475107 : if (gc_needed(av,2))
355 : {
356 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_translate(1), i = %ld/%ld", i,n);
357 0 : Q = gc_GEN(av, Q); R = Q+2;
358 : }
359 : }
360 : }
361 902101 : else if (equalim1(c))
362 : {
363 957346 : for (i=1; i<=n; i++)
364 : {
365 3180652 : for (k=n-i; k<n; k++) gel(R,k) = subii(gel(R,k), gel(R,k+1));
366 730688 : if (gc_needed(av,2))
367 : {
368 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_translate(-1), i = %ld/%ld", i,n);
369 0 : Q = gc_GEN(av, Q); R = Q+2;
370 : }
371 : }
372 : }
373 : else
374 : {
375 4706571 : for (i=1; i<=n; i++)
376 : {
377 26626774 : for (k=n-i; k<n; k++) gel(R,k) = addmulii_inplace(gel(R,k), c, gel(R,k+1));
378 4031128 : if (gc_needed(av,2))
379 : {
380 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_translate, i = %ld/%ld", i,n);
381 0 : Q = gc_GEN(av, Q); R = Q+2;
382 : }
383 : }
384 : }
385 3860390 : return gc_GEN(av, Q);
386 : }
387 :
388 : GEN
389 434 : Xpm1_powu(long n, long s, long v)
390 : {
391 : long d, k;
392 : GEN C;
393 434 : if (!n) return pol_1(v);
394 434 : d = (n + 1) >> 1;
395 434 : C = cgetg(n+3, t_POL);
396 434 : C[1] = evalsigne(1)| evalvarn(v);
397 434 : gel(C,2) = gen_1;
398 434 : gel(C,3) = utoipos(n);
399 1771 : for (k=2; k <= d; k++)
400 1337 : gel(C,k+2) = diviuexact(mului(n-k+1, gel(C,k+1)), k);
401 434 : if (s < 0)
402 0 : for (k = odd(n)? 0: 1; k <= d; k += 2)
403 0 : togglesign_safe(&gel(C,k+2));
404 434 : if (s > 0 || !odd(n))
405 2058 : for (k = d+1; k <= n; k++) gel(C,k+2) = gel(C,n-k+2);
406 : else
407 0 : for (k = d+1; k <= n; k++) gel(C,k+2) = negi(gel(C,n-k+2));
408 434 : return C;
409 : }
410 : /* return (x+u)^n */
411 : static GEN
412 0 : Z_XpN_powu(GEN u, long n, long v)
413 : {
414 : pari_sp av;
415 : long k;
416 : GEN B, C, V;
417 0 : if (!n) return pol_1(v);
418 0 : if (is_pm1(u))
419 0 : return Xpm1_powu(n, signe(u), v);
420 0 : av = avma;
421 0 : V = gpowers(u, n);
422 0 : B = vecbinomial(n);
423 0 : C = cgetg(n+3, t_POL);
424 0 : C[1] = evalsigne(1)| evalvarn(v);
425 0 : for (k=1; k <= n+1; k++)
426 0 : gel(C,k+1) = mulii(gel(V,n+2-k), gel(B,k));
427 0 : return gc_upto(av, C);
428 : }
429 :
430 : GEN
431 4000342 : ZX_Z_translate(GEN P, GEN c)
432 : {
433 4000342 : pari_sp av = avma;
434 : long n;
435 4000342 : if (!signe(c)) return ZX_copy(P);
436 3906562 : n = degpol(P);
437 3906559 : if (n < 220)
438 3906558 : return ZX_Z_translate_basecase(P, c);
439 : else
440 : {
441 1 : long d = n >> 1;
442 1 : GEN Q = ZX_Z_translate(RgX_shift_shallow(P, -d), c);
443 0 : GEN R = ZX_Z_translate(RgXn_red_shallow(P, d), c);
444 0 : GEN S = Z_XpN_powu(c, d, varn(P));
445 0 : return gc_upto(av, ZX_add(ZX_mul(Q, S), R));
446 : }
447 : }
448 :
449 : /* P(ax + b) */
450 : GEN
451 690868 : ZX_affine(GEN P, GEN a, GEN b)
452 : {
453 690868 : if (signe(b)) P = ZX_Z_translate(P, b);
454 690868 : return ZX_unscale(P, a);
455 : }
456 :
457 : GEN
458 524759 : ZX_Z_eval(GEN x, GEN y)
459 : {
460 524759 : long i = lg(x)-1, j;
461 524759 : pari_sp av = avma;
462 : GEN t, r;
463 :
464 524759 : if (i<=2) return (i==2)? icopy(gel(x,2)): gen_0;
465 486189 : if (!signe(y)) return icopy(gel(x,2));
466 :
467 486077 : t = gel(x,i); i--;
468 : #if 0 /* standard Horner's rule */
469 : for ( ; i>=2; i--)
470 : t = addii(mulii(t,y),gel(x,i));
471 : #endif
472 : /* specific attention to sparse polynomials */
473 2403141 : for ( ; i>=2; i = j-1)
474 : {
475 2100729 : for (j = i; !signe(gel(x,j)); j--)
476 183665 : if (j==2)
477 : {
478 3661 : if (i != j) y = powiu(y, i-j+1);
479 3661 : return gc_INT(av, mulii(t,y));
480 : }
481 1917064 : r = (i==j)? y: powiu(y, i-j+1);
482 1917064 : t = addii(mulii(t,r), gel(x,j));
483 1917064 : if (gc_needed(av,2))
484 : {
485 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_eval: i = %ld",i);
486 0 : t = gc_INT(av, t);
487 : }
488 : }
489 482416 : return gc_INT(av, t);
490 : }
491 :
492 : /* Return 2^(n degpol(P)) P(x >> n) */
493 : GEN
494 0 : ZX_rescale2n(GEN P, long n)
495 : {
496 0 : long i, l = lg(P), ni = n;
497 : GEN Q;
498 0 : if (l==2) return pol_0(varn(P));
499 0 : Q = cgetg(l,t_POL);
500 0 : gel(Q,l-1) = icopy(gel(P,l-1));
501 0 : for (i=l-2; i>=2; i--)
502 : {
503 0 : gel(Q,i) = shifti(gel(P,i), ni);
504 0 : ni += n;
505 : }
506 0 : Q[1] = P[1]; return Q;
507 : }
508 :
509 : /* Return h^deg(P) P(x / h), not memory clean. h integer, P ZX */
510 : GEN
511 36283 : ZX_rescale(GEN P, GEN h)
512 : {
513 36283 : long l = lg(P);
514 36283 : GEN Q = cgetg(l,t_POL);
515 36281 : if (l != 2)
516 : {
517 36281 : long i = l-1;
518 36281 : GEN hi = h;
519 36281 : gel(Q,i) = gel(P,i);
520 36281 : if (l != 3) { i--; gel(Q,i) = mulii(gel(P,i), h); }
521 114146 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
522 : }
523 36293 : Q[1] = P[1]; return Q;
524 : }
525 : /* Return h^(deg(P)-1) P(x / h), P!=0, h=lt(P), memory unclean; monic result */
526 : GEN
527 0 : ZX_rescale_lt(GEN P)
528 : {
529 0 : long l = lg(P);
530 0 : GEN Q = cgetg(l,t_POL);
531 0 : gel(Q,l-1) = gen_1;
532 0 : if (l != 3)
533 : {
534 0 : long i = l-1;
535 0 : GEN h = gel(P,i), hi = h;
536 0 : i--; gel(Q,i) = gel(P,i);
537 0 : if (l != 4) { i--; gel(Q,i) = mulii(gel(P,i), h); }
538 0 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
539 : }
540 0 : Q[1] = P[1]; return Q;
541 : }
542 :
543 : /*Eval x in 2^(k*BIL) in linear time*/
544 : static GEN
545 95689102 : ZX_eval2BILspec(GEN x, long k, long nx)
546 : {
547 95689102 : pari_sp av = avma;
548 95689102 : long i,j, lz = k*nx, ki;
549 95689102 : GEN pz = cgetipos(2+lz);
550 95679554 : GEN nz = cgetipos(2+lz);
551 4301952429 : for(i=0; i < lz; i++)
552 : {
553 4206260296 : *int_W(pz,i) = 0UL;
554 4206260296 : *int_W(nz,i) = 0UL;
555 : }
556 1287733126 : for(i=0, ki=0; i<nx; i++, ki+=k)
557 : {
558 1192040993 : GEN c = gel(x,i);
559 1192040993 : long lc = lgefint(c)-2;
560 1192040993 : if (signe(c)==0) continue;
561 970538673 : if (signe(c) > 0)
562 2355685080 : for (j=0; j<lc; j++) *int_W(pz,j+ki) = *int_W(c,j);
563 : else
564 359603481 : for (j=0; j<lc; j++) *int_W(nz,j+ki) = *int_W(c,j);
565 : }
566 95692133 : pz = int_normalize(pz,0);
567 95703105 : nz = int_normalize(nz,0); return gc_INT(av, subii(pz,nz));
568 : }
569 :
570 : static long
571 172366757 : ZX_expispec(GEN x, long nx)
572 : {
573 172366757 : long i, m = 0;
574 1560888729 : for(i = 0; i < nx; i++)
575 : {
576 1388547612 : long e = expi(gel(x,i));
577 1388521972 : if (e > m) m = e;
578 : }
579 172341117 : return m;
580 : }
581 :
582 : static GEN
583 57471739 : Z_mod2BIL_ZX(GEN x, long bs, long d, long vx)
584 : {
585 57471739 : long i, offset, lm = lgefint(x)-2, l = d+vx+3, sx = signe(x);
586 57471739 : GEN s1 = int2n(bs*BITS_IN_LONG), pol = cgetg(l, t_POL);
587 57451210 : int carry = 0;
588 57451210 : pol[1] = evalsigne(1);
589 78292972 : for (i=0; i<vx; i++) gel(pol,i+2) = gen_0;
590 1320255480 : for (offset=0; i <= d+vx; i++, offset += bs)
591 : {
592 1262944868 : pari_sp av = avma;
593 1262944868 : long lz = minss(bs, lm-offset);
594 1262831271 : GEN z = lz > 0 ?adduispec_offset(carry, x, offset, lz): utoi(carry);
595 1262707251 : if (lgefint(z) == 3+bs) { carry = 1; z = gen_0;}
596 : else
597 : {
598 1248928743 : carry = (lgefint(z) == 2+bs && (HIGHBIT & *int_W(z,bs-1)));
599 1248928743 : if (carry)
600 59549839 : z = gc_INT(av, (sx==-1)? subii(s1,z): subii(z,s1));
601 1189378904 : else if (sx==-1) togglesign(z);
602 : }
603 1262804270 : gel(pol,i+2) = z;
604 : }
605 57310612 : return ZX_renormalize(pol,l);
606 : }
607 :
608 : static GEN
609 18038535 : ZX_sqrspec_sqri(GEN x, long nx, long ex, long v)
610 : {
611 18038535 : long e = 2*ex + expu(nx) + 3;
612 18038459 : long N = divsBIL(e)+1;
613 18038364 : GEN z = sqri(ZX_eval2BILspec(x,N,nx));
614 18035856 : return Z_mod2BIL_ZX(z, N, nx*2-2, v);
615 : }
616 :
617 : static GEN
618 38181697 : ZX_mulspec_mulii(GEN x, GEN y, long nx, long ny, long ex, long ey, long v)
619 : {
620 38181697 : long e = ex + ey + expu(minss(nx,ny)) + 3;
621 38181034 : long N = divsBIL(e)+1;
622 38180813 : GEN z = mulii(ZX_eval2BILspec(x,N,nx), ZX_eval2BILspec(y,N,ny));
623 38176470 : return Z_mod2BIL_ZX(z, N, nx+ny-2, v);
624 : }
625 :
626 : INLINE GEN
627 315914007 : ZX_sqrspec_basecase_limb(GEN x, long a, long i)
628 : {
629 315914007 : pari_sp av = avma;
630 315914007 : GEN s = gen_0;
631 315914007 : long j, l = (i+1)>>1;
632 623299569 : for (j=a; j<l; j++)
633 : {
634 308225314 : GEN xj = gel(x,j), xx = gel(x,i-j);
635 308225314 : if (signe(xj) && signe(xx))
636 305407007 : s = addii(s, mulii(xj, xx));
637 : }
638 315074255 : s = shifti(s,1);
639 315343978 : if ((i&1) == 0)
640 : {
641 196338251 : GEN t = gel(x, i>>1);
642 196338251 : if (signe(t))
643 196105790 : s = addii(s, sqri(t));
644 : }
645 314927980 : return gc_INT(av,s);
646 : }
647 :
648 : static GEN
649 76662231 : ZX_sqrspec_basecase(GEN x, long nx, long v)
650 : {
651 : long i, lz, nz;
652 : GEN z;
653 :
654 76662231 : lz = (nx << 1) + 1; nz = lz-2;
655 76662231 : lz += v;
656 76662231 : z = cgetg(lz,t_POL); z[1] = evalsigne(1); z += 2;
657 76855172 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
658 272900614 : for (i=0; i<nx; i++)
659 196365145 : gel(z,i) = ZX_sqrspec_basecase_limb(x, 0, i);
660 196355873 : for ( ; i<nz; i++) gel(z,i) = ZX_sqrspec_basecase_limb(x, i-nx+1, i);
661 76511832 : z -= v+2; return z;
662 : }
663 :
664 : static GEN
665 43729890 : Z_sqrshiftspec_ZX(GEN x, long vx)
666 : {
667 43729890 : long i, nz = 2*vx+1;
668 43729890 : GEN z = cgetg(2+nz, t_POL);
669 43726839 : z[1] = evalsigne(1);
670 46612832 : for(i=2;i<nz+1;i++) gel(z,i) = gen_0;
671 43726839 : gel(z,nz+1) = sqri(x);
672 43692066 : return z;
673 : }
674 :
675 : static GEN
676 48369887 : Z_ZX_mulshiftspec(GEN x, GEN y, long ny, long vz)
677 : {
678 48369887 : long i, nz = vz+ny;
679 48369887 : GEN z = cgetg(2+nz, t_POL);
680 48369917 : z[1] = evalsigne(1);
681 191634944 : for (i=0; i<vz; i++) gel(z,i+2) = gen_0;
682 140694686 : for (i=0; i<ny; i++) gel(z,i+vz+2) = mulii(x, gel(y,i));
683 48362880 : return z;
684 : }
685 :
686 : GEN
687 139141740 : ZX_sqrspec(GEN x, long nx)
688 : {
689 : #ifdef PARI_KERNEL_GMP
690 80311662 : const long low[]={ 17, 32, 96, 112, 160, 128, 128, 160, 160, 160, 160, 160, 176, 192, 192, 192, 192, 192, 224, 224, 224, 240, 240, 240, 272, 288, 288, 240, 288, 304, 304, 304, 304, 304, 304, 352, 352, 368, 352, 352, 352, 368, 368, 432, 432, 496, 432, 496, 496};
691 80311662 : const long high[]={ 102860, 70254, 52783, 27086, 24623, 18500, 15289, 13899, 12635, 11487, 10442, 9493, 8630, 7845, 7132, 7132, 6484, 6484, 5894, 5894, 4428, 4428, 3660, 4428, 3660, 3660, 2749, 2499, 2272, 2066, 1282, 1282, 1166, 1166, 1166, 1166, 1166, 1166, 1166, 963, 963, 724, 658, 658, 658, 528, 528, 528, 528};
692 : #else
693 58830078 : const long low[]={ 17, 17, 32, 32, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 112, 112, 128, 112, 112, 112, 112, 112, 128, 128, 160, 160, 112, 128, 128, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 160, 176, 160, 160, 176, 160, 160, 176, 176, 208, 176, 176, 176, 192, 192, 176, 176, 224, 176, 224, 224, 176, 224, 224, 224, 176, 176, 176, 176, 176, 176, 176, 176, 224, 176, 176, 224, 224, 224, 224, 224, 224, 224, 240, 288, 240, 288, 288, 240, 288, 288, 240, 240, 304, 304};
694 58830078 : const long high[]={ 165657, 85008, 52783, 43622, 32774, 27086, 22385, 15289, 13899, 12635, 11487, 10442, 9493, 7845, 6484, 6484, 5894, 5894, 4871, 4871, 4428, 4026, 3660, 3660, 3660, 3327, 3327, 3024, 2749, 2749, 2272, 2749, 2499, 2499, 2272, 1878, 1878, 1878, 1707, 1552, 1552, 1552, 1552, 1552, 1411, 1411, 1411, 1282, 1282, 1282, 1282, 1282, 1166, 1166, 1166, 1166, 1166, 1166, 1166, 1166, 1060, 1060, 963, 963, 963, 963, 963, 963, 963, 963, 963, 963, 963, 876, 876, 876, 876, 796, 658, 724, 658, 724, 658, 658, 658, 658, 658, 658, 658, 658, 658, 658, 658, 658, 336, 658, 658, 592, 336, 336};
695 : #endif
696 139141740 : const long nblow = numberof(low);
697 : pari_sp av;
698 : long ex, vx;
699 : GEN z;
700 139141740 : if (!nx) return pol_0(0);
701 138322335 : vx = ZX_valspec(x,nx); nx-=vx; x+=vx;
702 138427134 : if (nx==1) return Z_sqrshiftspec_ZX(gel(x, 0), vx);
703 94696761 : av = avma;
704 94696761 : ex = ZX_expispec(x,nx);
705 94693165 : if (nx-2 < nblow && low[nx-2]<=ex && ex<=high[nx-2])
706 76663129 : z = ZX_sqrspec_basecase(x, nx, 2*vx);
707 : else
708 18030036 : z = ZX_sqrspec_sqri(x, nx, ex, 2*vx);
709 94624073 : return gc_upto(av, z);
710 : }
711 :
712 : GEN
713 139149151 : ZX_sqr(GEN x)
714 : {
715 139149151 : GEN z = ZX_sqrspec(x+2, lgpol(x));
716 139040204 : z[1] = x[1];
717 139040204 : return z;
718 : }
719 :
720 : GEN
721 124196828 : ZX_mulspec(GEN x, GEN y, long nx, long ny)
722 : {
723 : pari_sp av;
724 : long ex, ey, vx, vy, v;
725 124196828 : if (!nx || !ny) return pol_0(0);
726 99233747 : vx = ZX_valspec(x,nx); nx-=vx; x += vx;
727 99234374 : vy = ZX_valspec(y,ny); ny-=vy; y += vy;
728 99237394 : v = vx + vy;
729 99237394 : if (nx==1) return Z_ZX_mulshiftspec(gel(x,0), y, ny, v);
730 60997504 : if (ny==1) return Z_ZX_mulshiftspec(gel(y,0), x, nx, v);
731 50867409 : if (nx == 2 && ny == 2)
732 : {
733 12686174 : GEN a0 = gel(x,0), a1 = gel(x,1), A0, A1, A2;
734 12686174 : GEN b0 = gel(y,0), b1 = gel(y,1), z = cgetg(5 + v, t_POL);
735 : long i;
736 12686192 : z[1] = evalvarn(0) | evalsigne(1);
737 12686192 : A0 = mulii(a0, b0);
738 12685968 : A2 = mulii(a1, b1); av = avma;
739 12685886 : A1 = gc_INT(av, subii(addii(A0,A2),
740 : mulii(subii(a1,a0), subii(b1,b0))));
741 12686119 : i = 4 + v;
742 12686119 : gel(z,i--) = A2;
743 12686119 : gel(z,i--) = A1;
744 14777863 : gel(z,i--) = A0; while (i > 1) gel(z, i--) = gen_0;
745 12686119 : return z;
746 : }
747 : #if 0
748 : /* generically slower even when degrees differ a lot; sometimes about twice
749 : * faster when bitsize is moderate */
750 : if (DEBUGVAR)
751 : return RgX_mulspec(x - vx, y - vy, nx + vx, ny + vy);
752 : #endif
753 38181235 : av = avma;
754 38181235 : ex = ZX_expispec(x, nx);
755 38181139 : ey = ZX_expispec(y, ny);
756 38181224 : return gc_upto(av, ZX_mulspec_mulii(x,y,nx,ny,ex,ey,v));
757 : }
758 : GEN
759 115740035 : ZX_mul(GEN x, GEN y)
760 : {
761 : GEN z;
762 115740035 : if (x == y) return ZX_sqr(x);
763 115001353 : z = ZX_mulspec(x+2,y+2,lgpol(x),lgpol(y));
764 114995197 : z[1] = x[1];
765 114995197 : if (!signe(y)) z[1] &= VARNBITS;
766 114995197 : return z;
767 : }
768 :
769 : /* x,y two ZX in the same variable; assume y is monic */
770 : GEN
771 11654657 : ZX_rem(GEN x, GEN y)
772 : {
773 : long vx, dx, dy, dz, i, j, sx, lr;
774 : pari_sp av0, av;
775 : GEN z,p1,rem;
776 :
777 11654657 : vx = varn(x);
778 11654657 : dy = degpol(y);
779 11654548 : dx = degpol(x);
780 11654373 : if (dx < dy) return ZX_copy(x);
781 6179819 : if (!dy) return pol_0(vx); /* y is constant */
782 6177054 : av0 = avma; dz = dx-dy;
783 6177054 : z=cgetg(dz+3,t_POL); z[1] = x[1];
784 6177047 : x += 2; y += 2; z += 2;
785 :
786 6177047 : p1 = gel(x,dx);
787 6177047 : gel(z,dz) = icopy(p1);
788 31570455 : for (i=dx-1; i>=dy; i--)
789 : {
790 25393318 : av=avma; p1=gel(x,i);
791 278650773 : for (j=i-dy+1; j<=i && j<=dz; j++)
792 253294275 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
793 25356498 : gel(z,i-dy) = avma == av? icopy(p1): gc_INT(av, p1);
794 : }
795 6177137 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
796 9328296 : for (sx=0; ; i--)
797 : {
798 9328296 : p1 = gel(x,i);
799 65353265 : for (j=0; j<=i && j<=dz; j++)
800 56022900 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
801 9330365 : if (signe(p1)) { sx = 1; break; }
802 3533824 : if (!i) break;
803 3153467 : set_avma(av);
804 : }
805 6176898 : lr=i+3; rem -= lr;
806 6176898 : rem[0] = evaltyp(t_POL) | _evallg(lr);
807 6176898 : rem[1] = z[-1];
808 6176898 : p1 = gc_INT((pari_sp)rem, p1);
809 6177186 : rem += 2; gel(rem,i) = p1;
810 37319835 : for (i--; i>=0; i--)
811 : {
812 31141692 : av=avma; p1 = gel(x,i);
813 320977702 : for (j=0; j<=i && j<=dz; j++)
814 289900743 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
815 31076959 : gel(rem,i) = avma == av? icopy(p1): gc_INT(av, p1);
816 : }
817 6178143 : rem -= 2;
818 6178143 : if (!sx) (void)ZX_renormalize(rem, lr);
819 6178143 : return gc_upto(av0,rem);
820 : }
821 :
822 : /* return x(1) */
823 : GEN
824 4732 : ZX_eval1(GEN x)
825 : {
826 4732 : pari_sp av = avma;
827 4732 : long i = lg(x)-1;
828 : GEN s;
829 4732 : if (i < 2) return gen_0;
830 4207 : s = gel(x,i); i--;
831 4207 : if (i == 1) return icopy(s);
832 34671 : for ( ; i>=2; i--)
833 : {
834 31472 : GEN c = gel(x,i);
835 31472 : if (signe(c)) s = addii(s, c);
836 : }
837 3199 : return gc_INT(av,s);
838 : }
839 :
840 : /* reduce T mod X^n - 1. Shallow function */
841 : GEN
842 1735765 : ZX_mod_Xnm1(GEN T, ulong n)
843 : {
844 1735765 : long i, j, L = lg(T), l = n+2;
845 : GEN S;
846 1735765 : if (L <= l) return T;
847 1573066 : S = cgetg(l, t_POL);
848 1572769 : S[1] = T[1];
849 12137668 : for (i = 2; i < l; i++) gel(S,i) = gel(T,i);
850 7387028 : for (j = 2; i < L; i++) {
851 5821745 : gel(S,j) = addii(gel(S,j), gel(T,i));
852 5814259 : if (++j == l) j = 2;
853 : }
854 1565283 : return normalizepol_lg(S, l);
855 : }
856 :
857 : static GEN
858 0 : _ZX_mul(void* E, GEN x, GEN y)
859 0 : { (void) E; return ZX_mul(x, y); }
860 : static GEN
861 266 : _ZX_sqr(void *E, GEN x)
862 266 : { (void) E; return ZX_sqr(x); }
863 :
864 : static GEN
865 735 : _ZX_divrem(void * E, GEN x, GEN y, GEN *r)
866 735 : { (void) E; return RgX_divrem(x, y, r); }
867 :
868 : static GEN
869 0 : _ZX_add(void * E, GEN x, GEN y)
870 0 : { (void) E; return ZX_add(x, y); }
871 :
872 : static struct bb_ring ZX_ring = { _ZX_add,_ZX_mul,_ZX_sqr };
873 :
874 : GEN
875 518 : ZX_digits(GEN x, GEN T)
876 : {
877 518 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
878 518 : if (signe(x)==0) return(cgetg(1, t_VEC));
879 280 : return gen_digits(x, T, n, NULL, &ZX_ring, _ZX_divrem);
880 : }
881 :
882 : GEN
883 0 : ZXV_ZX_fromdigits(GEN x, GEN T)
884 0 : { return gen_fromdigits(x,T, NULL, &ZX_ring); }
885 :
886 : /*******************************************************************/
887 : /* */
888 : /* ZXV */
889 : /* */
890 : /*******************************************************************/
891 :
892 : int
893 28 : ZXV_equal(GEN V, GEN W)
894 : {
895 28 : long l = lg(V);
896 28 : if (l!=lg(W)) return 0;
897 28 : while (--l > 0)
898 28 : if (!ZX_equal(gel(V,l), gel(W,l))) return 0;
899 0 : return 1;
900 : }
901 :
902 : GEN
903 7476 : ZXV_Z_mul(GEN x, GEN y)
904 29904 : { pari_APPLY_same(ZX_Z_mul(gel(x,i), y)) }
905 :
906 : GEN
907 0 : ZXV_remi2n(GEN x, long N)
908 0 : { pari_APPLY_same(ZX_remi2n(gel(x,i), N)) }
909 :
910 : GEN
911 74242 : ZXV_dotproduct(GEN x, GEN y)
912 : {
913 74242 : pari_sp av = avma;
914 74242 : long i, lx = lg(x);
915 : GEN c;
916 74242 : if (lx == 1) return pol_0(varn(x));
917 74242 : c = ZX_mul(gel(x,1), gel(y,1));
918 364182 : for (i = 2; i < lx; i++)
919 : {
920 289940 : GEN t = ZX_mul(gel(x,i), gel(y,i));
921 289940 : if (signe(t)) c = ZX_add(c, t);
922 : }
923 74242 : return gc_upto(av, c);
924 : }
925 :
926 : /*******************************************************************/
927 : /* */
928 : /* ZXQM */
929 : /* */
930 : /*******************************************************************/
931 :
932 : GEN
933 2600127 : ZXn_mul(GEN x, GEN y, long n)
934 2600127 : { return RgXn_red_shallow(ZX_mul(x, y), n); }
935 :
936 : GEN
937 1575 : ZXn_sqr(GEN x, long n)
938 1575 : { return RgXn_red_shallow(ZX_sqr(x), n); }
939 :
940 : /*******************************************************************/
941 : /* */
942 : /* ZXQM */
943 : /* */
944 : /*******************************************************************/
945 :
946 : static long
947 5557926 : ZX_expi(GEN x)
948 : {
949 5557926 : if (signe(x)==0) return 0;
950 2883419 : if (typ(x)==t_INT) return expi(x);
951 1304771 : return ZX_expispec(x+2, lgpol(x));
952 : }
953 :
954 : static long
955 1092747 : ZXC_expi(GEN x)
956 : {
957 1092747 : long i, l = lg(x), m=0;
958 6650673 : for(i = 1; i < l; i++)
959 : {
960 5557926 : long e = ZX_expi(gel(x,i));
961 5557926 : if (e > m) m = e;
962 : }
963 1092747 : return m;
964 : }
965 :
966 : static long
967 177100 : ZXM_expi(GEN x)
968 : {
969 177100 : long i, l = lg(x), m=0;
970 1269847 : for(i = 1; i < l; i++)
971 : {
972 1092747 : long e = ZXC_expi(gel(x,i));
973 1092747 : if (e > m) m = e;
974 : }
975 177100 : return m;
976 : }
977 :
978 : static GEN
979 5557926 : ZX_eval2BIL(GEN x, long k)
980 : {
981 5557926 : if (signe(x)==0) return gen_0;
982 2883419 : if (typ(x)==t_INT) return x;
983 1304771 : return ZX_eval2BILspec(x+2, k, lgpol(x));
984 : }
985 :
986 : /*Eval x in 2^(k*BIL) in linear time*/
987 : static GEN
988 1092747 : ZXC_eval2BIL(GEN x, long k)
989 : {
990 1092747 : long i, lx = lg(x);
991 1092747 : GEN A = cgetg(lx, t_COL);
992 6650673 : for (i=1; i<lx; i++) gel(A,i) = ZX_eval2BIL(gel(x,i), k);
993 1092747 : return A;
994 : }
995 :
996 : static GEN
997 177100 : ZXM_eval2BIL(GEN x, long k)
998 : {
999 177100 : long i, lx = lg(x);
1000 177100 : GEN A = cgetg(lx, t_MAT);
1001 1269847 : for (i=1; i<lx; i++) gel(A,i) = ZXC_eval2BIL(gel(x,i), k);
1002 177100 : return A;
1003 : }
1004 :
1005 : static GEN
1006 65338 : Z_mod2BIL_ZXQ(GEN x, long bs, GEN T)
1007 : {
1008 65338 : pari_sp av = avma;
1009 65338 : long v = varn(T), d = 2*(degpol(T)-1);
1010 65338 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1011 65338 : setvarn(z, v);
1012 65338 : return gc_upto(av, ZX_rem(z, T));
1013 : }
1014 :
1015 : static GEN
1016 14896 : ZC_mod2BIL_ZXQC(GEN x, long bs, GEN T)
1017 : {
1018 14896 : long i, lx = lg(x);
1019 14896 : GEN A = cgetg(lx, t_COL);
1020 80234 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_ZXQ(gel(x,i), bs, T);
1021 14896 : return A;
1022 : }
1023 :
1024 : static GEN
1025 6769 : ZM_mod2BIL_ZXQM(GEN x, long bs, GEN T)
1026 : {
1027 6769 : long i, lx = lg(x);
1028 6769 : GEN A = cgetg(lx, t_MAT);
1029 21665 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_ZXQC(gel(x,i), bs, T);
1030 6769 : return A;
1031 : }
1032 :
1033 : GEN
1034 6629 : ZXQM_mul(GEN x, GEN y, GEN T)
1035 : {
1036 6629 : long d = degpol(T);
1037 : GEN z;
1038 6629 : pari_sp av = avma;
1039 6629 : if (d == 0)
1040 0 : z = ZM_mul(simplify_shallow(x),simplify_shallow(y));
1041 : else
1042 : {
1043 6629 : long e, N, ex = ZXM_expi(x), ey = ZXM_expi(y), n = lg(x)-1;
1044 6629 : e = ex + ey + expu(d) + expu(n) + 4;
1045 6629 : N = divsBIL(e)+1;
1046 6629 : z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1047 6629 : z = ZM_mod2BIL_ZXQM(z, N, T);
1048 : }
1049 6629 : return gc_upto(av, z);
1050 : }
1051 :
1052 : GEN
1053 140 : ZXQM_sqr(GEN x, GEN T)
1054 : {
1055 140 : long d = degpol(T);
1056 : GEN z;
1057 140 : pari_sp av = avma;
1058 140 : if (d == 0)
1059 0 : z = ZM_sqr(simplify_shallow(x));
1060 : else
1061 : {
1062 140 : long ex = ZXM_expi(x), d = degpol(T), n = lg(x)-1;
1063 140 : long e = 2*ex + expu(d) + expu(n) + 4;
1064 140 : long N = divsBIL(e)+1;
1065 140 : z = ZM_sqr(ZXM_eval2BIL(x,N));
1066 140 : z = ZM_mod2BIL_ZXQM(z, N, T);
1067 : }
1068 140 : return gc_upto(av, z);
1069 : }
1070 :
1071 : GEN
1072 5124 : QXQM_mul(GEN x, GEN y, GEN T)
1073 : {
1074 5124 : GEN dx, nx = Q_primitive_part(x, &dx);
1075 5124 : GEN dy, ny = Q_primitive_part(y, &dy);
1076 5124 : GEN z = ZXQM_mul(nx, ny, T);
1077 5124 : if (dx || dy)
1078 : {
1079 5124 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1080 5124 : if (!gequal1(d)) z = RgM_Rg_mul(z, d);
1081 : }
1082 5124 : return z;
1083 : }
1084 :
1085 : GEN
1086 7 : QXQM_sqr(GEN x, GEN T)
1087 : {
1088 7 : GEN dx, nx = Q_primitive_part(x, &dx);
1089 7 : GEN z = ZXQM_sqr(nx, T);
1090 7 : if (dx) z = RgM_Rg_mul(z, gsqr(dx));
1091 7 : return z;
1092 : }
1093 :
1094 : static GEN
1095 1199133 : Z_mod2BIL_Fq(GEN x, long bs, GEN T, GEN p)
1096 : {
1097 1199133 : pari_sp av = avma;
1098 1199133 : long v = get_FpX_var(T), d = 2*(get_FpX_degree(T)-1);
1099 1199133 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1100 1199133 : setvarn(z, v);
1101 1199133 : return gc_upto(av, FpX_rem(z, T, p));
1102 : }
1103 :
1104 : static GEN
1105 305402 : ZC_mod2BIL_FqC(GEN x, long bs, GEN T, GEN p)
1106 : {
1107 305402 : long i, lx = lg(x);
1108 305402 : GEN A = cgetg(lx, t_COL);
1109 1504535 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_Fq(gel(x,i), bs, T, p);
1110 305402 : return A;
1111 : }
1112 :
1113 : static GEN
1114 81851 : ZM_mod2BIL_FqM(GEN x, long bs, GEN T, GEN p)
1115 : {
1116 81851 : long i, lx = lg(x);
1117 81851 : GEN A = cgetg(lx, t_MAT);
1118 387253 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_FqC(gel(x,i), bs, T, p);
1119 81851 : return A;
1120 : }
1121 :
1122 : GEN
1123 81851 : FqM_mul_Kronecker(GEN x, GEN y, GEN T, GEN p)
1124 : {
1125 81851 : pari_sp av = avma;
1126 81851 : long ex = ZXM_expi(x), ey = ZXM_expi(y), d = get_FpX_degree(T), n = lg(x)-1;
1127 81851 : long e = ex + ey + expu(d) + expu(n) + 4;
1128 81851 : long N = divsBIL(e)+1;
1129 81851 : GEN z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1130 81851 : return gc_upto(av, ZM_mod2BIL_FqM(z, N, T, p));
1131 : }
1132 :
1133 : /*******************************************************************/
1134 : /* */
1135 : /* ZXX */
1136 : /* */
1137 : /*******************************************************************/
1138 :
1139 : void
1140 1022 : RgX_check_ZXX(GEN x, const char *s)
1141 : {
1142 1022 : long k = lg(x)-1;
1143 9891 : for ( ; k>1; k--) {
1144 8869 : GEN t = gel(x,k);
1145 8869 : switch(typ(t)) {
1146 7833 : case t_INT: break;
1147 1036 : case t_POL: if (RgX_is_ZX(t)) break;
1148 : /* fall through */
1149 0 : default: pari_err_TYPE(stack_strcat(s, " not in Z[X,Y]"),x);
1150 : }
1151 : }
1152 1022 : }
1153 :
1154 : /*Renormalize (in place) polynomial with t_INT or ZX coefficients.*/
1155 : GEN
1156 316466479 : ZXX_renormalize(GEN x, long lx)
1157 : {
1158 : long i;
1159 387657567 : for (i = lx-1; i>1; i--)
1160 370901536 : if (signe(gel(x,i))) break;
1161 316466479 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + (i+1)));
1162 316459819 : setlg(x, i+1); setsigne(x, i!=1); return x;
1163 : }
1164 :
1165 : GEN
1166 7877 : ZXX_evalx0(GEN y)
1167 : {
1168 7877 : long i, l = lg(y);
1169 7877 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1170 79425 : for(i=2; i<l; i++)
1171 : {
1172 71548 : GEN yi = gel(y,i);
1173 71548 : gel(z,i) = typ(yi)==t_INT? yi: constant_coeff(yi);
1174 : }
1175 7877 : return ZX_renormalize(z,l);
1176 : }
1177 :
1178 : long
1179 5599 : ZXX_max_lg(GEN x)
1180 : {
1181 5599 : long i, prec = 0, lx = lg(x);
1182 39091 : for (i=2; i<lx; i++)
1183 : {
1184 33492 : GEN p = gel(x,i);
1185 33492 : long l = (typ(p) == t_INT)? lgefint(p): ZX_max_lg(p);
1186 33492 : if (l > prec) prec = l;
1187 : }
1188 5599 : return prec;
1189 : }
1190 :
1191 : GEN
1192 8883 : ZXX_Z_mul(GEN y, GEN x)
1193 : {
1194 8883 : long i, l = lg(y);
1195 8883 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1196 230727 : for(i=2; i<l; i++)
1197 221844 : if(typ(gel(y,i))==t_INT)
1198 0 : gel(z,i) = mulii(gel(y,i),x);
1199 : else
1200 221844 : gel(z,i) = ZX_Z_mul(gel(y,i),x);
1201 8883 : return z;
1202 : }
1203 :
1204 : GEN
1205 826 : ZXX_shifti(GEN y, long n)
1206 : {
1207 826 : long i, l = lg(y);
1208 826 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1209 3157 : for(i=2; i<l; i++)
1210 2331 : if(typ(gel(y,i))==t_INT)
1211 42 : gel(z,i) = shifti(gel(y,i),n);
1212 : else
1213 2289 : gel(z,i) = ZX_shifti(gel(y,i), n);
1214 826 : return z;
1215 : }
1216 :
1217 : GEN
1218 0 : ZXX_Z_add_shallow(GEN x, GEN y)
1219 : {
1220 0 : long i, l = lg(x);
1221 : GEN z, a;
1222 0 : if (signe(x)==0) return scalarpol(y,varn(x));
1223 0 : z = cgetg(l,t_POL); z[1] = x[1];
1224 0 : a = gel(x,2);
1225 0 : gel(z, 2) = typ(a)==t_INT? addii(a,y): ZX_Z_add(a,y);
1226 0 : for(i=3; i<l; i++)
1227 0 : gel(z,i) = gel(x,i);
1228 0 : return z;
1229 : }
1230 :
1231 : GEN
1232 56469 : ZXX_Z_divexact(GEN y, GEN x)
1233 : {
1234 56469 : long i, l = lg(y);
1235 56469 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1236 467026 : for(i=2; i<l; i++)
1237 410557 : if(typ(gel(y,i))==t_INT)
1238 2989 : gel(z,i) = diviiexact(gel(y,i),x);
1239 : else
1240 407568 : gel(z,i) = ZX_Z_divexact(gel(y,i),x);
1241 56469 : return z;
1242 : }
1243 :
1244 : /* Kronecker substitution, ZXX -> ZX:
1245 : * P(X,Y) = sum_{0<=i<lP} P_i(X) * Y^i, where deg P_i < n.
1246 : * Returns P(X,X^(2n-1)) */
1247 : GEN
1248 10018057 : RgXX_to_Kronecker_spec(GEN P, long lP, long n)
1249 : {
1250 10018057 : long i, k, N = (n<<1) + 1;
1251 : GEN y;
1252 10018057 : if (!lP) return pol_0(0);
1253 9998197 : y = cgetg((N-2)*lP + 2, t_POL) + 2;
1254 44795346 : for (k=i=0; i<lP; i++)
1255 : {
1256 : long j;
1257 44798146 : GEN c = gel(P,i);
1258 44798146 : if (typ(c)!=t_POL)
1259 : {
1260 3252994 : gel(y,k++) = c;
1261 3252994 : j = 3;
1262 : }
1263 : else
1264 : {
1265 41545152 : long l = lg(c);
1266 41545152 : if (l-3 >= n)
1267 0 : pari_err_BUG("RgXX_to_Kronecker, P is not reduced mod Q");
1268 244844352 : for (j=2; j < l; j++) gel(y,k++) = gel(c,j);
1269 : }
1270 44798258 : if (i == lP-1) break;
1271 191846723 : for ( ; j < N; j++) gel(y,k++) = gen_0;
1272 : }
1273 9996047 : y-=2; setlg(y, k+2); y[1] = evalsigne(1); return y;
1274 : }
1275 : GEN
1276 3164 : RgXX_to_Kronecker_var(GEN P0, long n, long vx)
1277 : {
1278 3164 : long i, k, N, lP = lgpol(P0);
1279 : GEN P, y;
1280 3164 : if (!lP) return pol_0(vx);
1281 3164 : P = P0 + 2; N = (n<<1) + 1;
1282 3164 : y = cgetg((N-2)*lP + 2, t_POL) + 2;
1283 99519 : for (k=i=0; i<lP; i++)
1284 : {
1285 : long j;
1286 99519 : GEN c = gel(P,i);
1287 99519 : if (typ(c) != t_POL || varn(c) != vx)
1288 : {
1289 308 : gel(y,k++) = c;
1290 308 : j = 3;
1291 : }
1292 : else
1293 : {
1294 99211 : long l = lg(c);
1295 99211 : if (l-3 >= n)
1296 0 : pari_err_BUG("RgXX_to_Kronecker, P is not reduced mod Q");
1297 576611 : for (j=2; j < l; j++) gel(y,k++) = gel(c,j);
1298 : }
1299 99519 : if (i == lP-1) break;
1300 517433 : for ( ; j < N; j++) gel(y,k++) = gen_0;
1301 : }
1302 3164 : y-=2; setlg(y, k+2); y[1] = evalsigne(1) | evalvarn(vx); return y;
1303 : }
1304 :
1305 : /* shallow, n = deg(T) */
1306 : GEN
1307 1781622 : Kronecker_to_ZXX(GEN z, long n, long v)
1308 : {
1309 1781622 : long i,j,lx,l, N = (n<<1)+1;
1310 : GEN x, t;
1311 1781622 : l = lg(z); lx = (l-2) / (N-2);
1312 1781622 : x = cgetg(lx+3,t_POL);
1313 1781747 : x[1] = z[1];
1314 9040040 : for (i=2; i<lx+2; i++)
1315 : {
1316 7257487 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1317 108151549 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1318 7256708 : z += (N-2);
1319 7256708 : gel(x,i) = ZX_renormalize(t,N);
1320 : }
1321 1782553 : N = (l-2) % (N-2) + 2;
1322 1782553 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1323 5171214 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1324 1782057 : gel(x,i) = ZX_renormalize(t,N);
1325 1782045 : return ZXX_renormalize(x, i+1);
1326 : }
1327 :
1328 : GEN
1329 6753085 : RgXX_to_Kronecker(GEN P, long n)
1330 : {
1331 6753085 : GEN z = RgXX_to_Kronecker_spec(P+2, lgpol(P), n);
1332 6752372 : setvarn(z,varn(P)); return z;
1333 : }
1334 : GEN
1335 3052008 : ZXX_mul_Kronecker(GEN x, GEN y, long n)
1336 3052008 : { return ZX_mul(RgXX_to_Kronecker(x,n), RgXX_to_Kronecker(y,n)); }
1337 :
1338 : GEN
1339 5306 : ZXX_sqr_Kronecker(GEN x, long n)
1340 5306 : { return ZX_sqr(RgXX_to_Kronecker(x,n)); }
1341 :
1342 : /* shallow, n = deg(T) */
1343 : GEN
1344 24269 : Kronecker_to_ZXQX(GEN z, GEN T)
1345 : {
1346 24269 : long i,j,lx,l, N = (degpol(T)<<1)+1;
1347 : GEN x, t;
1348 24269 : l = lg(z); lx = (l-2) / (N-2);
1349 24269 : x = cgetg(lx+3,t_POL);
1350 24269 : x[1] = z[1];
1351 236278 : for (i=2; i<lx+2; i++)
1352 : {
1353 212009 : t = cgetg(N,t_POL); t[1] = T[1];
1354 1898200 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1355 212008 : z += (N-2);
1356 212008 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1357 : }
1358 24269 : N = (l-2) % (N-2) + 2;
1359 24269 : t = cgetg(N,t_POL); t[1] = T[1];
1360 57180 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1361 24269 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1362 24269 : return ZXX_renormalize(x, i+1);
1363 : }
1364 :
1365 : GEN
1366 1904 : ZXQX_sqr(GEN x, GEN T)
1367 : {
1368 1904 : pari_sp av = avma;
1369 1904 : long n = degpol(T);
1370 1904 : GEN z = ZXX_sqr_Kronecker(x, n);
1371 1904 : z = Kronecker_to_ZXQX(z, T);
1372 1904 : return gc_GEN(av, z);
1373 : }
1374 :
1375 : GEN
1376 22365 : ZXQX_mul(GEN x, GEN y, GEN T)
1377 : {
1378 22365 : pari_sp av = avma;
1379 22365 : long n = degpol(T);
1380 22365 : GEN z = ZXX_mul_Kronecker(x, y, n);
1381 22365 : z = Kronecker_to_ZXQX(z, T);
1382 22365 : return gc_GEN(av, z);
1383 : }
1384 :
1385 : GEN
1386 8988 : ZXQX_ZXQ_mul(GEN P, GEN U, GEN T)
1387 : {
1388 : long i, lP;
1389 : GEN res;
1390 8988 : res = cgetg_copy(P, &lP); res[1] = P[1];
1391 47108 : for(i=2; i<lP; i++)
1392 38121 : gel(res,i) = typ(gel(P,i))==t_POL? ZXQ_mul(U, gel(P,i), T)
1393 38121 : : gmul(U, gel(P,i));
1394 8987 : return ZXX_renormalize(res,lP);
1395 : }
1396 :
1397 : GEN
1398 3814400 : QX_mul(GEN x, GEN y)
1399 : {
1400 3814400 : GEN dx, nx = Q_primitive_part(x, &dx);
1401 3814400 : GEN dy, ny = Q_primitive_part(y, &dy);
1402 3814400 : GEN z = ZX_mul(nx, ny);
1403 3814400 : if (dx || dy)
1404 : {
1405 3718825 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1406 3718825 : return ZX_Q_mul(z, d);
1407 : } else
1408 95575 : return z;
1409 : }
1410 :
1411 : GEN
1412 110727 : QX_sqr(GEN x)
1413 : {
1414 110727 : GEN dx, nx = Q_primitive_part(x, &dx);
1415 110726 : GEN z = ZX_sqr(nx);
1416 110728 : if (dx)
1417 99285 : return ZX_Q_mul(z, gsqr(dx));
1418 : else
1419 11443 : return z;
1420 : }
1421 :
1422 : GEN
1423 1902116 : QX_ZX_rem(GEN x, GEN y)
1424 : {
1425 1902116 : pari_sp av = avma;
1426 1902116 : GEN d, nx = Q_primitive_part(x, &d);
1427 1902116 : GEN r = ZX_rem(nx, y);
1428 1902116 : if (d) r = ZX_Q_mul(r, d);
1429 1902116 : return gc_upto(av, r);
1430 : }
1431 :
1432 : GEN
1433 13377 : QXQX_mul(GEN x, GEN y, GEN T)
1434 : {
1435 13377 : GEN dx, nx = Q_primitive_part(x, &dx);
1436 13377 : GEN dy, ny = Q_primitive_part(y, &dy);
1437 13377 : GEN z = ZXQX_mul(nx, ny, T);
1438 13377 : if (dx || dy)
1439 : {
1440 8309 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1441 8309 : return ZXX_Q_mul(z, d);
1442 : } else
1443 5068 : return z;
1444 : }
1445 :
1446 : GEN
1447 1904 : QXQX_sqr(GEN x, GEN T)
1448 : {
1449 1904 : GEN dx, nx = Q_primitive_part(x, &dx);
1450 1904 : GEN z = ZXQX_sqr(nx, T);
1451 1904 : if (dx)
1452 588 : return ZXX_Q_mul(z, gsqr(dx));
1453 : else
1454 1316 : return z;
1455 : }
1456 :
1457 : GEN
1458 476 : QXQX_powers(GEN P, long n, GEN T)
1459 : {
1460 476 : GEN v = cgetg(n+2, t_VEC);
1461 : long i;
1462 476 : gel(v, 1) = pol_1(varn(T));
1463 476 : if (n==0) return v;
1464 476 : gel(v, 2) = gcopy(P);
1465 1554 : for (i = 2; i <= n; i++) gel(v,i+1) = QXQX_mul(P, gel(v,i), T);
1466 476 : return v;
1467 : }
1468 :
1469 : GEN
1470 3465 : QXQX_QXQ_mul(GEN P, GEN U, GEN T)
1471 : {
1472 : long i, lP;
1473 : GEN res;
1474 3465 : res = cgetg_copy(P, &lP); res[1] = P[1];
1475 22715 : for(i=2; i<lP; i++)
1476 19250 : gel(res,i) = typ(gel(P,i))==t_POL? QXQ_mul(U, gel(P,i), T)
1477 19250 : : gmul(U, gel(P,i));
1478 3465 : return ZXX_renormalize(res,lP);
1479 : }
|
