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 350659 : RgX_check_QX(GEN x, const char *s)
25 350659 : { if (!RgX_is_QX(x)) pari_err_TYPE(stack_strcat(s," [not in Q[X]]"), x); }
26 : void
27 2711172 : RgX_check_ZX(GEN x, const char *s)
28 2711172 : { if (!RgX_is_ZX(x)) pari_err_TYPE(stack_strcat(s," [not in Z[X]]"), x); }
29 : long
30 107128 : ZX_max_lg(GEN x)
31 : {
32 107128 : long i, l = 0, lx = lg(x);
33 647692 : for (i = 2; i < lx; i++) l = maxss(l, lgefint(gel(x,i)));
34 107127 : return l;
35 : }
36 :
37 : GEN
38 40627327 : ZX_add(GEN x, GEN y)
39 : {
40 : long lx,ly,i;
41 : GEN z;
42 40627327 : lx = lg(x); ly = lg(y); if (lx < ly) swapspec(x,y, lx,ly);
43 40627327 : z = cgetg(lx,t_POL); z[1] = x[1];
44 248153904 : for (i=2; i<ly; i++) gel(z,i) = addii(gel(x,i),gel(y,i));
45 79365526 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
46 40507432 : if (lx == ly) z = ZX_renormalize(z, lx);
47 40508235 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
48 39292134 : return z;
49 : }
50 :
51 : GEN
52 29598885 : ZX_sub(GEN x,GEN y)
53 : {
54 29598885 : long i, lx = lg(x), ly = lg(y);
55 : GEN z;
56 29598885 : if (lx >= ly)
57 : {
58 29144086 : z = cgetg(lx,t_POL); z[1] = x[1];
59 118400537 : for (i=2; i<ly; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
60 29105825 : if (lx == ly)
61 : {
62 23455125 : z = ZX_renormalize(z, lx);
63 23455256 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); z = pol_0(varn(x)); }
64 : }
65 : else
66 17425682 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
67 : }
68 : else
69 : {
70 454799 : z = cgetg(ly,t_POL); z[1] = y[1];
71 1201592 : for (i=2; i<lx; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
72 1410896 : for ( ; i<ly; i++) gel(z,i) = negi(gel(y,i));
73 : }
74 29597020 : return z;
75 : }
76 :
77 : GEN
78 598428 : ZX_neg(GEN x)
79 : {
80 598428 : long i, l = lg(x);
81 598428 : GEN y = cgetg(l,t_POL);
82 2465731 : y[1] = x[1]; for(i=2; i<l; i++) gel(y,i) = negi(gel(x,i));
83 598478 : return y;
84 : }
85 : GEN
86 5307859 : ZX_copy(GEN x)
87 : {
88 5307859 : long i, l = lg(x);
89 5307859 : GEN y = cgetg(l, t_POL);
90 5307865 : y[1] = x[1];
91 21798866 : for (i=2; i<l; i++)
92 : {
93 16491005 : GEN c = gel(x,i);
94 16491005 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
95 : }
96 5307861 : return y;
97 : }
98 :
99 : GEN
100 56107 : scalar_ZX(GEN x, long v)
101 : {
102 : GEN z;
103 56107 : if (!signe(x)) return pol_0(v);
104 4572 : z = cgetg(3, t_POL);
105 4572 : z[1] = evalsigne(1) | evalvarn(v);
106 4572 : gel(z,2) = icopy(x); return z;
107 : }
108 :
109 : GEN
110 23155 : scalar_ZX_shallow(GEN x, long v)
111 : {
112 : GEN z;
113 23155 : if (!signe(x)) return pol_0(v);
114 19991 : z = cgetg(3, t_POL);
115 19991 : z[1] = evalsigne(1) | evalvarn(v);
116 19991 : gel(z,2) = x; return z;
117 : }
118 :
119 : GEN
120 87405 : ZX_Z_add(GEN y, GEN x)
121 : {
122 : long lz, i;
123 87405 : GEN z = cgetg_copy(y, &lz);
124 87407 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
125 76669 : z[1] = y[1];
126 76669 : gel(z,2) = addii(gel(y,2),x);
127 1784395 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
128 76793 : if (lz==3) z = ZX_renormalize(z,lz);
129 76670 : return z;
130 : }
131 : GEN
132 94388 : ZX_Z_add_shallow(GEN y, GEN x)
133 : {
134 : long lz, i;
135 94388 : GEN z = cgetg_copy(y, &lz);
136 94389 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX_shallow(x,varn(y)); }
137 94158 : z[1] = y[1];
138 94158 : gel(z,2) = addii(gel(y,2),x);
139 552539 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
140 94150 : if (lz==3) z = ZX_renormalize(z,lz);
141 94150 : return z;
142 : }
143 :
144 : GEN
145 181417 : ZX_Z_sub(GEN y, GEN x)
146 : {
147 : long lz, i;
148 181417 : GEN z = cgetg_copy(y, &lz);
149 181417 : 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 180976 : z[1] = y[1];
159 180976 : gel(z,2) = subii(gel(y,2),x);
160 906362 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
161 180976 : if (lz==3) z = ZX_renormalize(z,lz);
162 180976 : return z;
163 : }
164 :
165 : GEN
166 2705437 : Z_ZX_sub(GEN x, GEN y)
167 : {
168 : long lz, i;
169 2705437 : GEN z = cgetg_copy(y, &lz);
170 2705474 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
171 2705474 : z[1] = y[1];
172 2705474 : gel(z,2) = subii(x, gel(y,2));
173 13539712 : for(i=3; i<lz; i++) gel(z,i) = negi(gel(y,i));
174 2705748 : if (lz==3) z = ZX_renormalize(z,lz);
175 2705537 : return z;
176 : }
177 :
178 : GEN
179 5077314 : ZX_Z_divexact(GEN y,GEN x)
180 : {
181 5077314 : long i, l = lg(y);
182 5077314 : GEN z = cgetg(l,t_POL); z[1] = y[1];
183 32437554 : for(i=2; i<l; i++) gel(z,i) = diviiexact(gel(y,i),x);
184 5075013 : return z;
185 : }
186 :
187 : GEN
188 977795 : ZX_divuexact(GEN y, ulong x)
189 : {
190 977795 : long i, l = lg(y);
191 977795 : GEN z = cgetg(l,t_POL); z[1] = y[1];
192 9684353 : for(i=2; i<l; i++) gel(z,i) = diviuexact(gel(y,i),x);
193 977795 : 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 17024213 : ZX_Z_mul(GEN y,GEN x)
207 : {
208 : GEN z;
209 : long i, l;
210 17024213 : if (!signe(x)) return pol_0(varn(y));
211 13979354 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
212 92893928 : for(i=2; i<l; i++) gel(z,i) = mulii(gel(y,i),x);
213 13974403 : return z;
214 : }
215 :
216 : GEN
217 253094 : ZX_mulu(GEN y, ulong x)
218 : {
219 : GEN z;
220 : long i, l;
221 253094 : if (!x) return pol_0(varn(y));
222 253094 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
223 1025517 : for(i=2; i<l; i++) gel(z,i) = mului(x,gel(y,i));
224 253094 : return z;
225 : }
226 :
227 : GEN
228 6419743 : ZX_shifti(GEN y, long n)
229 : {
230 : GEN z;
231 : long i, l;
232 6419743 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
233 30967493 : for(i=2; i<l; i++) gel(z,i) = shifti(gel(y,i),n);
234 6416721 : return ZX_renormalize(z,l);
235 : }
236 :
237 : GEN
238 106035 : ZX_remi2n(GEN y, long n)
239 : {
240 : GEN z;
241 : long i, l;
242 106035 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
243 4714246 : for(i=2; i<l; i++) gel(z,i) = remi2n(gel(y,i),n);
244 104604 : return ZX_renormalize(z,l);
245 : }
246 :
247 : GEN
248 45512 : ZXT_remi2n(GEN z, long n)
249 : {
250 45512 : if (typ(z) == t_POL)
251 39416 : return ZX_remi2n(z, n);
252 : else
253 : {
254 6096 : long i,l = lg(z);
255 6096 : GEN x = cgetg(l, t_VEC);
256 18317 : for (i=1; i<l; i++) gel(x,i) = ZXT_remi2n(gel(z,i), n);
257 6102 : return x;
258 : }
259 : }
260 :
261 : GEN
262 1323 : zx_to_ZX(GEN z)
263 : {
264 1323 : long i, l = lg(z);
265 1323 : GEN x = cgetg(l,t_POL);
266 100912 : for (i=2; i<l; i++) gel(x,i) = stoi(z[i]);
267 1323 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
268 : }
269 :
270 : GEN
271 7399056 : ZX_deriv(GEN x)
272 : {
273 7399056 : long i,lx = lg(x)-1;
274 : GEN y;
275 :
276 7399056 : if (lx<3) return pol_0(varn(x));
277 7316619 : y = cgetg(lx,t_POL);
278 38847744 : for (i=2; i<lx ; i++) gel(y,i) = mului(i-1,gel(x,i+1));
279 7315405 : y[1] = x[1]; return y;
280 : }
281 :
282 : int
283 1956759 : ZX_equal(GEN V, GEN W)
284 : {
285 1956759 : long i, l = lg(V);
286 1956759 : if (lg(W) != l) return 0;
287 5891941 : for (i = 2; i < l; i++)
288 4561458 : if (!equalii(gel(V,i), gel(W,i))) return 0;
289 1330483 : return 1;
290 : }
291 :
292 : static long
293 310072946 : ZX_valspec(GEN x, long nx)
294 : {
295 : long vx;
296 382333846 : for (vx = 0; vx<nx ; vx++)
297 382350689 : if (signe(gel(x,vx))) break;
298 310072946 : return vx;
299 : }
300 :
301 : long
302 330646 : ZX_val(GEN x)
303 : {
304 : long vx;
305 330646 : if (!signe(x)) return LONG_MAX;
306 264939 : for (vx = 0;; vx++)
307 264939 : if (signe(gel(x,2+vx))) break;
308 215013 : return vx;
309 : }
310 : long
311 28567286 : ZX_valrem(GEN x, GEN *Z)
312 : {
313 : long vx;
314 28567286 : if (!signe(x)) { *Z = pol_0(varn(x)); return LONG_MAX; }
315 49274739 : for (vx = 0;; vx++)
316 49274739 : if (signe(gel(x,2+vx))) break;
317 28567286 : *Z = RgX_shift_shallow(x, -vx);
318 28567282 : 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 3279415 : ZX_translate_basecase(GEN P, GEN c)
341 : {
342 3279415 : pari_sp av = avma;
343 : GEN Q, R;
344 : long i, k, n;
345 :
346 3279415 : if (!signe(P) || !signe(c)) return ZX_copy(P);
347 3238407 : Q = leafcopy(P);
348 3238466 : R = Q+2; n = degpol(P);
349 3238466 : if (equali1(c))
350 : {
351 17159029 : for (i=1; i<=n; i++)
352 : {
353 98009156 : for (k=n-i; k<n; k++) gel(R,k) = addii(gel(R,k), gel(R,k+1));
354 14525117 : if (gc_needed(av,2))
355 : {
356 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(1), i = %ld/%ld", i,n);
357 0 : Q = gerepilecopy(av, Q); R = Q+2;
358 : }
359 : }
360 : }
361 604557 : else if (equalim1(c))
362 : {
363 55248 : for (i=1; i<=n; i++)
364 : {
365 181403 : for (k=n-i; k<n; k++) gel(R,k) = subii(gel(R,k), gel(R,k+1));
366 45687 : if (gc_needed(av,2))
367 : {
368 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(-1), i = %ld/%ld", i,n);
369 0 : Q = gerepilecopy(av, Q); R = Q+2;
370 : }
371 : }
372 : }
373 : else
374 : {
375 4359184 : for (i=1; i<=n; i++)
376 : {
377 25467061 : for (k=n-i; k<n; k++) gel(R,k) = addmulii_inplace(gel(R,k), c, gel(R,k+1));
378 3764188 : if (gc_needed(av,2))
379 : {
380 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate, i = %ld/%ld", i,n);
381 0 : Q = gerepilecopy(av, Q); R = Q+2;
382 : }
383 : }
384 : }
385 3238148 : return gerepilecopy(av, Q);
386 : }
387 :
388 : static GEN
389 0 : Z_Xpm1_powu(long n, long s, long v)
390 : {
391 : long d, k;
392 : GEN C;
393 0 : if (!n) return pol_1(v);
394 0 : d = (n + 1) >> 1;
395 0 : C = cgetg(n+3, t_POL);
396 0 : C[1] = evalsigne(1)| evalvarn(v);
397 0 : gel(C,2) = gen_1;
398 0 : gel(C,3) = utoipos(n);
399 0 : for (k=2; k <= d; k++)
400 0 : gel(C,k+2) = diviuexact(mului(n-k+1, gel(C,k+1)), k);
401 0 : if (s < 0)
402 0 : for (k = odd(n)? 0: 1; k <= d; k += 2)
403 0 : togglesign_safe(&gel(C,k+2));
404 0 : if (s > 0 || !odd(n))
405 0 : 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 0 : 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 Z_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 gerepileupto(av, C);
428 : }
429 :
430 : GEN
431 3279422 : ZX_translate(GEN P, GEN c)
432 : {
433 3279422 : pari_sp av = avma;
434 3279422 : long n = degpol(P);
435 3279416 : if (n < 220)
436 3279416 : return ZX_translate_basecase(P, c);
437 : else
438 : {
439 0 : long d = n >> 1;
440 0 : GEN Q = ZX_translate(RgX_shift_shallow(P, -d), c);
441 0 : GEN R = ZX_translate(RgXn_red_shallow(P, d), c);
442 0 : GEN S = Z_XpN_powu(c, d, varn(P));
443 0 : return gerepileupto(av, ZX_add(ZX_mul(Q, S), R));
444 : }
445 : }
446 :
447 : /* P(ax + b) */
448 : GEN
449 991304 : ZX_affine(GEN P, GEN a, GEN b)
450 : {
451 991304 : if (signe(b)) P = ZX_translate(P, b);
452 991304 : return ZX_unscale(P, a);
453 : }
454 :
455 : GEN
456 531208 : ZX_Z_eval(GEN x, GEN y)
457 : {
458 531208 : long i = lg(x)-1, j;
459 531208 : pari_sp av = avma;
460 : GEN t, r;
461 :
462 531208 : if (i<=2) return (i==2)? icopy(gel(x,2)): gen_0;
463 493086 : if (!signe(y)) return icopy(gel(x,2));
464 :
465 492967 : t = gel(x,i); i--;
466 : #if 0 /* standard Horner's rule */
467 : for ( ; i>=2; i--)
468 : t = addii(mulii(t,y),gel(x,i));
469 : #endif
470 : /* specific attention to sparse polynomials */
471 2437298 : for ( ; i>=2; i = j-1)
472 : {
473 2120645 : for (j = i; !signe(gel(x,j)); j--)
474 176319 : if (j==2)
475 : {
476 2975 : if (i != j) y = powiu(y, i-j+1);
477 2975 : return gerepileuptoint(av, mulii(t,y));
478 : }
479 1944326 : r = (i==j)? y: powiu(y, i-j+1);
480 1944347 : t = addii(mulii(t,r), gel(x,j));
481 1944330 : if (gc_needed(av,2))
482 : {
483 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_eval: i = %ld",i);
484 0 : t = gerepileuptoint(av, t);
485 : }
486 : }
487 489997 : return gerepileuptoint(av, t);
488 : }
489 :
490 : /* Return 2^(n degpol(P)) P(x >> n) */
491 : GEN
492 0 : ZX_rescale2n(GEN P, long n)
493 : {
494 0 : long i, l = lg(P), ni = n;
495 : GEN Q;
496 0 : if (l==2) return pol_0(varn(P));
497 0 : Q = cgetg(l,t_POL);
498 0 : gel(Q,l-1) = icopy(gel(P,l-1));
499 0 : for (i=l-2; i>=2; i--)
500 : {
501 0 : gel(Q,i) = shifti(gel(P,i), ni);
502 0 : ni += n;
503 : }
504 0 : Q[1] = P[1]; return Q;
505 : }
506 :
507 : /* Return h^deg(P) P(x / h), not memory clean. h integer, P ZX */
508 : GEN
509 35632 : ZX_rescale(GEN P, GEN h)
510 : {
511 35632 : long l = lg(P);
512 35632 : GEN Q = cgetg(l,t_POL);
513 35635 : if (l != 2)
514 : {
515 35635 : long i = l-1;
516 35635 : GEN hi = h;
517 35635 : gel(Q,i) = gel(P,i);
518 35635 : if (l != 3) { i--; gel(Q,i) = mulii(gel(P,i), h); }
519 110282 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
520 : }
521 35640 : Q[1] = P[1]; return Q;
522 : }
523 : /* Return h^(deg(P)-1) P(x / h), P!=0, h=lt(P), memory unclean; monic result */
524 : GEN
525 0 : ZX_rescale_lt(GEN P)
526 : {
527 0 : long l = lg(P);
528 0 : GEN Q = cgetg(l,t_POL);
529 0 : gel(Q,l-1) = gen_1;
530 0 : if (l != 3)
531 : {
532 0 : long i = l-1;
533 0 : GEN h = gel(P,i), hi = h;
534 0 : i--; gel(Q,i) = gel(P,i);
535 0 : if (l != 4) { i--; gel(Q,i) = mulii(gel(P,i), h); }
536 0 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
537 : }
538 0 : Q[1] = P[1]; return Q;
539 : }
540 :
541 : /*Eval x in 2^(k*BIL) in linear time*/
542 : static GEN
543 92496440 : ZX_eval2BILspec(GEN x, long k, long nx)
544 : {
545 92496440 : pari_sp av = avma;
546 92496440 : long i,j, lz = k*nx, ki;
547 92496440 : GEN pz = cgetipos(2+lz);
548 92491781 : GEN nz = cgetipos(2+lz);
549 3274270308 : for(i=0; i < lz; i++)
550 : {
551 3181765671 : *int_W(pz,i) = 0UL;
552 3181765671 : *int_W(nz,i) = 0UL;
553 : }
554 997450374 : for(i=0, ki=0; i<nx; i++, ki+=k)
555 : {
556 904945737 : GEN c = gel(x,i);
557 904945737 : long lc = lgefint(c)-2;
558 904945737 : if (signe(c)==0) continue;
559 681044954 : if (signe(c) > 0)
560 1718638970 : for (j=0; j<lc; j++) *int_W(pz,j+ki) = *int_W(c,j);
561 : else
562 339512171 : for (j=0; j<lc; j++) *int_W(nz,j+ki) = *int_W(c,j);
563 : }
564 92504637 : pz = int_normalize(pz,0);
565 92513311 : nz = int_normalize(nz,0); return gerepileuptoint(av, subii(pz,nz));
566 : }
567 :
568 : static long
569 168702928 : ZX_expispec(GEN x, long nx)
570 : {
571 168702928 : long i, m = 0;
572 1269474007 : for(i = 0; i < nx; i++)
573 : {
574 1100805169 : long e = expi(gel(x,i));
575 1100771079 : if (e > m) m = e;
576 : }
577 168668838 : return m;
578 : }
579 :
580 : static GEN
581 55744280 : Z_mod2BIL_ZX(GEN x, long bs, long d, long vx)
582 : {
583 55744280 : long i, offset, lm = lgefint(x)-2, l = d+vx+3, sx = signe(x);
584 55744280 : GEN s1 = int2n(bs*BITS_IN_LONG), pol = cgetg(l, t_POL);
585 55724969 : int carry = 0;
586 55724969 : pol[1] = evalsigne(1);
587 76209125 : for (i=0; i<vx; i++) gel(pol,i+2) = gen_0;
588 993674655 : for (offset=0; i <= d+vx; i++, offset += bs)
589 : {
590 938081153 : pari_sp av = avma;
591 938081153 : long lz = minss(bs, lm-offset);
592 937955460 : GEN z = lz > 0 ?adduispec_offset(carry, x, offset, lz): utoi(carry);
593 937906320 : if (lgefint(z) == 3+bs) { carry = 1; z = gen_0;}
594 : else
595 : {
596 924454442 : carry = (lgefint(z) == 2+bs && (HIGHBIT & *int_W(z,bs-1)));
597 924454442 : if (carry)
598 57112563 : z = gerepileuptoint(av, (sx==-1)? subii(s1,z): subii(z,s1));
599 867341879 : else if (sx==-1) togglesign(z);
600 : }
601 937949686 : gel(pol,i+2) = z;
602 : }
603 55593502 : return ZX_renormalize(pol,l);
604 : }
605 :
606 : static GEN
607 17785670 : ZX_sqrspec_sqri(GEN x, long nx, long ex, long v)
608 : {
609 17785670 : long e = 2*ex + expu(nx) + 3;
610 17785579 : long N = divsBIL(e)+1;
611 17785489 : GEN z = sqri(ZX_eval2BILspec(x,N,nx));
612 17783174 : return Z_mod2BIL_ZX(z, N, nx*2-2, v);
613 : }
614 :
615 : static GEN
616 36671287 : ZX_mulspec_mulii(GEN x, GEN y, long nx, long ny, long ex, long ey, long v)
617 : {
618 36671287 : long e = ex + ey + expu(minss(nx,ny)) + 3;
619 36670829 : long N = divsBIL(e)+1;
620 36670852 : GEN z = mulii(ZX_eval2BILspec(x,N,nx), ZX_eval2BILspec(y,N,ny));
621 36667598 : return Z_mod2BIL_ZX(z, N, nx+ny-2, v);
622 : }
623 :
624 : INLINE GEN
625 315214906 : ZX_sqrspec_basecase_limb(GEN x, long a, long i)
626 : {
627 315214906 : pari_sp av = avma;
628 315214906 : GEN s = gen_0;
629 315214906 : long j, l = (i+1)>>1;
630 622559799 : for (j=a; j<l; j++)
631 : {
632 308261133 : GEN xj = gel(x,j), xx = gel(x,i-j);
633 308261133 : if (signe(xj) && signe(xx))
634 304156745 : s = addii(s, mulii(xj, xx));
635 : }
636 314298666 : s = shifti(s,1);
637 314591555 : if ((i&1) == 0)
638 : {
639 195736212 : GEN t = gel(x, i>>1);
640 195736212 : if (signe(t))
641 195451514 : s = addii(s, sqri(t));
642 : }
643 314169452 : return gerepileuptoint(av,s);
644 : }
645 :
646 : static GEN
647 76188333 : ZX_sqrspec_basecase(GEN x, long nx, long v)
648 : {
649 : long i, lz, nz;
650 : GEN z;
651 :
652 76188333 : lz = (nx << 1) + 1; nz = lz-2;
653 76188333 : lz += v;
654 76188333 : z = cgetg(lz,t_POL); z[1] = evalsigne(1); z += 2;
655 76382504 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
656 271833136 : for (i=0; i<nx; i++)
657 195798568 : gel(z,i) = ZX_sqrspec_basecase_limb(x, 0, i);
658 195760906 : for ( ; i<nz; i++) gel(z,i) = ZX_sqrspec_basecase_limb(x, i-nx+1, i);
659 76014226 : z -= v+2; return z;
660 : }
661 :
662 : static GEN
663 43593242 : Z_sqrshiftspec_ZX(GEN x, long vx)
664 : {
665 43593242 : long i, nz = 2*vx+1;
666 43593242 : GEN z = cgetg(2+nz, t_POL);
667 43590083 : z[1] = evalsigne(1);
668 46174789 : for(i=2;i<nz+1;i++) gel(z,i) = gen_0;
669 43590083 : gel(z,nz+1) = sqri(x);
670 43552441 : return z;
671 : }
672 :
673 : static GEN
674 36706486 : Z_ZX_mulshiftspec(GEN x, GEN y, long ny, long vz)
675 : {
676 36706486 : long i, nz = vz+ny;
677 36706486 : GEN z = cgetg(2+nz, t_POL);
678 36706507 : z[1] = evalsigne(1);
679 85170028 : for (i=0; i<vz; i++) gel(z,i+2) = gen_0;
680 112911504 : for (i=0; i<ny; i++) gel(z,i+vz+2) = mulii(x, gel(y,i));
681 36699161 : return z;
682 : }
683 :
684 : GEN
685 138234425 : ZX_sqrspec(GEN x, long nx)
686 : {
687 : #ifdef PARI_KERNEL_GMP
688 79490183 : 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};
689 79490183 : 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};
690 : #else
691 58744242 : 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};
692 58744242 : 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};
693 : #endif
694 138234425 : const long nblow = numberof(low);
695 : pari_sp av;
696 : long ex, vx;
697 : GEN z;
698 138234425 : if (!nx) return pol_0(0);
699 137444171 : vx = ZX_valspec(x,nx); nx-=vx; x+=vx;
700 137563484 : if (nx==1) return Z_sqrshiftspec_ZX(gel(x, 0), vx);
701 93970528 : av = avma;
702 93970528 : ex = ZX_expispec(x,nx);
703 93964792 : if (nx-2 < nblow && low[nx-2]<=ex && ex<=high[nx-2])
704 76189259 : z = ZX_sqrspec_basecase(x, nx, 2*vx);
705 : else
706 17775533 : z = ZX_sqrspec_sqri(x, nx, ex, 2*vx);
707 93879983 : return gerepileupto(av, z);
708 : }
709 :
710 : GEN
711 138238356 : ZX_sqr(GEN x)
712 : {
713 138238356 : GEN z = ZX_sqrspec(x+2, lgpol(x));
714 138119105 : z[1] = x[1];
715 138119105 : return z;
716 : }
717 :
718 : GEN
719 96053058 : ZX_mulspec(GEN x, GEN y, long nx, long ny)
720 : {
721 : pari_sp av;
722 : long ex, ey, vx, vy, v;
723 96053058 : if (!nx || !ny) return pol_0(0);
724 86270297 : vx = ZX_valspec(x,nx); nx-=vx; x += vx;
725 86270466 : vy = ZX_valspec(y,ny); ny-=vy; y += vy;
726 86273665 : v = vx + vy;
727 86273665 : if (nx==1) return Z_ZX_mulshiftspec(gel(x,0), y, ny, v);
728 59451722 : if (ny==1) return Z_ZX_mulshiftspec(gel(y,0), x, nx, v);
729 49567118 : if (nx == 2 && ny == 2)
730 : {
731 12895067 : GEN a0 = gel(x,0), a1 = gel(x,1), A0, A1, A2;
732 12895067 : GEN b0 = gel(y,0), b1 = gel(y,1), z = cgetg(5 + v, t_POL);
733 : long i;
734 12895114 : z[1] = evalvarn(0) | evalsigne(1);
735 12895114 : A0 = mulii(a0, b0);
736 12894829 : A2 = mulii(a1, b1); av = avma;
737 12894710 : A1 = gerepileuptoint(av, subii(addii(A0,A2),
738 : mulii(subii(a1,a0), subii(b1,b0))));
739 12894956 : i = 4 + v;
740 12894956 : gel(z,i--) = A2;
741 12894956 : gel(z,i--) = A1;
742 14975159 : gel(z,i--) = A0; while (i > 1) gel(z, i--) = gen_0;
743 12894956 : return z;
744 : }
745 : #if 0
746 : /* generically slower even when degrees differ a lot; sometimes about twice
747 : * faster when bitsize is moderate */
748 : if (DEBUGVAR)
749 : return RgX_mulspec(x - vx, y - vy, nx + vx, ny + vy);
750 : #endif
751 36672051 : av = avma;
752 36672051 : ex = ZX_expispec(x, nx);
753 36670412 : ey = ZX_expispec(y, ny);
754 36670538 : return gerepileupto(av, ZX_mulspec_mulii(x,y,nx,ny,ex,ey,v));
755 : }
756 : GEN
757 88793014 : ZX_mul(GEN x, GEN y)
758 : {
759 : GEN z;
760 88793014 : if (x == y) return ZX_sqr(x);
761 88056544 : z = ZX_mulspec(x+2,y+2,lgpol(x),lgpol(y));
762 88049571 : z[1] = x[1];
763 88049571 : if (!signe(y)) z[1] &= VARNBITS;
764 88049571 : return z;
765 : }
766 :
767 : /* x,y two ZX in the same variable; assume y is monic */
768 : GEN
769 10973404 : ZX_rem(GEN x, GEN y)
770 : {
771 : long vx, dx, dy, dz, i, j, sx, lr;
772 : pari_sp av0, av;
773 : GEN z,p1,rem;
774 :
775 10973404 : vx = varn(x);
776 10973404 : dy = degpol(y);
777 10973338 : dx = degpol(x);
778 10973293 : if (dx < dy) return ZX_copy(x);
779 6039641 : if (!dy) return pol_0(vx); /* y is constant */
780 6036876 : av0 = avma; dz = dx-dy;
781 6036876 : z=cgetg(dz+3,t_POL); z[1] = x[1];
782 6036907 : x += 2; y += 2; z += 2;
783 :
784 6036907 : p1 = gel(x,dx);
785 6036907 : gel(z,dz) = icopy(p1);
786 28571090 : for (i=dx-1; i>=dy; i--)
787 : {
788 22534083 : av=avma; p1=gel(x,i);
789 212743271 : for (j=i-dy+1; j<=i && j<=dz; j++)
790 190221613 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
791 22521658 : gel(z,i-dy) = avma == av? icopy(p1): gerepileuptoint(av, p1);
792 : }
793 6037007 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
794 8932996 : for (sx=0; ; i--)
795 : {
796 8932996 : p1 = gel(x,i);
797 57424675 : for (j=0; j<=i && j<=dz; j++)
798 48491855 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
799 8932820 : if (signe(p1)) { sx = 1; break; }
800 3272152 : if (!i) break;
801 2896012 : set_avma(av);
802 : }
803 6036808 : lr=i+3; rem -= lr;
804 6036808 : rem[0] = evaltyp(t_POL) | _evallg(lr);
805 6036808 : rem[1] = z[-1];
806 6036808 : p1 = gerepileuptoint((pari_sp)rem, p1);
807 6037033 : rem += 2; gel(rem,i) = p1;
808 33748124 : for (i--; i>=0; i--)
809 : {
810 27710170 : av=avma; p1 = gel(x,i);
811 249466510 : for (j=0; j<=i && j<=dz; j++)
812 221791579 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
813 27674931 : gel(rem,i) = avma == av? icopy(p1): gerepileuptoint(av, p1);
814 : }
815 6037954 : rem -= 2;
816 6037954 : if (!sx) (void)ZX_renormalize(rem, lr);
817 6037954 : return gerepileupto(av0,rem);
818 : }
819 :
820 : /* return x(1) */
821 : GEN
822 4718 : ZX_eval1(GEN x)
823 : {
824 4718 : pari_sp av = avma;
825 4718 : long i = lg(x)-1;
826 : GEN s;
827 4718 : if (i < 2) return gen_0;
828 4193 : s = gel(x,i); i--;
829 4193 : if (i == 1) return icopy(s);
830 34608 : for ( ; i>=2; i--)
831 : {
832 31423 : GEN c = gel(x,i);
833 31423 : if (signe(c)) s = addii(s, c);
834 : }
835 3185 : return gerepileuptoint(av,s);
836 : }
837 :
838 : /* reduce T mod X^n - 1. Shallow function */
839 : GEN
840 1735407 : ZX_mod_Xnm1(GEN T, ulong n)
841 : {
842 1735407 : long i, j, L = lg(T), l = n+2;
843 : GEN S;
844 1735407 : if (L <= l) return T;
845 1572758 : S = cgetg(l, t_POL);
846 1572494 : S[1] = T[1];
847 12135632 : for (i = 2; i < l; i++) gel(S,i) = gel(T,i);
848 7382000 : for (j = 2; i < L; i++) {
849 5818439 : gel(S,j) = addii(gel(S,j), gel(T,i));
850 5809506 : if (++j == l) j = 2;
851 : }
852 1563561 : return normalizepol_lg(S, l);
853 : }
854 :
855 : static GEN
856 0 : _ZX_mul(void* E, GEN x, GEN y)
857 0 : { (void) E; return ZX_mul(x, y); }
858 : static GEN
859 0 : _ZX_sqr(void *E, GEN x)
860 0 : { (void) E; return ZX_sqr(x); }
861 :
862 : static GEN
863 0 : _ZX_divrem(void * E, GEN x, GEN y, GEN *r)
864 0 : { (void) E; return RgX_divrem(x, y, r); }
865 :
866 : static GEN
867 0 : _ZX_add(void * E, GEN x, GEN y)
868 0 : { (void) E; return ZX_add(x, y); }
869 :
870 : static struct bb_ring ZX_ring = { _ZX_add,_ZX_mul,_ZX_sqr };
871 :
872 : GEN
873 0 : ZX_digits(GEN x, GEN T)
874 : {
875 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
876 0 : return gen_digits(x, T, n, NULL, &ZX_ring, _ZX_divrem);
877 : }
878 :
879 : GEN
880 0 : ZXV_ZX_fromdigits(GEN x, GEN T)
881 0 : { return gen_fromdigits(x,T, NULL, &ZX_ring); }
882 :
883 : /*******************************************************************/
884 : /* */
885 : /* ZXV */
886 : /* */
887 : /*******************************************************************/
888 :
889 : int
890 28 : ZXV_equal(GEN V, GEN W)
891 : {
892 28 : long l = lg(V);
893 28 : if (l!=lg(W)) return 0;
894 28 : while (--l > 0)
895 28 : if (!ZX_equal(gel(V,l), gel(W,l))) return 0;
896 0 : return 1;
897 : }
898 :
899 : GEN
900 7476 : ZXV_Z_mul(GEN x, GEN y)
901 29904 : { pari_APPLY_same(ZX_Z_mul(gel(x,i), y)) }
902 :
903 : GEN
904 0 : ZXV_remi2n(GEN x, long N)
905 0 : { pari_APPLY_same(ZX_remi2n(gel(x,i), N)) }
906 :
907 : GEN
908 82943 : ZXV_dotproduct(GEN x, GEN y)
909 : {
910 82943 : pari_sp av = avma;
911 82943 : long i, lx = lg(x);
912 : GEN c;
913 82943 : if (lx == 1) return pol_0(varn(x));
914 82943 : c = ZX_mul(gel(x,1), gel(y,1));
915 392175 : for (i = 2; i < lx; i++)
916 : {
917 309232 : GEN t = ZX_mul(gel(x,i), gel(y,i));
918 309232 : if (signe(t)) c = ZX_add(c, t);
919 : }
920 82943 : return gerepileupto(av, c);
921 : }
922 :
923 : GEN
924 122324 : ZXC_to_FlxC(GEN x, ulong p, long sv)
925 5044919 : { pari_APPLY_type(t_COL,typ(gel(x,i))==t_INT ?
926 : Z_to_Flx(gel(x,i), p, sv): ZX_to_Flx(gel(x,i), p)) }
927 :
928 : GEN
929 11592 : ZXM_to_FlxM(GEN x, ulong p, long sv)
930 102181 : { pari_APPLY_same(ZXC_to_FlxC(gel(x,i), p, sv)) }
931 :
932 : /*******************************************************************/
933 : /* */
934 : /* ZXQM */
935 : /* */
936 : /*******************************************************************/
937 :
938 : GEN
939 2598203 : ZXn_mul(GEN x, GEN y, long n)
940 2598203 : { return RgXn_red_shallow(ZX_mul(x, y), n); }
941 :
942 : GEN
943 1407 : ZXn_sqr(GEN x, long n)
944 1407 : { return RgXn_red_shallow(ZX_sqr(x), n); }
945 :
946 : /*******************************************************************/
947 : /* */
948 : /* ZXQM */
949 : /* */
950 : /*******************************************************************/
951 :
952 : static long
953 5624230 : ZX_expi(GEN x)
954 : {
955 5624230 : if (signe(x)==0) return 0;
956 2969420 : if (typ(x)==t_INT) return expi(x);
957 1384655 : return ZX_expispec(x+2, lgpol(x));
958 : }
959 :
960 : static long
961 1092397 : ZXC_expi(GEN x)
962 : {
963 1092397 : long i, l = lg(x), m=0;
964 6716627 : for(i = 1; i < l; i++)
965 : {
966 5624230 : long e = ZX_expi(gel(x,i));
967 5624230 : if (e > m) m = e;
968 : }
969 1092397 : return m;
970 : }
971 :
972 : static long
973 176162 : ZXM_expi(GEN x)
974 : {
975 176162 : long i, l = lg(x), m=0;
976 1268559 : for(i = 1; i < l; i++)
977 : {
978 1092397 : long e = ZXC_expi(gel(x,i));
979 1092397 : if (e > m) m = e;
980 : }
981 176162 : return m;
982 : }
983 :
984 : static GEN
985 5624230 : ZX_eval2BIL(GEN x, long k)
986 : {
987 5624230 : if (signe(x)==0) return gen_0;
988 2969420 : if (typ(x)==t_INT) return x;
989 1384655 : return ZX_eval2BILspec(x+2, k, lgpol(x));
990 : }
991 :
992 : /*Eval x in 2^(k*BIL) in linear time*/
993 : static GEN
994 1092397 : ZXC_eval2BIL(GEN x, long k)
995 : {
996 1092397 : long i, lx = lg(x);
997 1092397 : GEN A = cgetg(lx, t_COL);
998 6716627 : for (i=1; i<lx; i++) gel(A,i) = ZX_eval2BIL(gel(x,i), k);
999 1092397 : return A;
1000 : }
1001 :
1002 : static GEN
1003 176162 : ZXM_eval2BIL(GEN x, long k)
1004 : {
1005 176162 : long i, lx = lg(x);
1006 176162 : GEN A = cgetg(lx, t_MAT);
1007 1268559 : for (i=1; i<lx; i++) gel(A,i) = ZXC_eval2BIL(gel(x,i), k);
1008 176162 : return A;
1009 : }
1010 :
1011 : static GEN
1012 63112 : Z_mod2BIL_ZXQ(GEN x, long bs, GEN T)
1013 : {
1014 63112 : pari_sp av = avma;
1015 63112 : long v = varn(T), d = 2*(degpol(T)-1);
1016 63112 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1017 63112 : setvarn(z, v);
1018 63112 : return gerepileupto(av, ZX_rem(z, T));
1019 : }
1020 :
1021 : static GEN
1022 13692 : ZC_mod2BIL_ZXQC(GEN x, long bs, GEN T)
1023 : {
1024 13692 : long i, lx = lg(x);
1025 13692 : GEN A = cgetg(lx, t_COL);
1026 76804 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_ZXQ(gel(x,i), bs, T);
1027 13692 : return A;
1028 : }
1029 :
1030 : static GEN
1031 6279 : ZM_mod2BIL_ZXQM(GEN x, long bs, GEN T)
1032 : {
1033 6279 : long i, lx = lg(x);
1034 6279 : GEN A = cgetg(lx, t_MAT);
1035 19971 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_ZXQC(gel(x,i), bs, T);
1036 6279 : return A;
1037 : }
1038 :
1039 : GEN
1040 6139 : ZXQM_mul(GEN x, GEN y, GEN T)
1041 : {
1042 6139 : long d = degpol(T);
1043 : GEN z;
1044 6139 : pari_sp av = avma;
1045 6139 : if (d == 0)
1046 0 : z = ZM_mul(simplify_shallow(x),simplify_shallow(y));
1047 : else
1048 : {
1049 6139 : long e, N, ex = ZXM_expi(x), ey = ZXM_expi(y), n = lg(x)-1;
1050 6139 : e = ex + ey + expu(d) + expu(n) + 4;
1051 6139 : N = divsBIL(e)+1;
1052 6139 : z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1053 6139 : z = ZM_mod2BIL_ZXQM(z, N, T);
1054 : }
1055 6139 : return gerepileupto(av, z);
1056 : }
1057 :
1058 : GEN
1059 140 : ZXQM_sqr(GEN x, GEN T)
1060 : {
1061 140 : long d = degpol(T);
1062 : GEN z;
1063 140 : pari_sp av = avma;
1064 140 : if (d == 0)
1065 0 : z = ZM_sqr(simplify_shallow(x));
1066 : else
1067 : {
1068 140 : long ex = ZXM_expi(x), d = degpol(T), n = lg(x)-1;
1069 140 : long e = 2*ex + expu(d) + expu(n) + 4;
1070 140 : long N = divsBIL(e)+1;
1071 140 : z = ZM_sqr(ZXM_eval2BIL(x,N));
1072 140 : z = ZM_mod2BIL_ZXQM(z, N, T);
1073 : }
1074 140 : return gerepileupto(av, z);
1075 : }
1076 :
1077 : GEN
1078 4690 : QXQM_mul(GEN x, GEN y, GEN T)
1079 : {
1080 4690 : GEN dx, nx = Q_primitive_part(x, &dx);
1081 4690 : GEN dy, ny = Q_primitive_part(y, &dy);
1082 4690 : GEN z = ZXQM_mul(nx, ny, T);
1083 4690 : if (dx || dy)
1084 : {
1085 4690 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1086 4690 : if (!gequal1(d)) z = RgM_Rg_mul(z, d);
1087 : }
1088 4690 : return z;
1089 : }
1090 :
1091 : GEN
1092 7 : QXQM_sqr(GEN x, GEN T)
1093 : {
1094 7 : GEN dx, nx = Q_primitive_part(x, &dx);
1095 7 : GEN z = ZXQM_sqr(nx, T);
1096 7 : if (dx) z = RgM_Rg_mul(z, gsqr(dx));
1097 7 : return z;
1098 : }
1099 :
1100 : static GEN
1101 1235610 : Z_mod2BIL_Fq(GEN x, long bs, GEN T, GEN p)
1102 : {
1103 1235610 : pari_sp av = avma;
1104 1235610 : long v = get_FpX_var(T), d = 2*(get_FpX_degree(T)-1);
1105 1235610 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1106 1235610 : setvarn(z, v);
1107 1235610 : return gerepileupto(av, FpX_rem(z, T, p));
1108 : }
1109 :
1110 : static GEN
1111 306473 : ZC_mod2BIL_FqC(GEN x, long bs, GEN T, GEN p)
1112 : {
1113 306473 : long i, lx = lg(x);
1114 306473 : GEN A = cgetg(lx, t_COL);
1115 1542083 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_Fq(gel(x,i), bs, T, p);
1116 306473 : return A;
1117 : }
1118 :
1119 : static GEN
1120 81872 : ZM_mod2BIL_FqM(GEN x, long bs, GEN T, GEN p)
1121 : {
1122 81872 : long i, lx = lg(x);
1123 81872 : GEN A = cgetg(lx, t_MAT);
1124 388345 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_FqC(gel(x,i), bs, T, p);
1125 81872 : return A;
1126 : }
1127 :
1128 : GEN
1129 81872 : FqM_mul_Kronecker(GEN x, GEN y, GEN T, GEN p)
1130 : {
1131 81872 : pari_sp av = avma;
1132 81872 : long ex = ZXM_expi(x), ey = ZXM_expi(y), d= degpol(T), n = lg(x)-1;
1133 81872 : long e = ex + ey + expu(d) + expu(n) + 4;
1134 81872 : long N = divsBIL(e)+1;
1135 81872 : GEN z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1136 81872 : return gerepileupto(av, ZM_mod2BIL_FqM(z, N, T, p));
1137 : }
1138 :
1139 : /*******************************************************************/
1140 : /* */
1141 : /* ZXX */
1142 : /* */
1143 : /*******************************************************************/
1144 :
1145 : void
1146 994 : RgX_check_ZXX(GEN x, const char *s)
1147 : {
1148 994 : long k = lg(x)-1;
1149 9765 : for ( ; k>1; k--) {
1150 8771 : GEN t = gel(x,k);
1151 8771 : switch(typ(t)) {
1152 7749 : case t_INT: break;
1153 1022 : case t_POL: if (RgX_is_ZX(t)) break;
1154 : /* fall through */
1155 0 : default: pari_err_TYPE(stack_strcat(s, " not in Z[X,Y]"),x);
1156 : }
1157 : }
1158 994 : }
1159 :
1160 : /*Renormalize (in place) polynomial with t_INT or ZX coefficients.*/
1161 : GEN
1162 302221631 : ZXX_renormalize(GEN x, long lx)
1163 : {
1164 : long i;
1165 372636831 : for (i = lx-1; i>1; i--)
1166 358509855 : if (signe(gel(x,i))) break;
1167 302221631 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + (i+1)));
1168 302216002 : setlg(x, i+1); setsigne(x, i!=1); return x;
1169 : }
1170 :
1171 : GEN
1172 28548 : ZXX_evalx0(GEN y)
1173 : {
1174 28548 : long i, l = lg(y);
1175 28548 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1176 480693 : for(i=2; i<l; i++)
1177 : {
1178 452145 : GEN yi = gel(y,i);
1179 452145 : gel(z,i) = typ(yi)==t_INT? yi: constant_coeff(yi);
1180 : }
1181 28548 : return ZX_renormalize(z,l);
1182 : }
1183 :
1184 : long
1185 5380 : ZXX_max_lg(GEN x)
1186 : {
1187 5380 : long i, prec = 0, lx = lg(x);
1188 37578 : for (i=2; i<lx; i++)
1189 : {
1190 32198 : GEN p = gel(x,i);
1191 32198 : long l = (typ(p) == t_INT)? lgefint(p): ZX_max_lg(p);
1192 32198 : if (l > prec) prec = l;
1193 : }
1194 5380 : return prec;
1195 : }
1196 :
1197 : GEN
1198 8911 : ZXX_Z_mul(GEN y, GEN x)
1199 : {
1200 8911 : long i, l = lg(y);
1201 8911 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1202 230839 : for(i=2; i<l; i++)
1203 221928 : if(typ(gel(y,i))==t_INT)
1204 0 : gel(z,i) = mulii(gel(y,i),x);
1205 : else
1206 221928 : gel(z,i) = ZX_Z_mul(gel(y,i),x);
1207 8911 : return z;
1208 : }
1209 :
1210 : GEN
1211 0 : ZXX_Z_add_shallow(GEN x, GEN y)
1212 : {
1213 0 : long i, l = lg(x);
1214 : GEN z, a;
1215 0 : if (signe(x)==0) return scalarpol(y,varn(x));
1216 0 : z = cgetg(l,t_POL); z[1] = x[1];
1217 0 : a = gel(x,2);
1218 0 : gel(z, 2) = typ(a)==t_INT? addii(a,y): ZX_Z_add(a,y);
1219 0 : for(i=3; i<l; i++)
1220 0 : gel(z,i) = gel(x,i);
1221 0 : return z;
1222 : }
1223 :
1224 : GEN
1225 55209 : ZXX_Z_divexact(GEN y, GEN x)
1226 : {
1227 55209 : long i, l = lg(y);
1228 55209 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1229 463008 : for(i=2; i<l; i++)
1230 407799 : if(typ(gel(y,i))==t_INT)
1231 2989 : gel(z,i) = diviiexact(gel(y,i),x);
1232 : else
1233 404810 : gel(z,i) = ZX_Z_divexact(gel(y,i),x);
1234 55209 : return z;
1235 : }
1236 :
1237 : /* Kronecker substitution, ZXX -> ZX:
1238 : * P(X,Y) = sum_{0<=i<lP} P_i(X) * Y^i, where deg P_i < n.
1239 : * Returns P(X,X^(2n-1)) */
1240 : GEN
1241 10024775 : RgXX_to_Kronecker_spec(GEN P, long lP, long n)
1242 : {
1243 10024775 : long i, k, N = (n<<1) + 1;
1244 : GEN y;
1245 10024775 : if (!lP) return pol_0(0);
1246 10005476 : y = cgetg((N-2)*lP + 2, t_POL) + 2;
1247 45208139 : for (k=i=0; i<lP; i++)
1248 : {
1249 : long j;
1250 45210986 : GEN c = gel(P,i);
1251 45210986 : if (typ(c)!=t_POL)
1252 : {
1253 3208815 : gel(y,k++) = c;
1254 3208815 : j = 3;
1255 : }
1256 : else
1257 : {
1258 42002171 : long l = lg(c);
1259 42002171 : if (l-3 >= n)
1260 0 : pari_err_BUG("RgXX_to_Kronecker, P is not reduced mod Q");
1261 247838488 : for (j=2; j < l; j++) gel(y,k++) = gel(c,j);
1262 : }
1263 45211101 : if (i == lP-1) break;
1264 195771438 : for ( ; j < N; j++) gel(y,k++) = gen_0;
1265 : }
1266 10003984 : y-=2; setlg(y, k+2); y[1] = evalsigne(1); return y;
1267 : }
1268 :
1269 : /* shallow, n = deg(T) */
1270 : GEN
1271 1767321 : Kronecker_to_ZXX(GEN z, long n, long v)
1272 : {
1273 1767321 : long i,j,lx,l, N = (n<<1)+1;
1274 : GEN x, t;
1275 1767321 : l = lg(z); lx = (l-2) / (N-2);
1276 1767321 : x = cgetg(lx+3,t_POL);
1277 1767421 : x[1] = z[1];
1278 8928471 : for (i=2; i<lx+2; i++)
1279 : {
1280 7160432 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1281 107165200 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1282 7159563 : z += (N-2);
1283 7159563 : gel(x,i) = ZX_renormalize(t,N);
1284 : }
1285 1768039 : N = (l-2) % (N-2) + 2;
1286 1768039 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1287 5133932 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1288 1767710 : gel(x,i) = ZX_renormalize(t,N);
1289 1767645 : return ZXX_renormalize(x, i+1);
1290 : }
1291 :
1292 : GEN
1293 6676256 : RgXX_to_Kronecker(GEN P, long n)
1294 : {
1295 6676256 : GEN z = RgXX_to_Kronecker_spec(P+2, lgpol(P), n);
1296 6676083 : setvarn(z,varn(P)); return z;
1297 : }
1298 : GEN
1299 3026976 : ZXX_mul_Kronecker(GEN x, GEN y, long n)
1300 3026976 : { return ZX_mul(RgXX_to_Kronecker(x,n), RgXX_to_Kronecker(y,n)); }
1301 :
1302 : GEN
1303 5306 : ZXX_sqr_Kronecker(GEN x, long n)
1304 5306 : { return ZX_sqr(RgXX_to_Kronecker(x,n)); }
1305 :
1306 : /* shallow, n = deg(T) */
1307 : GEN
1308 24269 : Kronecker_to_ZXQX(GEN z, GEN T)
1309 : {
1310 24269 : long i,j,lx,l, N = (degpol(T)<<1)+1;
1311 : GEN x, t;
1312 24269 : l = lg(z); lx = (l-2) / (N-2);
1313 24269 : x = cgetg(lx+3,t_POL);
1314 24269 : x[1] = z[1];
1315 236274 : for (i=2; i<lx+2; i++)
1316 : {
1317 212005 : t = cgetg(N,t_POL); t[1] = T[1];
1318 1898640 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1319 212004 : z += (N-2);
1320 212004 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1321 : }
1322 24269 : N = (l-2) % (N-2) + 2;
1323 24269 : t = cgetg(N,t_POL); t[1] = T[1];
1324 57208 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1325 24269 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1326 24269 : return ZXX_renormalize(x, i+1);
1327 : }
1328 :
1329 : GEN
1330 1904 : ZXQX_sqr(GEN x, GEN T)
1331 : {
1332 1904 : pari_sp av = avma;
1333 1904 : long n = degpol(T);
1334 1904 : GEN z = ZXX_sqr_Kronecker(x, n);
1335 1904 : z = Kronecker_to_ZXQX(z, T);
1336 1904 : return gerepilecopy(av, z);
1337 : }
1338 :
1339 : GEN
1340 22365 : ZXQX_mul(GEN x, GEN y, GEN T)
1341 : {
1342 22365 : pari_sp av = avma;
1343 22365 : long n = degpol(T);
1344 22365 : GEN z = ZXX_mul_Kronecker(x, y, n);
1345 22365 : z = Kronecker_to_ZXQX(z, T);
1346 22365 : return gerepilecopy(av, z);
1347 : }
1348 :
1349 : GEN
1350 8988 : ZXQX_ZXQ_mul(GEN P, GEN U, GEN T)
1351 : {
1352 : long i, lP;
1353 : GEN res;
1354 8988 : res = cgetg_copy(P, &lP); res[1] = P[1];
1355 47106 : for(i=2; i<lP; i++)
1356 38118 : gel(res,i) = typ(gel(P,i))==t_POL? ZXQ_mul(U, gel(P,i), T)
1357 38118 : : gmul(U, gel(P,i));
1358 8988 : return ZXX_renormalize(res,lP);
1359 : }
1360 :
1361 : GEN
1362 3503182 : QX_mul(GEN x, GEN y)
1363 : {
1364 3503182 : GEN dx, nx = Q_primitive_part(x, &dx);
1365 3503182 : GEN dy, ny = Q_primitive_part(y, &dy);
1366 3503182 : GEN z = ZX_mul(nx, ny);
1367 3503182 : if (dx || dy)
1368 : {
1369 3410942 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1370 3410942 : return ZX_Q_mul(z, d);
1371 : } else
1372 92240 : return z;
1373 : }
1374 :
1375 : GEN
1376 92916 : QX_sqr(GEN x)
1377 : {
1378 92916 : GEN dx, nx = Q_primitive_part(x, &dx);
1379 92913 : GEN z = ZX_sqr(nx);
1380 92917 : if (dx)
1381 82538 : return ZX_Q_mul(z, gsqr(dx));
1382 : else
1383 10379 : return z;
1384 : }
1385 :
1386 : GEN
1387 1791401 : QX_ZX_rem(GEN x, GEN y)
1388 : {
1389 1791401 : pari_sp av = avma;
1390 1791401 : GEN d, nx = Q_primitive_part(x, &d);
1391 1791401 : GEN r = ZX_rem(nx, y);
1392 1791401 : if (d) r = ZX_Q_mul(r, d);
1393 1791401 : return gerepileupto(av, r);
1394 : }
1395 :
1396 : GEN
1397 13377 : QXQX_mul(GEN x, GEN y, GEN T)
1398 : {
1399 13377 : GEN dx, nx = Q_primitive_part(x, &dx);
1400 13377 : GEN dy, ny = Q_primitive_part(y, &dy);
1401 13377 : GEN z = ZXQX_mul(nx, ny, T);
1402 13377 : if (dx || dy)
1403 : {
1404 8309 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1405 8309 : return ZXX_Q_mul(z, d);
1406 : } else
1407 5068 : return z;
1408 : }
1409 :
1410 : GEN
1411 1904 : QXQX_sqr(GEN x, GEN T)
1412 : {
1413 1904 : GEN dx, nx = Q_primitive_part(x, &dx);
1414 1904 : GEN z = ZXQX_sqr(nx, T);
1415 1904 : if (dx)
1416 588 : return ZXX_Q_mul(z, gsqr(dx));
1417 : else
1418 1316 : return z;
1419 : }
1420 :
1421 : GEN
1422 476 : QXQX_powers(GEN P, long n, GEN T)
1423 : {
1424 476 : GEN v = cgetg(n+2, t_VEC);
1425 : long i;
1426 476 : gel(v, 1) = pol_1(varn(T));
1427 476 : if (n==0) return v;
1428 476 : gel(v, 2) = gcopy(P);
1429 1554 : for (i = 2; i <= n; i++) gel(v,i+1) = QXQX_mul(P, gel(v,i), T);
1430 476 : return v;
1431 : }
1432 :
1433 : GEN
1434 3465 : QXQX_QXQ_mul(GEN P, GEN U, GEN T)
1435 : {
1436 : long i, lP;
1437 : GEN res;
1438 3465 : res = cgetg_copy(P, &lP); res[1] = P[1];
1439 22715 : for(i=2; i<lP; i++)
1440 19250 : gel(res,i) = typ(gel(P,i))==t_POL? QXQ_mul(U, gel(P,i), T)
1441 19250 : : gmul(U, gel(P,i));
1442 3465 : return ZXX_renormalize(res,lP);
1443 : }
|
