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 324434 : RgX_check_QX(GEN x, const char *s)
25 324434 : { if (!RgX_is_QX(x)) pari_err_TYPE(stack_strcat(s," [not in Q[X]]"), x); }
26 : void
27 1825019 : RgX_check_ZX(GEN x, const char *s)
28 1825019 : { if (!RgX_is_ZX(x)) pari_err_TYPE(stack_strcat(s," [not in Z[X]]"), x); }
29 : long
30 106552 : ZX_max_lg(GEN x)
31 : {
32 106552 : long i, prec = 0, lx = lg(x);
33 :
34 645044 : for (i=2; i<lx; i++) { long l = lgefint(gel(x,i)); if (l > prec) prec = l; }
35 106552 : return prec;
36 : }
37 :
38 : GEN
39 40130878 : ZX_add(GEN x, GEN y)
40 : {
41 : long lx,ly,i;
42 : GEN z;
43 40130878 : lx = lg(x); ly = lg(y); if (lx < ly) swapspec(x,y, lx,ly);
44 40130878 : z = cgetg(lx,t_POL); z[1] = x[1];
45 246308502 : for (i=2; i<ly; i++) gel(z,i) = addii(gel(x,i),gel(y,i));
46 78748902 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
47 39927781 : if (lx == ly) z = ZX_renormalize(z, lx);
48 39943043 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
49 38801552 : return z;
50 : }
51 :
52 : GEN
53 29125485 : ZX_sub(GEN x,GEN y)
54 : {
55 29125485 : long i, lx = lg(x), ly = lg(y);
56 : GEN z;
57 29125485 : if (lx >= ly)
58 : {
59 28679611 : z = cgetg(lx,t_POL); z[1] = x[1];
60 116627006 : for (i=2; i<ly; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
61 28617320 : if (lx == ly)
62 : {
63 23073931 : z = ZX_renormalize(z, lx);
64 23084571 : if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); z = pol_0(varn(x)); }
65 : }
66 : else
67 17134998 : for ( ; i<lx; i++) gel(z,i) = icopy(gel(x,i));
68 : }
69 : else
70 : {
71 445874 : z = cgetg(ly,t_POL); z[1] = y[1];
72 1182272 : for (i=2; i<lx; i++) gel(z,i) = subii(gel(x,i),gel(y,i));
73 1390294 : for ( ; i<ly; i++) gel(z,i) = negi(gel(y,i));
74 : }
75 29130827 : return z;
76 : }
77 :
78 : GEN
79 595303 : ZX_neg(GEN x)
80 : {
81 595303 : long i, l = lg(x);
82 595303 : GEN y = cgetg(l,t_POL);
83 2446239 : y[1] = x[1]; for(i=2; i<l; i++) gel(y,i) = negi(gel(x,i));
84 595366 : return y;
85 : }
86 : GEN
87 5116952 : ZX_copy(GEN x)
88 : {
89 5116952 : long i, l = lg(x);
90 5116952 : GEN y = cgetg(l, t_POL);
91 5116961 : y[1] = x[1];
92 21178825 : for (i=2; i<l; i++)
93 : {
94 16061871 : GEN c = gel(x,i);
95 16061871 : gel(y,i) = lgefint(c) == 2? gen_0: icopy(c);
96 : }
97 5116954 : return y;
98 : }
99 :
100 : GEN
101 56030 : scalar_ZX(GEN x, long v)
102 : {
103 : GEN z;
104 56030 : if (!signe(x)) return pol_0(v);
105 4516 : z = cgetg(3, t_POL);
106 4516 : z[1] = evalsigne(1) | evalvarn(v);
107 4516 : gel(z,2) = icopy(x); return z;
108 : }
109 :
110 : GEN
111 23038 : scalar_ZX_shallow(GEN x, long v)
112 : {
113 : GEN z;
114 23038 : if (!signe(x)) return pol_0(v);
115 19776 : z = cgetg(3, t_POL);
116 19776 : z[1] = evalsigne(1) | evalvarn(v);
117 19776 : gel(z,2) = x; return z;
118 : }
119 :
120 : GEN
121 87833 : ZX_Z_add(GEN y, GEN x)
122 : {
123 : long lz, i;
124 87833 : GEN z = cgetg_copy(y, &lz);
125 87833 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
126 77095 : z[1] = y[1];
127 77095 : gel(z,2) = addii(gel(y,2),x);
128 1856459 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
129 77713 : if (lz==3) z = ZX_renormalize(z,lz);
130 77098 : return z;
131 : }
132 : GEN
133 93707 : ZX_Z_add_shallow(GEN y, GEN x)
134 : {
135 : long lz, i;
136 93707 : GEN z = cgetg_copy(y, &lz);
137 93716 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX_shallow(x,varn(y)); }
138 93527 : z[1] = y[1];
139 93527 : gel(z,2) = addii(gel(y,2),x);
140 548498 : for(i=3; i<lz; i++) gel(z,i) = gel(y,i);
141 93500 : if (lz==3) z = ZX_renormalize(z,lz);
142 93510 : return z;
143 : }
144 :
145 : GEN
146 180753 : ZX_Z_sub(GEN y, GEN x)
147 : {
148 : long lz, i;
149 180753 : GEN z = cgetg_copy(y, &lz);
150 180753 : if (lz == 2)
151 : { /* scalarpol(negi(x), v) */
152 441 : long v = varn(y);
153 441 : set_avma((pari_sp)(z + 2));
154 441 : if (!signe(x)) return pol_0(v);
155 441 : z = cgetg(3,t_POL);
156 441 : z[1] = evalvarn(v) | evalsigne(1);
157 441 : gel(z,2) = negi(x); return z;
158 : }
159 180312 : z[1] = y[1];
160 180312 : gel(z,2) = subii(gel(y,2),x);
161 901606 : for(i=3; i<lz; i++) gel(z,i) = icopy(gel(y,i));
162 180315 : if (lz==3) z = ZX_renormalize(z,lz);
163 180312 : return z;
164 : }
165 :
166 : GEN
167 2688148 : Z_ZX_sub(GEN x, GEN y)
168 : {
169 : long lz, i;
170 2688148 : GEN z = cgetg_copy(y, &lz);
171 2688249 : if (lz == 2) { set_avma((pari_sp)(z + 2)); return scalar_ZX(x,varn(y)); }
172 2688249 : z[1] = y[1];
173 2688249 : gel(z,2) = subii(x, gel(y,2));
174 13473466 : for(i=3; i<lz; i++) gel(z,i) = negi(gel(y,i));
175 2688662 : if (lz==3) z = ZX_renormalize(z,lz);
176 2688328 : return z;
177 : }
178 :
179 : GEN
180 5088873 : ZX_Z_divexact(GEN y,GEN x)
181 : {
182 5088873 : long i, l = lg(y);
183 5088873 : GEN z = cgetg(l,t_POL); z[1] = y[1];
184 32654440 : for(i=2; i<l; i++) gel(z,i) = diviiexact(gel(y,i),x);
185 5085404 : return z;
186 : }
187 :
188 : GEN
189 977795 : ZX_divuexact(GEN y, ulong x)
190 : {
191 977795 : long i, l = lg(y);
192 977795 : GEN z = cgetg(l,t_POL); z[1] = y[1];
193 9684353 : for(i=2; i<l; i++) gel(z,i) = diviuexact(gel(y,i),x);
194 977795 : return z;
195 : }
196 :
197 : GEN
198 112 : zx_z_divexact(GEN y, long x)
199 : {
200 112 : long i, l = lg(y);
201 112 : GEN z = cgetg(l,t_VECSMALL); z[1] = y[1];
202 595 : for (i=2; i<l; i++) z[i] = y[i]/x;
203 112 : return z;
204 : }
205 :
206 : GEN
207 16854736 : ZX_Z_mul(GEN y,GEN x)
208 : {
209 : GEN z;
210 : long i, l;
211 16854736 : if (!signe(x)) return pol_0(varn(y));
212 13811806 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
213 92346814 : for(i=2; i<l; i++) gel(z,i) = mulii(gel(y,i),x);
214 13802531 : return z;
215 : }
216 :
217 : GEN
218 252619 : ZX_mulu(GEN y, ulong x)
219 : {
220 : GEN z;
221 : long i, l;
222 252619 : if (!x) return pol_0(varn(y));
223 252619 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
224 1024336 : for(i=2; i<l; i++) gel(z,i) = mului(x,gel(y,i));
225 252619 : return z;
226 : }
227 :
228 : GEN
229 6324120 : ZX_shifti(GEN y, long n)
230 : {
231 : GEN z;
232 : long i, l;
233 6324120 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
234 30594841 : for(i=2; i<l; i++) gel(z,i) = shifti(gel(y,i),n);
235 6304676 : return ZX_renormalize(z,l);
236 : }
237 :
238 : GEN
239 107190 : ZX_remi2n(GEN y, long n)
240 : {
241 : GEN z;
242 : long i, l;
243 107190 : l = lg(y); z = cgetg(l,t_POL); z[1] = y[1];
244 4734020 : for(i=2; i<l; i++) gel(z,i) = remi2n(gel(y,i),n);
245 100537 : return ZX_renormalize(z,l);
246 : }
247 :
248 : GEN
249 46383 : ZXT_remi2n(GEN z, long n)
250 : {
251 46383 : if (typ(z) == t_POL)
252 39992 : return ZX_remi2n(z, n);
253 : else
254 : {
255 6391 : long i,l = lg(z);
256 6391 : GEN x = cgetg(l, t_VEC);
257 19181 : for (i=1; i<l; i++) gel(x,i) = ZXT_remi2n(gel(z,i), n);
258 6390 : return x;
259 : }
260 : }
261 :
262 : GEN
263 1323 : zx_to_ZX(GEN z)
264 : {
265 1323 : long i, l = lg(z);
266 1323 : GEN x = cgetg(l,t_POL);
267 100912 : for (i=2; i<l; i++) gel(x,i) = stoi(z[i]);
268 1323 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
269 : }
270 :
271 : GEN
272 5393999 : ZX_deriv(GEN x)
273 : {
274 5393999 : long i,lx = lg(x)-1;
275 : GEN y;
276 :
277 5393999 : if (lx<3) return pol_0(varn(x));
278 5311779 : y = cgetg(lx,t_POL);
279 30627392 : for (i=2; i<lx ; i++) gel(y,i) = mului(i-1,gel(x,i+1));
280 5310025 : y[1] = x[1]; return y;
281 : }
282 :
283 : int
284 2053071 : ZX_equal(GEN V, GEN W)
285 : {
286 2053071 : long i, l = lg(V);
287 2053071 : if (lg(W) != l) return 0;
288 5988388 : for (i = 2; i < l; i++)
289 4652768 : if (!equalii(gel(V,i), gel(W,i))) return 0;
290 1335620 : return 1;
291 : }
292 :
293 : static long
294 300638543 : ZX_valspec(GEN x, long nx)
295 : {
296 : long vx;
297 372829465 : for (vx = 0; vx<nx ; vx++)
298 372844290 : if (signe(gel(x,vx))) break;
299 300638543 : return vx;
300 : }
301 :
302 : long
303 328829 : ZX_val(GEN x)
304 : {
305 : long vx;
306 328829 : if (!signe(x)) return LONG_MAX;
307 263351 : for (vx = 0;; vx++)
308 263351 : if (signe(gel(x,2+vx))) break;
309 213509 : return vx;
310 : }
311 : long
312 24066046 : ZX_valrem(GEN x, GEN *Z)
313 : {
314 : long vx;
315 24066046 : if (!signe(x)) { *Z = pol_0(varn(x)); return LONG_MAX; }
316 44758415 : for (vx = 0;; vx++)
317 44758415 : if (signe(gel(x,2+vx))) break;
318 24066046 : *Z = RgX_shift_shallow(x, -vx);
319 24066066 : return vx;
320 : }
321 :
322 : GEN
323 21 : ZX_div_by_X_1(GEN a, GEN *r)
324 : {
325 21 : long l = lg(a), i;
326 21 : GEN a0, z0, z = cgetg(l-1, t_POL);
327 21 : z[1] = a[1];
328 21 : a0 = a + l-1;
329 21 : z0 = z + l-2; *z0 = *a0--;
330 4732 : for (i=l-3; i>1; i--) /* z[i] = a[i+1] + z[i+1] */
331 : {
332 4711 : GEN t = addii(gel(a0--,0), gel(z0--,0));
333 4711 : gel(z0,0) = t;
334 : }
335 21 : if (r) *r = addii(gel(a0,0), gel(z0,0));
336 21 : return z;
337 : }
338 :
339 : /* return P(X + c) using destructive Horner, optimize for c = 1,-1 */
340 : static GEN
341 1278450 : ZX_translate_basecase(GEN P, GEN c)
342 : {
343 1278450 : pari_sp av = avma;
344 : GEN Q, R;
345 : long i, k, n;
346 :
347 1278450 : if (!signe(P) || !signe(c)) return ZX_copy(P);
348 1237582 : Q = leafcopy(P);
349 1237599 : R = Q+2; n = degpol(P);
350 1237600 : if (equali1(c))
351 : {
352 5960005 : for (i=1; i<=n; i++)
353 : {
354 35242216 : for (k=n-i; k<n; k++) gel(R,k) = addii(gel(R,k), gel(R,k+1));
355 4992111 : if (gc_needed(av,2))
356 : {
357 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(1), i = %ld/%ld", i,n);
358 0 : Q = gerepilecopy(av, Q); R = Q+2;
359 : }
360 : }
361 : }
362 269702 : else if (equalim1(c))
363 : {
364 55555 : for (i=1; i<=n; i++)
365 : {
366 182238 : for (k=n-i; k<n; k++) gel(R,k) = subii(gel(R,k), gel(R,k+1));
367 45932 : if (gc_needed(av,2))
368 : {
369 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate(-1), i = %ld/%ld", i,n);
370 0 : Q = gerepilecopy(av, Q); R = Q+2;
371 : }
372 : }
373 : }
374 : else
375 : {
376 1607070 : for (i=1; i<=n; i++)
377 : {
378 7600542 : for (k=n-i; k<n; k++) gel(R,k) = addmulii_inplace(gel(R,k), c, gel(R,k+1));
379 1346991 : if (gc_needed(av,2))
380 : {
381 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_translate, i = %ld/%ld", i,n);
382 0 : Q = gerepilecopy(av, Q); R = Q+2;
383 : }
384 : }
385 : }
386 1237379 : return gerepilecopy(av, Q);
387 : }
388 :
389 : static GEN
390 0 : Z_Xpm1_powu(long n, long s, long v)
391 : {
392 : long d, k;
393 : GEN C;
394 0 : if (!n) return pol_1(v);
395 0 : d = (n + 1) >> 1;
396 0 : C = cgetg(n+3, t_POL);
397 0 : C[1] = evalsigne(1)| evalvarn(v);
398 0 : gel(C,2) = gen_1;
399 0 : gel(C,3) = utoipos(n);
400 0 : for (k=2; k <= d; k++)
401 0 : gel(C,k+2) = diviuexact(mului(n-k+1, gel(C,k+1)), k);
402 0 : if (s < 0)
403 0 : for (k = odd(n)? 0: 1; k <= d; k += 2)
404 0 : togglesign_safe(&gel(C,k+2));
405 0 : if (s > 0 || !odd(n))
406 0 : for (k = d+1; k <= n; k++) gel(C,k+2) = gel(C,n-k+2);
407 : else
408 0 : for (k = d+1; k <= n; k++) gel(C,k+2) = negi(gel(C,n-k+2));
409 0 : return C;
410 : }
411 : /* return (x+u)^n */
412 : static GEN
413 0 : Z_XpN_powu(GEN u, long n, long v)
414 : {
415 : pari_sp av;
416 : long k;
417 : GEN B, C, V;
418 0 : if (!n) return pol_1(v);
419 0 : if (is_pm1(u))
420 0 : return Z_Xpm1_powu(n, signe(u), v);
421 0 : av = avma;
422 0 : V = gpowers(u, n);
423 0 : B = vecbinomial(n);
424 0 : C = cgetg(n+3, t_POL);
425 0 : C[1] = evalsigne(1)| evalvarn(v);
426 0 : for (k=1; k <= n+1; k++)
427 0 : gel(C,k+1) = mulii(gel(V,n+2-k), gel(B,k));
428 0 : return gerepileupto(av, C);
429 : }
430 :
431 : GEN
432 1278446 : ZX_translate(GEN P, GEN c)
433 : {
434 1278446 : pari_sp av = avma;
435 1278446 : long n = degpol(P);
436 1278451 : if (n < 220)
437 1278451 : return ZX_translate_basecase(P, c);
438 : else
439 : {
440 0 : long d = n >> 1;
441 0 : GEN Q = ZX_translate(RgX_shift_shallow(P, -d), c);
442 0 : GEN R = ZX_translate(RgXn_red_shallow(P, d), c);
443 0 : GEN S = Z_XpN_powu(c, d, varn(P));
444 0 : return gerepileupto(av, ZX_add(ZX_mul(Q, S), R));
445 : }
446 : }
447 :
448 : /* P(ax + b) */
449 : GEN
450 981188 : ZX_affine(GEN P, GEN a, GEN b)
451 : {
452 981188 : if (signe(b)) P = ZX_translate(P, b);
453 981188 : return ZX_unscale(P, a);
454 : }
455 :
456 : GEN
457 506979 : ZX_Z_eval(GEN x, GEN y)
458 : {
459 506979 : long i = lg(x)-1, j;
460 506979 : pari_sp av = avma;
461 : GEN t, r;
462 :
463 506979 : if (i<=2) return (i==2)? icopy(gel(x,2)): gen_0;
464 469522 : if (!signe(y)) return icopy(gel(x,2));
465 :
466 469403 : t = gel(x,i); i--;
467 : #if 0 /* standard Horner's rule */
468 : for ( ; i>=2; i--)
469 : t = addii(mulii(t,y),gel(x,i));
470 : #endif
471 : /* specific attention to sparse polynomials */
472 2367790 : for ( ; i>=2; i = j-1)
473 : {
474 2074528 : for (j = i; !signe(gel(x,j)); j--)
475 176133 : if (j==2)
476 : {
477 2947 : if (i != j) y = powiu(y, i-j+1);
478 2947 : return gerepileuptoint(av, mulii(t,y));
479 : }
480 1898395 : r = (i==j)? y: powiu(y, i-j+1);
481 1898399 : t = addii(mulii(t,r), gel(x,j));
482 1898382 : if (gc_needed(av,2))
483 : {
484 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_Z_eval: i = %ld",i);
485 0 : t = gerepileuptoint(av, t);
486 : }
487 : }
488 466448 : return gerepileuptoint(av, t);
489 : }
490 :
491 : /* Return 2^(n degpol(P)) P(x >> n) */
492 : GEN
493 0 : ZX_rescale2n(GEN P, long n)
494 : {
495 0 : long i, l = lg(P), ni = n;
496 : GEN Q;
497 0 : if (l==2) return pol_0(varn(P));
498 0 : Q = cgetg(l,t_POL);
499 0 : gel(Q,l-1) = icopy(gel(P,l-1));
500 0 : for (i=l-2; i>=2; i--)
501 : {
502 0 : gel(Q,i) = shifti(gel(P,i), ni);
503 0 : ni += n;
504 : }
505 0 : Q[1] = P[1]; return Q;
506 : }
507 :
508 : /* Return h^deg(P) P(x / h), not memory clean. h integer, P ZX */
509 : GEN
510 35614 : ZX_rescale(GEN P, GEN h)
511 : {
512 35614 : long l = lg(P);
513 35614 : GEN Q = cgetg(l,t_POL);
514 35615 : if (l != 2)
515 : {
516 35615 : long i = l-1;
517 35615 : GEN hi = h;
518 35615 : gel(Q,i) = gel(P,i);
519 35615 : if (l != 3) { i--; gel(Q,i) = mulii(gel(P,i), h); }
520 110225 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
521 : }
522 35614 : Q[1] = P[1]; return Q;
523 : }
524 : /* Return h^(deg(P)-1) P(x / h), P!=0, h=lt(P), memory unclean; monic result */
525 : GEN
526 0 : ZX_rescale_lt(GEN P)
527 : {
528 0 : long l = lg(P);
529 0 : GEN Q = cgetg(l,t_POL);
530 0 : gel(Q,l-1) = gen_1;
531 0 : if (l != 3)
532 : {
533 0 : long i = l-1;
534 0 : GEN h = gel(P,i), hi = h;
535 0 : i--; gel(Q,i) = gel(P,i);
536 0 : if (l != 4) { i--; gel(Q,i) = mulii(gel(P,i), h); }
537 0 : for (i--; i>=2; i--) { hi = mulii(hi,h); gel(Q,i) = mulii(gel(P,i), hi); }
538 : }
539 0 : Q[1] = P[1]; return Q;
540 : }
541 :
542 : /*Eval x in 2^(k*BIL) in linear time*/
543 : static GEN
544 88604033 : ZX_eval2BILspec(GEN x, long k, long nx)
545 : {
546 88604033 : pari_sp av = avma;
547 88604033 : long i,j, lz = k*nx, ki;
548 88604033 : GEN pz = cgetipos(2+lz);
549 88601911 : GEN nz = cgetipos(2+lz);
550 3139177507 : for(i=0; i < lz; i++)
551 : {
552 3050557372 : *int_W(pz,i) = 0UL;
553 3050557372 : *int_W(nz,i) = 0UL;
554 : }
555 981298503 : for(i=0, ki=0; i<nx; i++, ki+=k)
556 : {
557 892678368 : GEN c = gel(x,i);
558 892678368 : long lc = lgefint(c)-2;
559 892678368 : if (signe(c)==0) continue;
560 668965185 : if (signe(c) > 0)
561 1675722031 : for (j=0; j<lc; j++) *int_W(pz,j+ki) = *int_W(c,j);
562 : else
563 301762107 : for (j=0; j<lc; j++) *int_W(nz,j+ki) = *int_W(c,j);
564 : }
565 88620135 : pz = int_normalize(pz,0);
566 88621579 : nz = int_normalize(nz,0); return gerepileuptoint(av, subii(pz,nz));
567 : }
568 :
569 : static long
570 162392812 : ZX_expispec(GEN x, long nx)
571 : {
572 162392812 : long i, m = 0;
573 1245418817 : for(i = 0; i < nx; i++)
574 : {
575 1083082060 : long e = expi(gel(x,i));
576 1083026005 : if (e > m) m = e;
577 : }
578 162336757 : return m;
579 : }
580 :
581 : static GEN
582 53580740 : Z_mod2BIL_ZX(GEN x, long bs, long d, long vx)
583 : {
584 53580740 : long i, offset, lm = lgefint(x)-2, l = d+vx+3, sx = signe(x);
585 53580740 : GEN s1 = int2n(bs*BITS_IN_LONG), pol = cgetg(l, t_POL);
586 53568596 : int carry = 0;
587 53568596 : pol[1] = evalsigne(1);
588 74073778 : for (i=0; i<vx; i++) gel(pol,i+2) = gen_0;
589 979572199 : for (offset=0; i <= d+vx; i++, offset += bs)
590 : {
591 926388033 : pari_sp av = avma;
592 926388033 : long lz = minss(bs, lm-offset);
593 927386763 : GEN z = lz > 0 ?adduispec_offset(carry, x, offset, lz): utoi(carry);
594 925686269 : if (lgefint(z) == 3+bs) { carry = 1; z = gen_0;}
595 : else
596 : {
597 912264436 : carry = (lgefint(z) == 2+bs && (HIGHBIT & *int_W(z,bs-1)));
598 912264436 : if (carry)
599 53766895 : z = gerepileuptoint(av, (sx==-1)? subii(s1,z): subii(z,s1));
600 858497541 : else if (sx==-1) togglesign(z);
601 : }
602 926003603 : gel(pol,i+2) = z;
603 : }
604 53184166 : return ZX_renormalize(pol,l);
605 : }
606 :
607 : static GEN
608 17331726 : ZX_sqrspec_sqri(GEN x, long nx, long ex, long v)
609 : {
610 17331726 : long e = 2*ex + expu(nx) + 3;
611 17331518 : long N = divsBIL(e)+1;
612 17331462 : GEN z = sqri(ZX_eval2BILspec(x,N,nx));
613 17330684 : return Z_mod2BIL_ZX(z, N, nx*2-2, v);
614 : }
615 :
616 : static GEN
617 35191388 : ZX_mulspec_mulii(GEN x, GEN y, long nx, long ny, long ex, long ey, long v)
618 : {
619 35191388 : long e = ex + ey + expu(minss(nx,ny)) + 3;
620 35190683 : long N = divsBIL(e)+1;
621 35191021 : GEN z = mulii(ZX_eval2BILspec(x,N,nx), ZX_eval2BILspec(y,N,ny));
622 35191094 : return Z_mod2BIL_ZX(z, N, nx+ny-2, v);
623 : }
624 :
625 : INLINE GEN
626 307307908 : ZX_sqrspec_basecase_limb(GEN x, long a, long i)
627 : {
628 307307908 : pari_sp av = avma;
629 307307908 : GEN s = gen_0;
630 307307908 : long j, l = (i+1)>>1;
631 611111182 : for (j=a; j<l; j++)
632 : {
633 304550499 : GEN xj = gel(x,j), xx = gel(x,i-j);
634 304550499 : if (signe(xj) && signe(xx))
635 300454215 : s = addii(s, mulii(xj, xx));
636 : }
637 306560683 : s = shifti(s,1);
638 306684304 : if ((i&1) == 0)
639 : {
640 190580671 : GEN t = gel(x, i>>1);
641 190580671 : if (signe(t))
642 190287560 : s = addii(s, sqri(t));
643 : }
644 306254996 : return gerepileuptoint(av,s);
645 : }
646 :
647 : static GEN
648 73770048 : ZX_sqrspec_basecase(GEN x, long nx, long v)
649 : {
650 : long i, lz, nz;
651 : GEN z;
652 :
653 73770048 : lz = (nx << 1) + 1; nz = lz-2;
654 73770048 : lz += v;
655 73770048 : z = cgetg(lz,t_POL); z[1] = evalsigne(1); z += 2;
656 73970274 : for (i=0; i<v; i++) gel(z++, 0) = gen_0;
657 264345814 : for (i=0; i<nx; i++)
658 190654739 : gel(z,i) = ZX_sqrspec_basecase_limb(x, 0, i);
659 190655829 : for ( ; i<nz; i++) gel(z,i) = ZX_sqrspec_basecase_limb(x, i-nx+1, i);
660 73694638 : z -= v+2; return z;
661 : }
662 :
663 : static GEN
664 41996582 : Z_sqrshiftspec_ZX(GEN x, long vx)
665 : {
666 41996582 : long i, nz = 2*vx+1;
667 41996582 : GEN z = cgetg(2+nz, t_POL);
668 41991136 : z[1] = evalsigne(1);
669 44549721 : for(i=2;i<nz+1;i++) gel(z,i) = gen_0;
670 41991136 : gel(z,nz+1) = sqri(x);
671 41935759 : return z;
672 : }
673 :
674 : static GEN
675 36656018 : Z_ZX_mulshiftspec(GEN x, GEN y, long ny, long vz)
676 : {
677 36656018 : long i, nz = vz+ny;
678 36656018 : GEN z = cgetg(2+nz, t_POL);
679 36656003 : z[1] = evalsigne(1);
680 85054924 : for (i=0; i<vz; i++) gel(z,i+2) = gen_0;
681 113066792 : for (i=0; i<ny; i++) gel(z,i+vz+2) = mulii(x, gel(y,i));
682 36645442 : return z;
683 : }
684 :
685 : GEN
686 133742171 : ZX_sqrspec(GEN x, long nx)
687 : {
688 : #ifdef PARI_KERNEL_GMP
689 76679876 : 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};
690 76679876 : 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};
691 : #else
692 57062295 : 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};
693 57062295 : 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};
694 : #endif
695 133742171 : const long nblow = numberof(low);
696 : pari_sp av;
697 : long ex, vx;
698 : GEN z;
699 133742171 : if (!nx) return pol_0(0);
700 132982705 : vx = ZX_valspec(x,nx); nx-=vx; x+=vx;
701 133104991 : if (nx==1) return Z_sqrshiftspec_ZX(gel(x, 0), vx);
702 91106786 : av = avma;
703 91106786 : ex = ZX_expispec(x,nx);
704 91091713 : if (nx-2 < nblow && low[nx-2]<=ex && ex<=high[nx-2])
705 73774376 : z = ZX_sqrspec_basecase(x, nx, 2*vx);
706 : else
707 17317337 : z = ZX_sqrspec_sqri(x, nx, ex, 2*vx);
708 91059992 : return gerepileupto(av, z);
709 : }
710 :
711 : GEN
712 133752961 : ZX_sqr(GEN x)
713 : {
714 133752961 : GEN z = ZX_sqrspec(x+2, lgpol(x));
715 133614741 : z[1] = x[1];
716 133614741 : return z;
717 : }
718 :
719 : GEN
720 93579932 : ZX_mulspec(GEN x, GEN y, long nx, long ny)
721 : {
722 : pari_sp av;
723 : long ex, ey, vx, vy, v;
724 93579932 : if (!nx || !ny) return pol_0(0);
725 83774224 : vx = ZX_valspec(x,nx); nx-=vx; x += vx;
726 83774132 : vy = ZX_valspec(y,ny); ny-=vy; y += vy;
727 83777344 : v = vx + vy;
728 83777344 : if (nx==1) return Z_ZX_mulshiftspec(gel(x,0), y, ny, v);
729 56937802 : if (ny==1) return Z_ZX_mulshiftspec(gel(y,0), x, nx, v);
730 47121146 : if (nx == 2 && ny == 2)
731 : {
732 11929988 : GEN a0 = gel(x,0), a1 = gel(x,1), A0, A1, A2;
733 11929988 : GEN b0 = gel(y,0), b1 = gel(y,1), z = cgetg(5 + v, t_POL);
734 : long i;
735 11930026 : z[1] = evalvarn(0) | evalsigne(1);
736 11930026 : A0 = mulii(a0, b0);
737 11929677 : A2 = mulii(a1, b1); av = avma;
738 11929595 : A1 = gerepileuptoint(av, subii(addii(A0,A2),
739 : mulii(subii(a1,a0), subii(b1,b0))));
740 11929887 : i = 4 + v;
741 11929887 : gel(z,i--) = A2;
742 11929887 : gel(z,i--) = A1;
743 14006996 : gel(z,i--) = A0; while (i > 1) gel(z, i--) = gen_0;
744 11929887 : return z;
745 : }
746 : #if 0
747 : /* generically slower even when degrees differ a lot; sometimes about twice
748 : * faster when bitsize is moderate */
749 : if (DEBUGVAR)
750 : return RgX_mulspec(x - vx, y - vy, nx + vx, ny + vy);
751 : #endif
752 35191158 : av = avma;
753 35191158 : ex = ZX_expispec(x, nx);
754 35190012 : ey = ZX_expispec(y, ny);
755 35190543 : return gerepileupto(av, ZX_mulspec_mulii(x,y,nx,ny,ex,ey,v));
756 : }
757 : GEN
758 86317968 : ZX_mul(GEN x, GEN y)
759 : {
760 : GEN z;
761 86317968 : if (x == y) return ZX_sqr(x);
762 85582037 : z = ZX_mulspec(x+2,y+2,lgpol(x),lgpol(y));
763 85573386 : z[1] = x[1];
764 85573386 : if (!signe(y)) z[1] &= VARNBITS;
765 85573386 : return z;
766 : }
767 :
768 : /* x,y two ZX in the same variable; assume y is monic */
769 : GEN
770 10769032 : ZX_rem(GEN x, GEN y)
771 : {
772 : long vx, dx, dy, dz, i, j, sx, lr;
773 : pari_sp av0, av;
774 : GEN z,p1,rem;
775 :
776 10769032 : vx = varn(x);
777 10769032 : dy = degpol(y);
778 10768970 : dx = degpol(x);
779 10768967 : if (dx < dy) return ZX_copy(x);
780 6009508 : if (!dy) return pol_0(vx); /* y is constant */
781 6006813 : av0 = avma; dz = dx-dy;
782 6006813 : z=cgetg(dz+3,t_POL); z[1] = x[1];
783 6006853 : x += 2; y += 2; z += 2;
784 :
785 6006853 : p1 = gel(x,dx);
786 6006853 : gel(z,dz) = icopy(p1);
787 28556934 : for (i=dx-1; i>=dy; i--)
788 : {
789 22550019 : av=avma; p1=gel(x,i);
790 212353867 : for (j=i-dy+1; j<=i && j<=dz; j++)
791 189817325 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
792 22536542 : gel(z,i-dy) = avma == av? icopy(p1): gerepileuptoint(av, p1);
793 : }
794 6006915 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
795 8860289 : for (sx=0; ; i--)
796 : {
797 8860289 : p1 = gel(x,i);
798 57195914 : for (j=0; j<=i && j<=dz; j++)
799 48335901 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
800 8860013 : if (signe(p1)) { sx = 1; break; }
801 3226040 : if (!i) break;
802 2853476 : set_avma(av);
803 : }
804 6006537 : lr=i+3; rem -= lr;
805 6006537 : rem[0] = evaltyp(t_POL) | evallg(lr);
806 6006502 : rem[1] = z[-1];
807 6006502 : p1 = gerepileuptoint((pari_sp)rem, p1);
808 6004117 : rem += 2; gel(rem,i) = p1;
809 33684620 : for (i--; i>=0; i--)
810 : {
811 27677878 : av=avma; p1 = gel(x,i);
812 248990561 : for (j=0; j<=i && j<=dz; j++)
813 221351127 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
814 27639434 : gel(rem,i) = avma == av? icopy(p1): gerepileuptoint(av, p1);
815 : }
816 6006742 : rem -= 2;
817 6006742 : if (!sx) (void)ZX_renormalize(rem, lr);
818 6006742 : return gerepileupto(av0,rem);
819 : }
820 :
821 : /* return x(1) */
822 : GEN
823 4718 : ZX_eval1(GEN x)
824 : {
825 4718 : pari_sp av = avma;
826 4718 : long i = lg(x)-1;
827 : GEN s;
828 4718 : if (i < 2) return gen_0;
829 4193 : s = gel(x,i); i--;
830 4193 : if (i == 1) return icopy(s);
831 34608 : for ( ; i>=2; i--)
832 : {
833 31423 : GEN c = gel(x,i);
834 31423 : if (signe(c)) s = addii(s, c);
835 : }
836 3185 : return gerepileuptoint(av,s);
837 : }
838 :
839 : /* reduce T mod X^n - 1. Shallow function */
840 : GEN
841 1790573 : ZX_mod_Xnm1(GEN T, ulong n)
842 : {
843 1790573 : long i, j, L = lg(T), l = n+2;
844 : GEN S;
845 1790573 : if (L <= l) return T;
846 1599556 : S = cgetg(l, t_POL);
847 1599081 : S[1] = T[1];
848 12293896 : for (i = 2; i < l; i++) gel(S,i) = gel(T,i);
849 7465914 : for (j = 2; i < L; i++) {
850 5873418 : gel(S,j) = addii(gel(S,j), gel(T,i));
851 5866833 : if (++j == l) j = 2;
852 : }
853 1592496 : return normalizepol_lg(S, l);
854 : }
855 :
856 : static GEN
857 0 : _ZX_mul(void* E, GEN x, GEN y)
858 0 : { (void) E; return ZX_mul(x, y); }
859 : static GEN
860 0 : _ZX_sqr(void *E, GEN x)
861 0 : { (void) E; return ZX_sqr(x); }
862 :
863 : static GEN
864 0 : _ZX_divrem(void * E, GEN x, GEN y, GEN *r)
865 0 : { (void) E; return RgX_divrem(x, y, r); }
866 :
867 : static GEN
868 0 : _ZX_add(void * E, GEN x, GEN y)
869 0 : { (void) E; return ZX_add(x, y); }
870 :
871 : static struct bb_ring ZX_ring = { _ZX_add,_ZX_mul,_ZX_sqr };
872 :
873 : GEN
874 0 : ZX_digits(GEN x, GEN T)
875 : {
876 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
877 0 : return gen_digits(x, T, n, NULL, &ZX_ring, _ZX_divrem);
878 : }
879 :
880 : GEN
881 0 : ZXV_ZX_fromdigits(GEN x, GEN T)
882 0 : { return gen_fromdigits(x,T, NULL, &ZX_ring); }
883 :
884 : /*******************************************************************/
885 : /* */
886 : /* ZXV */
887 : /* */
888 : /*******************************************************************/
889 :
890 : int
891 98 : ZXV_equal(GEN V, GEN W)
892 : {
893 98 : long l = lg(V);
894 98 : if (l!=lg(W)) return 0;
895 98 : while (--l > 0)
896 98 : if (!ZX_equal(gel(V,l), gel(W,l))) return 0;
897 0 : return 1;
898 : }
899 :
900 : GEN
901 7476 : ZXV_Z_mul(GEN x, GEN y)
902 29904 : { pari_APPLY_same(ZX_Z_mul(gel(x,i), y)) }
903 :
904 : GEN
905 0 : ZXV_remi2n(GEN x, long N)
906 0 : { pari_APPLY_same(ZX_remi2n(gel(x,i), N)) }
907 :
908 : GEN
909 82943 : ZXV_dotproduct(GEN x, GEN y)
910 : {
911 82943 : pari_sp av = avma;
912 82943 : long i, lx = lg(x);
913 : GEN c;
914 82943 : if (lx == 1) return pol_0(varn(x));
915 82943 : c = ZX_mul(gel(x,1), gel(y,1));
916 392175 : for (i = 2; i < lx; i++)
917 : {
918 309232 : GEN t = ZX_mul(gel(x,i), gel(y,i));
919 309232 : if (signe(t)) c = ZX_add(c, t);
920 : }
921 82943 : return gerepileupto(av, c);
922 : }
923 :
924 : GEN
925 122292 : ZXC_to_FlxC(GEN x, ulong p, long sv)
926 5044782 : { pari_APPLY_type(t_COL,typ(gel(x,i))==t_INT ?
927 : Z_to_Flx(gel(x,i), p, sv): ZX_to_Flx(gel(x,i), p)) }
928 :
929 : GEN
930 11579 : ZXM_to_FlxM(GEN x, ulong p, long sv)
931 102097 : { pari_APPLY_same(ZXC_to_FlxC(gel(x,i), p, sv)) }
932 :
933 : /*******************************************************************/
934 : /* */
935 : /* ZXQM */
936 : /* */
937 : /*******************************************************************/
938 :
939 : GEN
940 2598183 : ZXn_mul(GEN x, GEN y, long n)
941 2598183 : { return RgXn_red_shallow(ZX_mul(x, y), n); }
942 :
943 : GEN
944 1407 : ZXn_sqr(GEN x, long n)
945 1407 : { return RgXn_red_shallow(ZX_sqr(x), n); }
946 :
947 : /*******************************************************************/
948 : /* */
949 : /* ZXQM */
950 : /* */
951 : /*******************************************************************/
952 :
953 : static long
954 3065436 : ZX_expi(GEN x)
955 : {
956 3065436 : if (signe(x)==0) return 0;
957 1385107 : if (typ(x)==t_INT) return expi(x);
958 907245 : return ZX_expispec(x+2, lgpol(x));
959 : }
960 :
961 : static long
962 713417 : ZXC_expi(GEN x)
963 : {
964 713417 : long i, l = lg(x), m=0;
965 3778853 : for(i = 1; i < l; i++)
966 : {
967 3065436 : long e = ZX_expi(gel(x,i));
968 3065436 : if (e > m) m = e;
969 : }
970 713417 : return m;
971 : }
972 :
973 : static long
974 137354 : ZXM_expi(GEN x)
975 : {
976 137354 : long i, l = lg(x), m=0;
977 850771 : for(i = 1; i < l; i++)
978 : {
979 713417 : long e = ZXC_expi(gel(x,i));
980 713417 : if (e > m) m = e;
981 : }
982 137354 : return m;
983 : }
984 :
985 : static GEN
986 3065436 : ZX_eval2BIL(GEN x, long k)
987 : {
988 3065436 : if (signe(x)==0) return gen_0;
989 1385107 : if (typ(x)==t_INT) return x;
990 907245 : return ZX_eval2BILspec(x+2, k, lgpol(x));
991 : }
992 :
993 : /*Eval x in 2^(k*BIL) in linear time*/
994 : static GEN
995 713417 : ZXC_eval2BIL(GEN x, long k)
996 : {
997 713417 : long i, lx = lg(x);
998 713417 : GEN A = cgetg(lx, t_COL);
999 3778853 : for (i=1; i<lx; i++) gel(A,i) = ZX_eval2BIL(gel(x,i), k);
1000 713417 : return A;
1001 : }
1002 :
1003 : static GEN
1004 137354 : ZXM_eval2BIL(GEN x, long k)
1005 : {
1006 137354 : long i, lx = lg(x);
1007 137354 : GEN A = cgetg(lx, t_MAT);
1008 850771 : for (i=1; i<lx; i++) gel(A,i) = ZXC_eval2BIL(gel(x,i), k);
1009 137354 : return A;
1010 : }
1011 :
1012 : static GEN
1013 61292 : Z_mod2BIL_ZXQ(GEN x, long bs, GEN T)
1014 : {
1015 61292 : pari_sp av = avma;
1016 61292 : long v = varn(T), d = 2*(degpol(T)-1);
1017 61292 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1018 61292 : setvarn(z, v);
1019 61292 : return gerepileupto(av, ZX_rem(z, T));
1020 : }
1021 :
1022 : static GEN
1023 12768 : ZC_mod2BIL_ZXQC(GEN x, long bs, GEN T)
1024 : {
1025 12768 : long i, lx = lg(x);
1026 12768 : GEN A = cgetg(lx, t_COL);
1027 74060 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_ZXQ(gel(x,i), bs, T);
1028 12768 : return A;
1029 : }
1030 :
1031 : static GEN
1032 5957 : ZM_mod2BIL_ZXQM(GEN x, long bs, GEN T)
1033 : {
1034 5957 : long i, lx = lg(x);
1035 5957 : GEN A = cgetg(lx, t_MAT);
1036 18725 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_ZXQC(gel(x,i), bs, T);
1037 5957 : return A;
1038 : }
1039 :
1040 : GEN
1041 5817 : ZXQM_mul(GEN x, GEN y, GEN T)
1042 : {
1043 5817 : long d = degpol(T);
1044 : GEN z;
1045 5817 : pari_sp av = avma;
1046 5817 : if (d == 0)
1047 0 : z = ZM_mul(simplify_shallow(x),simplify_shallow(y));
1048 : else
1049 : {
1050 5817 : long e, N, ex = ZXM_expi(x), ey = ZXM_expi(y), n = lg(x)-1;
1051 5817 : e = ex + ey + expu(d) + expu(n) + 4;
1052 5817 : N = divsBIL(e)+1;
1053 5817 : z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1054 5817 : z = ZM_mod2BIL_ZXQM(z, N, T);
1055 : }
1056 5817 : return gerepileupto(av, z);
1057 : }
1058 :
1059 : GEN
1060 140 : ZXQM_sqr(GEN x, GEN T)
1061 : {
1062 140 : long d = degpol(T);
1063 : GEN z;
1064 140 : pari_sp av = avma;
1065 140 : if (d == 0)
1066 0 : z = ZM_sqr(simplify_shallow(x));
1067 : else
1068 : {
1069 140 : long ex = ZXM_expi(x), d = degpol(T), n = lg(x)-1;
1070 140 : long e = 2*ex + expu(d) + expu(n) + 4;
1071 140 : long N = divsBIL(e)+1;
1072 140 : z = ZM_sqr(ZXM_eval2BIL(x,N));
1073 140 : z = ZM_mod2BIL_ZXQM(z, N, T);
1074 : }
1075 140 : return gerepileupto(av, z);
1076 : }
1077 :
1078 : GEN
1079 4396 : QXQM_mul(GEN x, GEN y, GEN T)
1080 : {
1081 4396 : GEN dx, nx = Q_primitive_part(x, &dx);
1082 4396 : GEN dy, ny = Q_primitive_part(y, &dy);
1083 4396 : GEN z = ZXQM_mul(nx, ny, T);
1084 4396 : if (dx || dy)
1085 : {
1086 4396 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1087 4396 : if (!gequal1(d)) z = RgM_Rg_mul(z, d);
1088 : }
1089 4396 : return z;
1090 : }
1091 :
1092 : GEN
1093 7 : QXQM_sqr(GEN x, GEN T)
1094 : {
1095 7 : GEN dx, nx = Q_primitive_part(x, &dx);
1096 7 : GEN z = ZXQM_sqr(nx, T);
1097 7 : if (dx) z = RgM_Rg_mul(z, gsqr(dx));
1098 7 : return z;
1099 : }
1100 :
1101 : static GEN
1102 1002608 : Z_mod2BIL_Fq(GEN x, long bs, GEN T, GEN p)
1103 : {
1104 1002608 : pari_sp av = avma;
1105 1002608 : long v = get_FpX_var(T), d = 2*(get_FpX_degree(T)-1);
1106 1002608 : GEN z = Z_mod2BIL_ZX(x, bs, d, 0);
1107 1002608 : setvarn(z, v);
1108 1002608 : return gerepileupto(av, FpX_rem(z, T, p));
1109 : }
1110 :
1111 : static GEN
1112 240617 : ZC_mod2BIL_FqC(GEN x, long bs, GEN T, GEN p)
1113 : {
1114 240617 : long i, lx = lg(x);
1115 240617 : GEN A = cgetg(lx, t_COL);
1116 1243225 : for (i=1; i<lx; i++) gel(A,i) = Z_mod2BIL_Fq(gel(x,i), bs, T, p);
1117 240617 : return A;
1118 : }
1119 :
1120 : static GEN
1121 62790 : ZM_mod2BIL_FqM(GEN x, long bs, GEN T, GEN p)
1122 : {
1123 62790 : long i, lx = lg(x);
1124 62790 : GEN A = cgetg(lx, t_MAT);
1125 303407 : for (i=1; i<lx; i++) gel(A,i) = ZC_mod2BIL_FqC(gel(x,i), bs, T, p);
1126 62790 : return A;
1127 : }
1128 :
1129 : GEN
1130 62790 : FqM_mul_Kronecker(GEN x, GEN y, GEN T, GEN p)
1131 : {
1132 62790 : pari_sp av = avma;
1133 62790 : long ex = ZXM_expi(x), ey = ZXM_expi(y), d= degpol(T), n = lg(x)-1;
1134 62790 : long e = ex + ey + expu(d) + expu(n) + 4;
1135 62790 : long N = divsBIL(e)+1;
1136 62790 : GEN z = ZM_mul(ZXM_eval2BIL(x,N), ZXM_eval2BIL(y,N));
1137 62790 : return gerepileupto(av, ZM_mod2BIL_FqM(z, N, T, p));
1138 : }
1139 :
1140 : /*******************************************************************/
1141 : /* */
1142 : /* ZXX */
1143 : /* */
1144 : /*******************************************************************/
1145 :
1146 : void
1147 994 : RgX_check_ZXX(GEN x, const char *s)
1148 : {
1149 994 : long k = lg(x)-1;
1150 9765 : for ( ; k>1; k--) {
1151 8771 : GEN t = gel(x,k);
1152 8771 : switch(typ(t)) {
1153 7749 : case t_INT: break;
1154 1022 : case t_POL: if (RgX_is_ZX(t)) break;
1155 : /* fall through */
1156 0 : default: pari_err_TYPE(stack_strcat(s, " not in Z[X,Y]"),x);
1157 : }
1158 : }
1159 994 : }
1160 :
1161 : /*Renormalize (in place) polynomial with t_INT or ZX coefficients.*/
1162 : GEN
1163 298950709 : ZXX_renormalize(GEN x, long lx)
1164 : {
1165 : long i;
1166 368547218 : for (i = lx-1; i>1; i--)
1167 354712402 : if (signe(gel(x,i))) break;
1168 298950709 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + (i+1)));
1169 299137912 : setlg(x, i+1); setsigne(x, i!=1); return x;
1170 : }
1171 :
1172 : GEN
1173 28548 : ZXX_evalx0(GEN y)
1174 : {
1175 28548 : long i, l = lg(y);
1176 28548 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1177 480693 : for(i=2; i<l; i++)
1178 : {
1179 452145 : GEN yi = gel(y,i);
1180 452145 : gel(z,i) = typ(yi)==t_INT? yi: constant_coeff(yi);
1181 : }
1182 28548 : return ZX_renormalize(z,l);
1183 : }
1184 :
1185 : long
1186 5292 : ZXX_max_lg(GEN x)
1187 : {
1188 5292 : long i, prec = 0, lx = lg(x);
1189 37226 : for (i=2; i<lx; i++)
1190 : {
1191 31934 : GEN p = gel(x,i);
1192 31934 : long l = (typ(p) == t_INT)? lgefint(p): ZX_max_lg(p);
1193 31934 : if (l > prec) prec = l;
1194 : }
1195 5292 : return prec;
1196 : }
1197 :
1198 : GEN
1199 8869 : ZXX_Z_mul(GEN y, GEN x)
1200 : {
1201 8869 : long i, l = lg(y);
1202 8869 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1203 230671 : for(i=2; i<l; i++)
1204 221802 : if(typ(gel(y,i))==t_INT)
1205 0 : gel(z,i) = mulii(gel(y,i),x);
1206 : else
1207 221802 : gel(z,i) = ZX_Z_mul(gel(y,i),x);
1208 8869 : return z;
1209 : }
1210 :
1211 : GEN
1212 0 : ZXX_Z_add_shallow(GEN x, GEN y)
1213 : {
1214 0 : long i, l = lg(x);
1215 : GEN z, a;
1216 0 : if (signe(x)==0) return scalarpol(y,varn(x));
1217 0 : z = cgetg(l,t_POL); z[1] = x[1];
1218 0 : a = gel(x,2);
1219 0 : gel(z, 2) = typ(a)==t_INT? addii(a,y): ZX_Z_add(a,y);
1220 0 : for(i=3; i<l; i++)
1221 0 : gel(z,i) = gel(x,i);
1222 0 : return z;
1223 : }
1224 :
1225 : GEN
1226 55209 : ZXX_Z_divexact(GEN y, GEN x)
1227 : {
1228 55209 : long i, l = lg(y);
1229 55209 : GEN z = cgetg(l,t_POL); z[1] = y[1];
1230 463008 : for(i=2; i<l; i++)
1231 407799 : if(typ(gel(y,i))==t_INT)
1232 2989 : gel(z,i) = diviiexact(gel(y,i),x);
1233 : else
1234 404810 : gel(z,i) = ZX_Z_divexact(gel(y,i),x);
1235 55209 : return z;
1236 : }
1237 :
1238 : /* Kronecker substitution, ZXX -> ZX:
1239 : * P(X,Y) = sum_{0<=i<lP} P_i(X) * Y^i, where deg P_i < n.
1240 : * Returns P(X,X^(2n-1)) */
1241 : GEN
1242 9994026 : RgXX_to_Kronecker_spec(GEN P, long lP, long n)
1243 : {
1244 9994026 : long i, k, N = (n<<1) + 1;
1245 : GEN y;
1246 9994026 : if (!lP) return pol_0(0);
1247 9974725 : y = cgetg((N-2)*lP + 2, t_POL) + 2;
1248 45082576 : for (k=i=0; i<lP; i++)
1249 : {
1250 : long j;
1251 45086634 : GEN c = gel(P,i);
1252 45086634 : if (typ(c)!=t_POL)
1253 : {
1254 3195276 : gel(y,k++) = c;
1255 3195276 : j = 3;
1256 : }
1257 : else
1258 : {
1259 41891358 : long l = lg(c);
1260 41891358 : if (l-3 >= n)
1261 0 : pari_err_BUG("RgXX_to_Kronecker, P is not reduced mod Q");
1262 247390259 : for (j=2; j < l; j++) gel(y,k++) = gel(c,j);
1263 : }
1264 45088180 : if (i == lP-1) break;
1265 195531873 : for ( ; j < N; j++) gel(y,k++) = gen_0;
1266 : }
1267 9974035 : y-=2; setlg(y, k+2); y[1] = evalsigne(1); return y;
1268 : }
1269 :
1270 : /* shallow, n = deg(T) */
1271 : GEN
1272 1764378 : Kronecker_to_ZXX(GEN z, long n, long v)
1273 : {
1274 1764378 : long i,j,lx,l, N = (n<<1)+1;
1275 : GEN x, t;
1276 1764378 : l = lg(z); lx = (l-2) / (N-2);
1277 1764378 : x = cgetg(lx+3,t_POL);
1278 1764376 : x[1] = z[1];
1279 8916708 : for (i=2; i<lx+2; i++)
1280 : {
1281 7151982 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1282 107069381 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1283 7151434 : z += (N-2);
1284 7151434 : gel(x,i) = ZX_renormalize(t,N);
1285 : }
1286 1764726 : N = (l-2) % (N-2) + 2;
1287 1764726 : t = cgetg(N,t_POL); t[1] = evalvarn(v);
1288 5125520 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1289 1764750 : gel(x,i) = ZX_renormalize(t,N);
1290 1764619 : return ZXX_renormalize(x, i+1);
1291 : }
1292 :
1293 : GEN
1294 6649375 : RgXX_to_Kronecker(GEN P, long n)
1295 : {
1296 6649375 : GEN z = RgXX_to_Kronecker_spec(P+2, lgpol(P), n);
1297 6649955 : setvarn(z,varn(P)); return z;
1298 : }
1299 : GEN
1300 3021385 : ZXX_mul_Kronecker(GEN x, GEN y, long n)
1301 3021385 : { return ZX_mul(RgXX_to_Kronecker(x,n), RgXX_to_Kronecker(y,n)); }
1302 :
1303 : GEN
1304 5306 : ZXX_sqr_Kronecker(GEN x, long n)
1305 5306 : { return ZX_sqr(RgXX_to_Kronecker(x,n)); }
1306 :
1307 : /* shallow, n = deg(T) */
1308 : GEN
1309 24143 : Kronecker_to_ZXQX(GEN z, GEN T)
1310 : {
1311 24143 : long i,j,lx,l, N = (degpol(T)<<1)+1;
1312 : GEN x, t;
1313 24143 : l = lg(z); lx = (l-2) / (N-2);
1314 24143 : x = cgetg(lx+3,t_POL);
1315 24143 : x[1] = z[1];
1316 235350 : for (i=2; i<lx+2; i++)
1317 : {
1318 211207 : t = cgetg(N,t_POL); t[1] = T[1];
1319 1889820 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1320 211207 : z += (N-2);
1321 211207 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1322 : }
1323 24143 : N = (l-2) % (N-2) + 2;
1324 24143 : t = cgetg(N,t_POL); t[1] = T[1];
1325 56790 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
1326 24143 : gel(x,i) = ZX_rem(ZX_renormalize(t,N), T);
1327 24143 : return ZXX_renormalize(x, i+1);
1328 : }
1329 :
1330 : GEN
1331 1904 : ZXQX_sqr(GEN x, GEN T)
1332 : {
1333 1904 : pari_sp av = avma;
1334 1904 : long n = degpol(T);
1335 1904 : GEN z = ZXX_sqr_Kronecker(x, n);
1336 1904 : z = Kronecker_to_ZXQX(z, T);
1337 1904 : return gerepilecopy(av, z);
1338 : }
1339 :
1340 : GEN
1341 22239 : ZXQX_mul(GEN x, GEN y, GEN T)
1342 : {
1343 22239 : pari_sp av = avma;
1344 22239 : long n = degpol(T);
1345 22239 : GEN z = ZXX_mul_Kronecker(x, y, n);
1346 22239 : z = Kronecker_to_ZXQX(z, T);
1347 22239 : return gerepilecopy(av, z);
1348 : }
1349 :
1350 : GEN
1351 8862 : ZXQX_ZXQ_mul(GEN P, GEN U, GEN T)
1352 : {
1353 : long i, lP;
1354 : GEN res;
1355 8862 : res = cgetg_copy(P, &lP); res[1] = P[1];
1356 46618 : for(i=2; i<lP; i++)
1357 37757 : gel(res,i) = typ(gel(P,i))==t_POL? ZXQ_mul(U, gel(P,i), T)
1358 37757 : : gmul(U, gel(P,i));
1359 8861 : return ZXX_renormalize(res,lP);
1360 : }
1361 :
1362 : GEN
1363 3446496 : QX_mul(GEN x, GEN y)
1364 : {
1365 3446496 : GEN dx, nx = Q_primitive_part(x, &dx);
1366 3446496 : GEN dy, ny = Q_primitive_part(y, &dy);
1367 3446496 : GEN z = ZX_mul(nx, ny);
1368 3446496 : if (dx || dy)
1369 : {
1370 3357991 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1371 3357991 : return ZX_Q_mul(z, d);
1372 : } else
1373 88505 : return z;
1374 : }
1375 :
1376 : GEN
1377 44099 : QX_sqr(GEN x)
1378 : {
1379 44099 : GEN dx, nx = Q_primitive_part(x, &dx);
1380 44099 : GEN z = ZX_sqr(nx);
1381 44099 : if (dx)
1382 33733 : return ZX_Q_mul(z, gsqr(dx));
1383 : else
1384 10366 : return z;
1385 : }
1386 :
1387 : GEN
1388 1685681 : QX_ZX_rem(GEN x, GEN y)
1389 : {
1390 1685681 : pari_sp av = avma;
1391 1685681 : GEN d, nx = Q_primitive_part(x, &d);
1392 1685681 : GEN r = ZX_rem(nx, y);
1393 1685681 : if (d) r = ZX_Q_mul(r, d);
1394 1685681 : return gerepileupto(av, r);
1395 : }
1396 :
1397 : GEN
1398 13377 : QXQX_mul(GEN x, GEN y, GEN T)
1399 : {
1400 13377 : GEN dx, nx = Q_primitive_part(x, &dx);
1401 13377 : GEN dy, ny = Q_primitive_part(y, &dy);
1402 13377 : GEN z = ZXQX_mul(nx, ny, T);
1403 13377 : if (dx || dy)
1404 : {
1405 8309 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
1406 8309 : return ZXX_Q_mul(z, d);
1407 : } else
1408 5068 : return z;
1409 : }
1410 :
1411 : GEN
1412 1904 : QXQX_sqr(GEN x, GEN T)
1413 : {
1414 1904 : GEN dx, nx = Q_primitive_part(x, &dx);
1415 1904 : GEN z = ZXQX_sqr(nx, T);
1416 1904 : if (dx)
1417 588 : return ZXX_Q_mul(z, gsqr(dx));
1418 : else
1419 1316 : return z;
1420 : }
1421 :
1422 : GEN
1423 476 : QXQX_powers(GEN P, long n, GEN T)
1424 : {
1425 476 : GEN v = cgetg(n+2, t_VEC);
1426 : long i;
1427 476 : gel(v, 1) = pol_1(varn(T));
1428 476 : if (n==0) return v;
1429 476 : gel(v, 2) = gcopy(P);
1430 1554 : for (i = 2; i <= n; i++) gel(v,i+1) = QXQX_mul(P, gel(v,i), T);
1431 476 : return v;
1432 : }
1433 :
1434 : GEN
1435 3465 : QXQX_QXQ_mul(GEN P, GEN U, GEN T)
1436 : {
1437 : long i, lP;
1438 : GEN res;
1439 3465 : res = cgetg_copy(P, &lP); res[1] = P[1];
1440 22715 : for(i=2; i<lP; i++)
1441 19250 : gel(res,i) = typ(gel(P,i))==t_POL? QXQ_mul(U, gel(P,i), T)
1442 19250 : : gmul(U, gel(P,i));
1443 3465 : return ZXX_renormalize(res,lP);
1444 : }
|
