Line data Source code
1 : /* Copyright (C) 2012 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 : /* Not so fast arithmetic with polynomials over FpX */
19 :
20 : /*******************************************************************/
21 : /* */
22 : /* FpXX */
23 : /* */
24 : /*******************************************************************/
25 : /*Polynomials whose coefficients are either polynomials or integers*/
26 :
27 : static GEN
28 53770 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
29 :
30 : static ulong
31 1251301 : to_FlxqX(GEN P, GEN Q, GEN T, GEN p, GEN *pt_P, GEN *pt_Q, GEN *pt_T)
32 : {
33 1251301 : ulong pp = uel(p,2);
34 1251301 : long v = get_FpX_var(T);
35 1251292 : *pt_P = ZXX_to_FlxX(P, pp, v);
36 1251247 : if (pt_Q) *pt_Q = ZXX_to_FlxX(Q, pp, v);
37 1251166 : *pt_T = ZXT_to_FlxT(T, pp);
38 1251265 : return pp;
39 : }
40 :
41 : static GEN
42 126 : ZXX_copy(GEN a) { return gcopy(a); }
43 :
44 : GEN
45 39934 : FpXX_red(GEN z, GEN p)
46 : {
47 : GEN res;
48 39934 : long i, l = lg(z);
49 39934 : res = cgetg(l,t_POL); res[1] = z[1];
50 281451 : for (i=2; i<l; i++)
51 : {
52 241517 : GEN zi = gel(z,i), c;
53 241517 : if (typ(zi)==t_INT)
54 14672 : c = modii(zi,p);
55 : else
56 : {
57 226845 : pari_sp av = avma;
58 226845 : c = FpX_red(zi,p);
59 226845 : switch(lg(c)) {
60 7 : case 2: set_avma(av); c = gen_0; break;
61 20516 : case 3: c = gerepilecopy(av, gel(c,2)); break;
62 : }
63 : }
64 241517 : gel(res,i) = c;
65 : }
66 39934 : return FpXX_renormalize(res,lg(res));
67 : }
68 : GEN
69 446122 : FpXX_add(GEN x, GEN y, GEN p)
70 : {
71 : long i,lz;
72 : GEN z;
73 446122 : long lx=lg(x);
74 446122 : long ly=lg(y);
75 446122 : if (ly>lx) swapspec(x,y, lx,ly);
76 446122 : lz = lx; z = cgetg(lz, t_POL); z[1]=x[1];
77 9338130 : for (i=2; i<ly; i++) gel(z,i) = Fq_add(gel(x,i), gel(y,i), NULL, p);
78 1671109 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
79 446122 : return FpXX_renormalize(z, lz);
80 : }
81 : GEN
82 30887 : FpXX_sub(GEN x, GEN y, GEN p)
83 : {
84 : long i,lz;
85 : GEN z;
86 30887 : long lx=lg(x);
87 30887 : long ly=lg(y);
88 30887 : if (ly <= lx)
89 : {
90 15528 : lz = lx; z = cgetg(lz, t_POL); z[1]=x[1];
91 197919 : for (i=2; i<ly; i++) gel(z,i) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
92 33422 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
93 : }
94 : else
95 : {
96 15359 : lz = ly; z = cgetg(lz, t_POL); z[1]=x[1];
97 94153 : for (i=2; i<lx; i++) gel(z,i) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
98 56661 : for ( ; i<ly; i++) gel(z,i) = Fq_neg(gel(y,i), NULL, p);
99 : }
100 30887 : return FpXX_renormalize(z, lz);
101 : }
102 :
103 : static GEN
104 152221 : FpXX_subspec(GEN x, GEN y, GEN p, long nx, long ny)
105 : {
106 : long i,lz;
107 : GEN z;
108 152221 : if (ny <= nx)
109 : {
110 152221 : lz = nx+2; z = cgetg(lz, t_POL);
111 3714434 : for (i=0; i<ny; i++) gel(z,i+2) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
112 152221 : for ( ; i<nx; i++) gel(z,i+2) = gcopy(gel(x,i));
113 : }
114 : else
115 : {
116 0 : lz = ny+2; z = cgetg(lz, t_POL);
117 0 : for (i=0; i<nx; i++) gel(z,i+2) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
118 0 : for ( ; i<ny; i++) gel(z,i+2) = Fq_neg(gel(y,i), NULL, p);
119 : }
120 152221 : z[1] = 0; return FpXX_renormalize(z, lz);
121 : }
122 :
123 : GEN
124 2073 : FpXX_neg(GEN x, GEN p)
125 : {
126 2073 : long i, lx = lg(x);
127 2073 : GEN y = cgetg(lx,t_POL);
128 2073 : y[1] = x[1];
129 53964 : for(i=2; i<lx; i++) gel(y,i) = Fq_neg(gel(x,i), NULL, p);
130 2073 : return FpXX_renormalize(y, lx);
131 : }
132 :
133 : GEN
134 56452 : FpXX_Fp_mul(GEN P, GEN u, GEN p)
135 : {
136 : long i, lP;
137 56452 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
138 484283 : for(i=2; i<lP; i++)
139 : {
140 427831 : GEN x = gel(P,i);
141 427831 : gel(res,i) = typ(x)==t_INT? Fp_mul(x,u,p): FpX_Fp_mul(x,u,p);
142 : }
143 56452 : return FpXX_renormalize(res,lP);
144 : }
145 :
146 : GEN
147 7116 : FpXX_mulu(GEN P, ulong u, GEN p)
148 : {
149 : long i, lP;
150 7116 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
151 52335 : for(i=2; i<lP; i++)
152 : {
153 45219 : GEN x = gel(P,i);
154 45219 : gel(res,i) = typ(x)==t_INT? Fp_mulu(x,u,p): FpX_mulu(x,u,p);
155 : }
156 7116 : return FpXX_renormalize(res,lP);
157 : }
158 :
159 : GEN
160 2240 : FpXX_halve(GEN P, GEN p)
161 : {
162 : long i, lP;
163 2240 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
164 7770 : for(i=2; i<lP; i++)
165 : {
166 5530 : GEN x = gel(P,i);
167 5530 : gel(res,i) = typ(x)==t_INT? Fp_halve(x,p): FpX_halve(x,p);
168 : }
169 2240 : return FpXX_renormalize(res,lP);
170 : }
171 :
172 : GEN
173 13405 : FpXX_deriv(GEN P, GEN p)
174 : {
175 13405 : long i, l = lg(P)-1;
176 : GEN res;
177 :
178 13405 : if (l < 3) return pol_0(varn(P));
179 13055 : res = cgetg(l, t_POL);
180 13055 : res[1] = P[1];
181 80960 : for (i=2; i<l ; i++)
182 : {
183 67905 : GEN x = gel(P,i+1);
184 67905 : gel(res,i) = typ(x)==t_INT? Fp_mulu(x,i-1,p): FpX_mulu(x,i-1,p);
185 : }
186 13055 : return FpXX_renormalize(res, l);
187 : }
188 :
189 : GEN
190 0 : FpXX_integ(GEN P, GEN p)
191 : {
192 0 : long i, l = lg(P);
193 : GEN res;
194 :
195 0 : if (l == 2) return pol_0(varn(P));
196 0 : res = cgetg(l+1, t_POL);
197 0 : res[1] = P[1];
198 0 : gel(res,2) = gen_0;
199 0 : for (i=3; i<=l ; i++)
200 : {
201 0 : GEN x = gel(P,i-1);
202 0 : if (signe(x))
203 : {
204 0 : GEN i1 = Fp_inv(utoi(i-2), p);
205 0 : gel(res,i) = typ(x)==t_INT? Fp_mul(x,i1,p): FpX_Fp_mul(x,i1,p);
206 : } else
207 0 : gel(res,i) = gen_0;
208 : }
209 0 : return FpXX_renormalize(res, l+1);
210 : }
211 :
212 : /*******************************************************************/
213 : /* */
214 : /* (Fp[X]/(Q))[Y] */
215 : /* */
216 : /*******************************************************************/
217 :
218 : static GEN
219 1332728 : get_FpXQX_red(GEN T, GEN *B)
220 : {
221 1332728 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
222 94922 : *B = gel(T,1); return gel(T,2);
223 : }
224 :
225 : GEN
226 52 : random_FpXQX(long d1, long v, GEN T, GEN p)
227 : {
228 52 : long dT = get_FpX_degree(T), vT = get_FpX_var(T);
229 52 : long i, d = d1+2;
230 52 : GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
231 284 : for (i=2; i<d; i++) gel(y,i) = random_FpX(dT, vT, p);
232 52 : return FpXQX_renormalize(y,d);
233 : }
234 :
235 : /*Not stack clean*/
236 : GEN
237 1891094 : Kronecker_to_FpXQX(GEN Z, GEN T, GEN p)
238 : {
239 1891094 : long i,j,lx,l, N = (get_FpX_degree(T)<<1) + 1;
240 1891072 : GEN x, t = cgetg(N,t_POL), z = FpX_red(Z, p);
241 1890177 : t[1] = evalvarn(get_FpX_var(T));
242 1890212 : l = lg(z); lx = (l-2) / (N-2);
243 1890212 : x = cgetg(lx+3,t_POL);
244 1890748 : x[1] = z[1];
245 29935944 : for (i=2; i<lx+2; i++)
246 : {
247 238540326 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
248 28045137 : z += (N-2);
249 28045137 : gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
250 : }
251 1890807 : N = (l-2) % (N-2) + 2;
252 3393955 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
253 1890807 : gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
254 1890870 : return FpXQX_renormalize(x, i+1);
255 : }
256 :
257 : GEN
258 1937270 : FpXQX_red(GEN z, GEN T, GEN p)
259 : {
260 1937270 : long i, l = lg(z);
261 1937270 : GEN res = cgetg(l,t_POL); res[1] = z[1];
262 16302134 : for(i=2;i<l;i++)
263 14368425 : if (typ(gel(z,i)) == t_INT)
264 160149 : gel(res,i) = modii(gel(z,i),p);
265 : else
266 14208276 : gel(res,i) = FpXQ_red(gel(z,i),T,p);
267 1933709 : return FpXQX_renormalize(res,l);
268 : }
269 :
270 : GEN
271 0 : FpXQXV_red(GEN x, GEN T, GEN p)
272 0 : { pari_APPLY_type(t_VEC, FpXQX_red(gel(x,i), T, p)) }
273 :
274 : GEN
275 0 : FpXQXT_red(GEN x, GEN T, GEN p)
276 : {
277 0 : if (typ(x) == t_POL)
278 0 : return FpXQX_red(x, T, p);
279 : else
280 0 : pari_APPLY_type(t_VEC, FpXQXT_red(gel(x,i), T, p))
281 : }
282 :
283 : static GEN
284 2191 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
285 :
286 : GEN
287 532 : FpXQX_to_mod(GEN z, GEN T, GEN p)
288 : {
289 532 : long i, l = lg(z);
290 : GEN x;
291 532 : if (l == 2)
292 : {
293 7 : x = cgetg(3, t_POL); x[1] = z[1];
294 7 : p = icopy(p); T = FpX_to_mod_raw(T, p);
295 7 : gel(x,2) = mkpolmod(mkintmod(gen_0, p), T);
296 7 : return x;
297 : }
298 525 : x = cgetg(l, t_POL); x[1] = z[1];
299 525 : p = icopy(p); T = FpX_to_mod_raw(T, p);
300 6692 : for (i=2; i<l; i++)
301 : {
302 6167 : GEN zi = gel(z,i);
303 6167 : gel(x,i) = typ(zi) == t_POL? mkpolmod(FpX_to_mod_raw(zi, p), T)
304 6167 : : to_intmod(zi, p);
305 : }
306 525 : return normalizepol_lg(x,l);
307 : }
308 :
309 : static GEN
310 0 : FpXQX_to_mod_raw(GEN z, GEN T, GEN p)
311 : {
312 0 : long i, l = lg(z);
313 : GEN x;
314 :
315 0 : if (l == 2)
316 : {
317 0 : x = cgetg(3, t_POL); x[1] = z[1];
318 0 : p = icopy(p);
319 0 : gel(x,2) = mkpolmod(mkintmod(gen_0, p), T);
320 0 : return x;
321 : }
322 0 : x = cgetg(l, t_POL); x[1] = z[1];
323 0 : for (i=2; i<l; i++)
324 : {
325 0 : GEN zi = gel(z,i);
326 0 : gel(x,i) = typ(zi) == t_POL? mkpolmod(FpX_to_mod_raw(zi, p), T)
327 0 : : to_intmod(zi, p);
328 : }
329 0 : return normalizepol_lg(x,l);
330 : }
331 :
332 : INLINE GEN
333 0 : FqX_to_mod_raw(GEN f, GEN T, GEN p)
334 0 : { return T?FpXQX_to_mod_raw(f, T, p): FpX_to_mod_raw(f, p); }
335 :
336 : static GEN
337 0 : FqXC_to_mod_raw(GEN x, GEN T, GEN p)
338 0 : { pari_APPLY_type(t_COL, FqX_to_mod_raw(gel(x,i), T, p)) }
339 :
340 : GEN
341 14 : FqXC_to_mod(GEN z, GEN T, GEN p)
342 : {
343 : GEN x;
344 14 : long i,l = lg(z);
345 14 : if (!T) return FpXC_to_mod(z, p);
346 0 : x = cgetg(l, t_COL);
347 0 : if (l == 1) return x;
348 0 : p = icopy(p);
349 0 : T = FpX_to_mod_raw(T, p);
350 0 : for (i=1; i<l; i++)
351 0 : gel(x,i) = FqX_to_mod_raw(gel(z, i), T, p);
352 0 : return x;
353 : }
354 :
355 : GEN
356 0 : FqXM_to_mod(GEN z, GEN T, GEN p)
357 : {
358 : GEN x;
359 0 : long i,l = lg(z);
360 0 : if (!T) return FpXM_to_mod(z, p);
361 0 : x = cgetg(l, t_MAT);
362 0 : if (l == 1) return x;
363 0 : p = icopy(p);
364 0 : T = FpX_to_mod_raw(T, p);
365 0 : for (i=1; i<l; i++)
366 0 : gel(x,i) = FqXC_to_mod_raw(gel(z, i), T, p);
367 0 : return x;
368 : }
369 :
370 : static int
371 3639930 : ZXX_is_ZX_spec(GEN a,long na)
372 : {
373 : long i;
374 3951408 : for(i=0;i<na;i++)
375 3890855 : if(typ(gel(a,i))!=t_INT) return 0;
376 60553 : return 1;
377 : }
378 :
379 : static int
380 246030 : ZXX_is_ZX(GEN a) { return ZXX_is_ZX_spec(a+2,lgpol(a)); }
381 :
382 : static GEN
383 140824 : FpXX_FpX_mulspec(GEN P, GEN U, GEN p, long v, long lU)
384 : {
385 140824 : long i, lP =lg(P);
386 : GEN res;
387 140824 : res = cgetg(lP, t_POL); res[1] = P[1];
388 7696178 : for(i=2; i<lP; i++)
389 : {
390 7555354 : GEN Pi = gel(P,i);
391 7555354 : gel(res,i) = typ(Pi)==t_INT? FpX_Fp_mulspec(U, Pi, p, lU):
392 7540677 : FpX_mulspec(U, Pi+2, p, lU, lgpol(Pi));
393 7555354 : setvarn(gel(res,i),v);
394 : }
395 140824 : return FpXQX_renormalize(res,lP);
396 : }
397 :
398 : GEN
399 125645 : FpXX_FpX_mul(GEN P, GEN U, GEN p)
400 125645 : { return FpXX_FpX_mulspec(P,U+2,p,varn(U),lgpol(U)); }
401 :
402 : static GEN
403 15179 : FpXY_FpY_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
404 : {
405 15179 : pari_sp av = avma;
406 15179 : long v = fetch_var();
407 15179 : GEN z = RgXY_swapspec(x,get_FpX_degree(T)-1,v,lx);
408 15179 : z = FpXX_FpX_mulspec(z,y,p,v,ly);
409 15179 : z = RgXY_swapspec(z+2,lx+ly+3,get_FpX_var(T),lgpol(z));
410 15179 : (void)delete_var(); return gerepilecopy(av,z);
411 : }
412 :
413 : static GEN
414 1696824 : FpXQX_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
415 : {
416 1696824 : pari_sp av = avma;
417 : GEN z, kx, ky;
418 : long n;
419 1696824 : if (ZXX_is_ZX_spec(y,ly))
420 : {
421 15692 : if (ZXX_is_ZX_spec(x,lx))
422 8271 : return FpX_mulspec(x,y,p,lx,ly);
423 : else
424 7421 : return FpXY_FpY_mulspec(x,y,T,p,lx,ly);
425 1681331 : } else if (ZXX_is_ZX_spec(x,lx))
426 7758 : return FpXY_FpY_mulspec(y,x,T,p,ly,lx);
427 1673823 : n = get_FpX_degree(T);
428 1673847 : kx = RgXX_to_Kronecker_spec(x, lx, n);
429 1673913 : ky = RgXX_to_Kronecker_spec(y, ly, n);
430 1674045 : z = Kronecker_to_FpXQX(ZX_mul(ky,kx), T, p);
431 1673684 : return gerepileupto(av, z);
432 : }
433 :
434 : GEN
435 1384532 : FpXQX_mul(GEN x, GEN y, GEN T, GEN p)
436 : {
437 1384532 : GEN z = FpXQX_mulspec(x+2,y+2,T,p,lgpol(x),lgpol(y));
438 1385077 : setvarn(z,varn(x)); return z;
439 : }
440 :
441 : GEN
442 186666 : FpXQX_sqr(GEN x, GEN T, GEN p)
443 : {
444 186666 : pari_sp av = avma;
445 : GEN z, kx;
446 186666 : if (ZXX_is_ZX(x)) return ZX_sqr(x);
447 179492 : kx= RgXX_to_Kronecker(x, get_FpX_degree(T));
448 179492 : z = Kronecker_to_FpXQX(ZX_sqr(kx), T, p);
449 179492 : return gerepileupto(av, z);
450 : }
451 :
452 : GEN
453 553161 : FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p)
454 : {
455 : long i, lP;
456 : GEN res;
457 553161 : res = cgetg_copy(P, &lP); res[1] = P[1];
458 2065501 : for(i=2; i<lP; i++)
459 2747380 : gel(res,i) = typ(gel(P,i))==t_INT? FpX_Fp_mul(U, gel(P,i), p):
460 1234997 : FpXQ_mul(U, gel(P,i), T,p);
461 552975 : return FpXQX_renormalize(res,lP);
462 : }
463 :
464 : /* x and y in Z[Y][X]. Assume T irreducible mod p */
465 : static GEN
466 173365 : FpXQX_divrem_basecase(GEN x, GEN y, GEN T, GEN p, GEN *pr)
467 : {
468 : long vx, dx, dy, dy1, dz, i, j, sx, lr;
469 : pari_sp av0, av, tetpil;
470 : GEN z,p1,rem,lead;
471 :
472 173365 : if (!signe(y)) pari_err_INV("FpX_divrem",y);
473 173365 : vx=varn(x); dy=degpol(y); dx=degpol(x);
474 173365 : if (dx < dy)
475 : {
476 185 : if (pr)
477 : {
478 135 : av0 = avma; x = FpXQX_red(x, T, p);
479 135 : if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: pol_0(vx); }
480 135 : if (pr == ONLY_REM) return x;
481 135 : *pr = x;
482 : }
483 185 : return pol_0(vx);
484 : }
485 173180 : lead = leading_coeff(y);
486 173180 : if (!dy) /* y is constant */
487 : {
488 1408 : if (pr && pr != ONLY_DIVIDES)
489 : {
490 1044 : if (pr == ONLY_REM) return pol_0(vx);
491 7 : *pr = pol_0(vx);
492 : }
493 371 : if (gequal1(lead)) return FpXQX_red(x,T,p);
494 355 : av0 = avma; x = FqX_Fq_mul(x, Fq_inv(lead, T,p), T,p);
495 355 : return gerepileupto(av0,x);
496 : }
497 171772 : av0 = avma; dz = dx-dy;
498 171772 : lead = gequal1(lead)? NULL: gclone(Fq_inv(lead,T,p));
499 171772 : set_avma(av0);
500 171772 : z = cgetg(dz+3,t_POL); z[1] = x[1];
501 171772 : x += 2; y += 2; z += 2;
502 178250 : for (dy1=dy-1; dy1>=0 && !signe(gel(y, dy1)); dy1--);
503 :
504 171772 : p1 = gel(x,dx); av = avma;
505 171772 : gel(z,dz) = lead? gerepileupto(av, Fq_mul(p1,lead, T, p)): gcopy(p1);
506 521131 : for (i=dx-1; i>=dy; i--)
507 : {
508 349359 : av=avma; p1=gel(x,i);
509 1164344 : for (j=i-dy1; j<=i && j<=dz; j++)
510 814985 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
511 349359 : if (lead) p1 = Fq_mul(p1, lead, NULL,p);
512 349359 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil,Fq_red(p1,T,p));
513 : }
514 171772 : if (!pr) { guncloneNULL(lead); return z-2; }
515 :
516 168613 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
517 180656 : for (sx=0; ; i--)
518 : {
519 180656 : p1 = gel(x,i);
520 692303 : for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
521 511647 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
522 180656 : tetpil=avma; p1 = Fq_red(p1, T, p); if (signe(p1)) { sx = 1; break; }
523 13968 : if (!i) break;
524 12043 : set_avma(av);
525 : }
526 168613 : if (pr == ONLY_DIVIDES)
527 : {
528 0 : guncloneNULL(lead);
529 0 : if (sx) return gc_NULL(av0);
530 0 : return gc_const((pari_sp)rem, z-2);
531 : }
532 168613 : lr=i+3; rem -= lr;
533 168613 : rem[0] = evaltyp(t_POL) | _evallg(lr);
534 168613 : rem[1] = z[-1];
535 168613 : p1 = gerepile((pari_sp)rem,tetpil,p1);
536 168613 : rem += 2; gel(rem,i) = p1;
537 1479834 : for (i--; i>=0; i--)
538 : {
539 1311221 : av=avma; p1 = gel(x,i);
540 4109436 : for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
541 2798215 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j), NULL,p), NULL,p);
542 1311221 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, Fq_red(p1, T, p));
543 : }
544 168613 : rem -= 2;
545 168613 : guncloneNULL(lead);
546 168613 : if (!sx) (void)FpXQX_renormalize(rem, lr);
547 168613 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
548 13546 : *pr = rem; return z-2;
549 : }
550 :
551 : static GEN
552 752 : FpXQX_addmulmul(GEN u, GEN v, GEN x, GEN y, GEN T, GEN p)
553 : {
554 752 : return FpXX_add(FpXQX_mul(u, x, T, p),FpXQX_mul(v, y, T, p), p);
555 : }
556 :
557 : static GEN
558 376 : FpXQXM_FpXQX_mul2(GEN M, GEN x, GEN y, GEN T, GEN p)
559 : {
560 376 : GEN res = cgetg(3, t_COL);
561 376 : gel(res, 1) = FpXQX_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, T, p);
562 376 : gel(res, 2) = FpXQX_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, T, p);
563 376 : return res;
564 : }
565 :
566 : static GEN
567 161 : FpXQXM_mul2(GEN A, GEN B, GEN T, GEN p)
568 : {
569 161 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
570 161 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
571 161 : GEN M1 = FpXQX_mul(FpXX_add(A11,A22, p), FpXX_add(B11,B22, p), T, p);
572 161 : GEN M2 = FpXQX_mul(FpXX_add(A21,A22, p), B11, T, p);
573 161 : GEN M3 = FpXQX_mul(A11, FpXX_sub(B12,B22, p), T, p);
574 161 : GEN M4 = FpXQX_mul(A22, FpXX_sub(B21,B11, p), T, p);
575 161 : GEN M5 = FpXQX_mul(FpXX_add(A11,A12, p), B22, T, p);
576 161 : GEN M6 = FpXQX_mul(FpXX_sub(A21,A11, p), FpXX_add(B11,B12, p), T, p);
577 161 : GEN M7 = FpXQX_mul(FpXX_sub(A12,A22, p), FpXX_add(B21,B22, p), T, p);
578 161 : GEN T1 = FpXX_add(M1,M4, p), T2 = FpXX_sub(M7,M5, p);
579 161 : GEN T3 = FpXX_sub(M1,M2, p), T4 = FpXX_add(M3,M6, p);
580 161 : retmkmat22(FpXX_add(T1,T2, p), FpXX_add(M3,M5, p),
581 : FpXX_add(M2,M4, p), FpXX_add(T3,T4, p));
582 : }
583 : /* Return [0,1;1,-q]*M */
584 : static GEN
585 161 : FpXQX_FpXQXM_qmul(GEN q, GEN M, GEN T, GEN p)
586 : {
587 161 : GEN u = FpXQX_mul(gcoeff(M,2,1), q, T, p);
588 161 : GEN v = FpXQX_mul(gcoeff(M,2,2), q, T, p);
589 161 : retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
590 : FpXX_sub(gcoeff(M,1,1), u, p), FpXX_sub(gcoeff(M,1,2), v, p));
591 : }
592 :
593 : static GEN
594 0 : matid2_FpXQXM(long v)
595 0 : { retmkmat22(pol_1(v),pol_0(v),pol_0(v),pol_1(v)); }
596 :
597 : static GEN
598 0 : matJ2_FpXQXM(long v)
599 0 : { retmkmat22(pol_0(v),pol_1(v),pol_1(v),pol_0(v)); }
600 :
601 : static GEN
602 19504 : FpXX_shift(GEN a, long n) { return RgX_shift_shallow(a, n); }
603 :
604 : INLINE GEN
605 8442 : FpXXn_red(GEN a, long n) { return RgXn_red_shallow(a, n); }
606 :
607 : /* Fast resultant formula from William Hart in Flint <http://flintlib.org/> */
608 :
609 : struct FpXQX_res
610 : {
611 : GEN res, lc;
612 : long deg0, deg1, off;
613 : };
614 :
615 : INLINE void
616 0 : FpXQX_halfres_update(long da, long db, long dr, GEN T, GEN p, struct FpXQX_res *res)
617 : {
618 0 : if (dr >= 0)
619 : {
620 0 : if (!ZX_equal1(res->lc))
621 : {
622 0 : res->lc = FpXQ_powu(res->lc, da - dr, T, p);
623 0 : res->res = FpXQ_mul(res->res, res->lc, T, p);
624 : }
625 0 : if (both_odd(da + res->off, db + res->off))
626 0 : res->res = FpX_neg(res->res, p);
627 : } else
628 : {
629 0 : if (db == 0)
630 : {
631 0 : if (!ZX_equal1(res->lc))
632 : {
633 0 : res->lc = FpXQ_powu(res->lc, da, T, p);
634 0 : res->res = FpXQ_mul(res->res, res->lc, T, p);
635 : }
636 : } else
637 0 : res->res = pol_0(get_FpX_var(T));
638 : }
639 0 : }
640 :
641 : static GEN
642 275 : FpXQX_halfres_basecase(GEN a, GEN b, GEN T, GEN p, GEN *pa, GEN *pb, struct FpXQX_res *res)
643 : {
644 275 : pari_sp av=avma;
645 : GEN u,u1,v,v1, M;
646 275 : long vx = varn(a), vT = get_FpX_var(T), n = lgpol(a)>>1;
647 275 : u1 = v = pol_0(vx);
648 275 : u = v1 = pol_1(vx);
649 2837 : while (lgpol(b)>n)
650 : {
651 : GEN r, q;
652 2562 : q = FpXQX_divrem(a,b, T, p, &r);
653 2562 : if (res)
654 : {
655 0 : long da = degpol(a), db=degpol(b), dr = degpol(r);
656 0 : res->lc = to_ZX(gel(b,db+2),vT);
657 0 : if (dr >= n)
658 0 : FpXQX_halfres_update(da, db, dr, T, p, res);
659 : else
660 : {
661 0 : res->deg0 = da;
662 0 : res->deg1 = db;
663 : }
664 : }
665 2562 : a = b; b = r; swap(u,u1); swap(v,v1);
666 2562 : u1 = FpXX_sub(u1, FpXQX_mul(u, q, T, p), p);
667 2562 : v1 = FpXX_sub(v1, FpXQX_mul(v, q, T, p), p);
668 2562 : if (gc_needed(av,2))
669 : {
670 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_halfgcd (d = %ld)",degpol(b));
671 0 : gerepileall(av,res ? 8: 6, &a,&b,&u1,&v1,&u,&v,&res->res,&res->lc);
672 : }
673 : }
674 275 : M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
675 275 : return gc_all(av, res ? 5: 3, &M, pa, pb, &res->res, &res->lc);
676 : }
677 :
678 : static GEN FpXQX_halfres_i(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res);
679 :
680 : static GEN
681 215 : FpXQX_halfres_split(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res)
682 : {
683 215 : pari_sp av = avma;
684 : GEN Q, R, S, V1, V2;
685 : GEN x1, y1, r, q;
686 215 : long l = lgpol(x), n = l>>1, k, vT = get_FpX_var(T);
687 215 : if (lgpol(y) <= n)
688 0 : { *a = RgX_copy(x); *b = RgX_copy(y); return matid2_FpXQXM(varn(x)); }
689 215 : if (res)
690 : {
691 0 : res->lc = to_ZX(leading_coeff(y), vT);
692 0 : res->deg0 -= n;
693 0 : res->deg1 -= n;
694 0 : res->off += n;
695 : }
696 215 : R = FpXQX_halfres_i(FpXX_shift(x,-n),FpXX_shift(y,-n), T, p, a, b, res);
697 215 : if (res)
698 : {
699 0 : res->off -= n;
700 0 : res->deg0 += n;
701 0 : res->deg1 += n;
702 : }
703 215 : V1 = FpXQXM_FpXQX_mul2(R, Flxn_red(x,n), Flxn_red(y,n), T, p);
704 215 : x1 = FpXX_add(FpXX_shift(*a,n), gel(V1,1), p);
705 215 : y1 = FpXX_add(FpXX_shift(*b,n), gel(V1,2), p);
706 215 : if (lgpol(y1) <= n)
707 : {
708 54 : *a = x1; *b = y1;
709 54 : return gc_all(av, res ? 5: 3, &R, a, b, &res->res, &res->lc);
710 : }
711 161 : k = 2*n-degpol(y1);
712 161 : q = FpXQX_divrem(x1, y1, T, p, &r);
713 161 : if (res)
714 : {
715 0 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
716 0 : if (dy1 < degpol(y))
717 0 : FpXQX_halfres_update(res->deg0, res->deg1, dy1, T, p, res);
718 0 : res->lc = to_ZX(leading_coeff(y1), vT);
719 0 : res->deg0 = dx1;
720 0 : res->deg1 = dy1;
721 0 : if (dr >= n)
722 : {
723 0 : FpXQX_halfres_update(dx1, dy1, dr, T, p, res);
724 0 : res->deg0 = dy1;
725 0 : res->deg1 = dr;
726 : }
727 0 : res->deg0 -= k;
728 0 : res->deg1 -= k;
729 0 : res->off += k;
730 : }
731 161 : S = FpXQX_halfres_i(FpXX_shift(y1,-k), FpXX_shift(r,-k), T, p, a, b, res);
732 161 : if (res)
733 : {
734 0 : res->deg0 += k;
735 0 : res->deg1 += k;
736 0 : res->off -= k;
737 : }
738 161 : Q = FpXQXM_mul2(S,FpXQX_FpXQXM_qmul(q, R, T, p), T, p);
739 161 : V2 = FpXQXM_FpXQX_mul2(S, FpXXn_red(y1,k), FpXXn_red(r,k), T, p);
740 161 : *a = FpXX_add(FpXX_shift(*a,k), gel(V2,1), p);
741 161 : *b = FpXX_add(FpXX_shift(*b,k), gel(V2,2), p);
742 161 : return gc_all(av, res ? 5: 3, &Q, a, b, &res->res, &res->lc);
743 : }
744 :
745 : static GEN
746 490 : FpXQX_halfres_i(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res)
747 : {
748 490 : if (lgpol(x) < FpXQX_HALFGCD_LIMIT)
749 275 : return FpXQX_halfres_basecase(x, y, T, p, a, b, res);
750 215 : return FpXQX_halfres_split(x, y, T, p, a, b, res);
751 : }
752 :
753 : static GEN
754 114 : FpXQX_halfgcd_all_i(GEN x, GEN y, GEN T, GEN p, GEN *pa, GEN *pb)
755 : {
756 : GEN a, b;
757 114 : GEN R = FpXQX_halfres_i(x, y, T, p, &a, &b, NULL);
758 114 : if (pa) *pa = a;
759 114 : if (pb) *pb = b;
760 114 : return R;
761 : }
762 :
763 : /* Return M in GL_2(Fp[X]/(T)[Y]) such that:
764 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
765 : */
766 :
767 : GEN
768 114 : FpXQX_halfgcd_all(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b)
769 : {
770 114 : pari_sp av = avma;
771 : GEN R,q,r;
772 114 : if (lgefint(p)==3)
773 : {
774 0 : ulong pp = to_FlxqX(x, y, T, p, &x, &y, &T);
775 0 : R = FlxXM_to_ZXXM(FlxqX_halfgcd(x, y, T, pp));
776 0 : if (a) *a = Flx_to_ZX(*a);
777 0 : if (b) *b = Flx_to_ZX(*b);
778 0 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
779 : }
780 114 : if (!signe(x))
781 : {
782 0 : if (a) *a = RgX_copy(y);
783 0 : if (b) *b = RgX_copy(x);
784 0 : return matJ2_FpXQXM(varn(x));
785 : }
786 114 : if (degpol(y)<degpol(x)) return FpXQX_halfgcd_all_i(x, y, T, p, a, b);
787 26 : q = FpXQX_divrem(y, x, T, p, &r);
788 26 : R = FpXQX_halfgcd_all_i(x, r, T, p, a, b);
789 26 : gcoeff(R,1,1) = FpXX_sub(gcoeff(R,1,1),
790 26 : FpXQX_mul(q, gcoeff(R,1,2), T, p), p);
791 26 : gcoeff(R,2,1) = FpXX_sub(gcoeff(R,2,1),
792 26 : FpXQX_mul(q, gcoeff(R,2,2), T, p), p);
793 26 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
794 : }
795 :
796 : GEN
797 44 : FpXQX_halfgcd(GEN x, GEN y, GEN T, GEN p)
798 44 : { return FpXQX_halfgcd_all(x, y, T, p, NULL, NULL); }
799 :
800 : static GEN
801 3896 : FpXQX_gcd_basecase(GEN a, GEN b, GEN T, GEN p)
802 : {
803 3896 : pari_sp av = avma, av0=avma;
804 39080 : while (signe(b))
805 : {
806 : GEN c;
807 35184 : if (gc_needed(av0,2))
808 : {
809 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_gcd (d = %ld)",degpol(b));
810 0 : gerepileall(av0,2, &a,&b);
811 : }
812 35184 : av = avma; c = FpXQX_rem(a, b, T, p); a=b; b=c;
813 : }
814 3896 : return gc_const(av, a);
815 : }
816 :
817 : GEN
818 14545 : FpXQX_gcd(GEN x, GEN y, GEN T, GEN p)
819 : {
820 14545 : pari_sp av = avma;
821 14545 : if (lgefint(p) == 3)
822 : {
823 : GEN Pl, Ql, Tl, U;
824 10540 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
825 10540 : U = FlxqX_gcd(Pl, Ql, Tl, pp);
826 10540 : return gerepileupto(av, FlxX_to_ZXX(U));
827 : }
828 4005 : x = FpXQX_red(x, T, p);
829 4005 : y = FpXQX_red(y, T, p);
830 4005 : if (!signe(x)) return gerepileupto(av, y);
831 3959 : while (lgpol(y)>=FpXQX_GCD_LIMIT)
832 : {
833 63 : if (lgpol(y)<=(lgpol(x)>>1))
834 : {
835 0 : GEN r = FpXQX_rem(x, y, T, p);
836 0 : x = y; y = r;
837 : }
838 63 : (void) FpXQX_halfgcd_all(x,y, T, p, &x, &y);
839 63 : if (gc_needed(av,2))
840 : {
841 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_gcd (y = %ld)",degpol(y));
842 0 : gerepileall(av,2,&x,&y);
843 : }
844 : }
845 3896 : return gerepileupto(av, FpXQX_gcd_basecase(x, y, T, p));
846 : }
847 :
848 : static GEN
849 0 : FpXQX_extgcd_basecase(GEN a, GEN b, GEN T, GEN p, GEN *ptu, GEN *ptv)
850 : {
851 0 : pari_sp av=avma;
852 : GEN u,v,d,d1,v1;
853 0 : long vx = varn(a);
854 0 : d = a; d1 = b;
855 0 : v = pol_0(vx); v1 = pol_1(vx);
856 0 : while (signe(d1))
857 : {
858 0 : GEN r, q = FpXQX_divrem(d, d1, T, p, &r);
859 0 : v = FpXX_sub(v,FpXQX_mul(q,v1,T, p),p);
860 0 : u=v; v=v1; v1=u;
861 0 : u=r; d=d1; d1=u;
862 0 : if (gc_needed(av,2))
863 : {
864 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_extgcd (d = %ld)",degpol(d));
865 0 : gerepileall(av,5, &d,&d1,&u,&v,&v1);
866 : }
867 : }
868 0 : if (ptu) *ptu = FpXQX_div(FpXX_sub(d,FpXQX_mul(b,v, T, p), p), a, T, p);
869 0 : *ptv = v; return d;
870 : }
871 :
872 : static GEN
873 0 : FpXQX_extgcd_halfgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
874 : {
875 : GEN u,v;
876 0 : GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
877 0 : long i, n = 0, vs = varn(x);
878 0 : while (lgpol(y) >= FpXQX_EXTGCD_LIMIT)
879 : {
880 0 : if (lgpol(y)<=(lgpol(x)>>1))
881 : {
882 0 : GEN r, q = FpXQX_divrem(x, y, T, p, &r);
883 0 : x = y; y = r;
884 0 : gel(V,++n) = mkmat22(pol_0(vs),pol_1(vs),pol_1(vs),FpXX_neg(q,p));
885 : } else
886 0 : gel(V,++n) = FpXQX_halfgcd_all(x, y, T, p, &x, &y);
887 : }
888 0 : y = FpXQX_extgcd_basecase(x,y, T, p, &u,&v);
889 0 : for (i = n; i>1; i--)
890 : {
891 0 : GEN R = gel(V,i);
892 0 : GEN u1 = FpXQX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), T, p);
893 0 : GEN v1 = FpXQX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), T, p);
894 0 : u = u1; v = v1;
895 : }
896 : {
897 0 : GEN R = gel(V,1);
898 0 : if (ptu)
899 0 : *ptu = FpXQX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), T, p);
900 0 : *ptv = FpXQX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), T, p);
901 : }
902 0 : return y;
903 : }
904 :
905 : /* x and y in Z[Y][X], return lift(gcd(x mod T,p, y mod T,p)). Set u and v st
906 : * ux + vy = gcd (mod T,p) */
907 : GEN
908 133609 : FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
909 : {
910 133609 : pari_sp av = avma;
911 : GEN d;
912 133609 : if (lgefint(p) == 3)
913 : {
914 : GEN Pl, Ql, Tl, Dl;
915 133609 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
916 133613 : Dl = FlxqX_extgcd(Pl, Ql, Tl, pp, ptu, ptv);
917 133616 : if (ptu) *ptu = FlxX_to_ZXX(*ptu);
918 133614 : *ptv = FlxX_to_ZXX(*ptv);
919 133614 : d = FlxX_to_ZXX(Dl);
920 : }
921 : else
922 : {
923 0 : x = FpXQX_red(x, T, p);
924 0 : y = FpXQX_red(y, T, p);
925 0 : if (lgpol(y)>=FpXQX_EXTGCD_LIMIT)
926 0 : d = FpXQX_extgcd_halfgcd(x, y, T, p, ptu, ptv);
927 : else
928 0 : d = FpXQX_extgcd_basecase(x, y, T, p, ptu, ptv);
929 : }
930 133608 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
931 : }
932 :
933 : static GEN
934 0 : FpXQX_halfres(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, GEN *r)
935 : {
936 : struct FpXQX_res res;
937 : GEN V;
938 0 : long dB, vT = get_FpX_var(T);
939 :
940 0 : res.res = *r;
941 0 : res.lc = to_ZX(leading_coeff(y),vT);
942 0 : res.deg0 = degpol(x);
943 0 : res.deg1 = degpol(y);
944 0 : res.off = 0;
945 0 : V = FpXQX_halfres_i(x, y, T, p, a, b, &res);
946 0 : dB = degpol(*b);
947 0 : if (dB < degpol(y))
948 0 : FpXQX_halfres_update(res.deg0, res.deg1, dB, T, p, &res);
949 0 : *r = res.res;
950 0 : return V;
951 : }
952 :
953 : /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab), with R=A%B, a=deg(A) ...*/
954 : static GEN
955 28 : FpXQX_resultant_basecase(GEN a, GEN b, GEN T, GEN p)
956 : {
957 28 : pari_sp av = avma;
958 28 : long vT = get_FpX_var(T), da,db,dc;
959 28 : GEN c,lb, res = pol_1(vT);
960 :
961 28 : if (!signe(a) || !signe(b)) return pol_0(vT);
962 :
963 28 : da = degpol(a);
964 28 : db = degpol(b);
965 28 : if (db > da)
966 : {
967 0 : swapspec(a,b, da,db);
968 0 : if (both_odd(da,db)) res = FpX_neg(res, p);
969 : }
970 28 : if (!da) return pol_1(vT); /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
971 98 : while (db)
972 : {
973 70 : lb = to_ZX(gel(b,db+2),vT);
974 70 : c = FpXQX_rem(a,b, T,p);
975 70 : a = b; b = c; dc = degpol(c);
976 70 : if (dc < 0) { set_avma(av); return pol_0(vT); }
977 :
978 70 : if (both_odd(da,db)) res = FpX_neg(res, p);
979 70 : if (!ZX_equal1(lb)) res = FpXQ_mul(res, FpXQ_powu(lb, da - dc, T, p), T, p);
980 70 : if (gc_needed(av,2))
981 : {
982 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_resultant (da = %ld)",da);
983 0 : gerepileall(av,3, &a,&b,&res);
984 : }
985 70 : da = db; /* = degpol(a) */
986 70 : db = dc; /* = degpol(b) */
987 : }
988 28 : res = FpXQ_mul(res, FpXQ_powu(gel(b,2), da, T, p), T, p);
989 28 : return gerepileupto(av, res);
990 : }
991 :
992 : GEN
993 63 : FpXQX_resultant(GEN x, GEN y, GEN T, GEN p)
994 : {
995 63 : pari_sp av = avma;
996 63 : long dx, dy, vT = get_FpX_var(T);
997 63 : GEN res = pol_1(vT);
998 63 : if (!signe(x) || !signe(y)) return pol_0(vT);
999 63 : if (lgefint(p) == 3)
1000 : {
1001 35 : pari_sp av = avma;
1002 : GEN Pl, Ql, Tl, R;
1003 35 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
1004 35 : R = FlxqX_resultant(Pl, Ql, Tl, pp);
1005 35 : return gerepileupto(av, Flx_to_ZX(R));
1006 : }
1007 :
1008 28 : dx = degpol(x); dy = degpol(y);
1009 28 : if (dx < dy)
1010 : {
1011 14 : swap(x,y);
1012 14 : if (both_odd(dx, dy))
1013 0 : res = Fp_neg(res, p);
1014 : }
1015 28 : while (lgpol(y) >= FpXQX_GCD_LIMIT)
1016 : {
1017 0 : if (lgpol(y)<=(lgpol(x)>>1))
1018 : {
1019 0 : GEN r = FpXQX_rem(x, y, T, p);
1020 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
1021 0 : GEN ly = FpX_red(gel(y,dy+2),p);
1022 0 : if (!ZX_equal1(ly)) res = FpXQ_mul(res, FpXQ_powu(ly, dx - dr, T, p), T, p);
1023 0 : if (both_odd(dx, dy))
1024 0 : res = Fp_neg(res, p);
1025 0 : x = y; y = r;
1026 : }
1027 0 : (void) FpXQX_halfres(x, y, T, p, &x, &y, &res);
1028 0 : if (gc_needed(av,2))
1029 : {
1030 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_resultant (y = %ld)",degpol(y));
1031 0 : gerepileall(av,3,&x,&y,&res);
1032 : }
1033 : }
1034 28 : return gerepileupto(av, FpXQ_mul(res, FpXQX_resultant_basecase(x, y, T, p), T, p));
1035 : }
1036 :
1037 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
1038 : GEN
1039 35 : FpXQX_disc(GEN P, GEN T, GEN p)
1040 : {
1041 35 : pari_sp av = avma;
1042 35 : GEN L, dP = FpXX_deriv(P, p), D = FpXQX_resultant(P, dP, T, p);
1043 : long dd;
1044 35 : if (!signe(D)) return pol_0(get_FpX_var(T));
1045 35 : dd = degpol(P) - 2 - degpol(dP); /* >= -1; > -1 iff p | deg(P) */
1046 35 : L = leading_coeff(P);
1047 35 : if (dd && !gequal1(L))
1048 0 : D = (dd == -1)? FpXQ_div(D,L,T,p): FpXQ_mul(D, FpXQ_powu(L, dd, T, p), T, p);
1049 35 : if (degpol(P) & 2) D = FpX_neg(D, p);
1050 35 : return gerepileupto(av, D);
1051 : }
1052 :
1053 : GEN
1054 396 : FpXQX_dotproduct(GEN x, GEN y, GEN T, GEN p)
1055 : {
1056 396 : long i, l = minss(lg(x), lg(y));
1057 : pari_sp av;
1058 : GEN c;
1059 396 : if (l == 2) return gen_0;
1060 396 : av = avma; c = gmul(gel(x,2),gel(y,2));
1061 1642 : for (i=3; i<l; i++) c = gadd(c, gmul(gel(x,i),gel(y,i)));
1062 396 : return gerepileupto(av, Fq_red(c,T,p));
1063 : }
1064 :
1065 : /***********************************************************************/
1066 : /** **/
1067 : /** Barrett reduction **/
1068 : /** **/
1069 : /***********************************************************************/
1070 :
1071 : /* Return new lgpol */
1072 : static long
1073 316761 : ZXX_lgrenormalizespec(GEN x, long lx)
1074 : {
1075 : long i;
1076 317271 : for (i = lx-1; i>=0; i--)
1077 317271 : if (signe(gel(x,i))) break;
1078 316761 : return i+1;
1079 : }
1080 :
1081 : static GEN
1082 2918 : FpXQX_invBarrett_basecase(GEN S, GEN T, GEN p)
1083 : {
1084 2918 : long i, l=lg(S)-1, lr = l-1, k;
1085 2918 : GEN r=cgetg(lr, t_POL); r[1]=S[1];
1086 2918 : gel(r,2) = gen_1;
1087 24894 : for (i=3; i<lr; i++)
1088 : {
1089 21976 : pari_sp av = avma;
1090 21976 : GEN u = gel(S,l-i+2);
1091 221103 : for (k=3; k<i; k++)
1092 199127 : u = Fq_add(u, Fq_mul(gel(S,l-i+k), gel(r,k), NULL, p), NULL, p);
1093 21976 : gel(r,i) = gerepileupto(av, Fq_red(Fq_neg(u, NULL, p), T, p));
1094 : }
1095 2918 : return FpXQX_renormalize(r,lr);
1096 : }
1097 :
1098 : INLINE GEN
1099 306360 : FpXX_recipspec(GEN x, long l, long n)
1100 : {
1101 306360 : return RgX_recipspec_shallow(x, l, n);
1102 : }
1103 :
1104 : static GEN
1105 537 : FpXQX_invBarrett_Newton(GEN S, GEN T, GEN p)
1106 : {
1107 537 : pari_sp av = avma;
1108 537 : long nold, lx, lz, lq, l = degpol(S), i, lQ;
1109 537 : GEN q, y, z, x = cgetg(l+2, t_POL) + 2;
1110 537 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1111 38210 : for (i=0;i<l;i++) gel(x,i) = gen_0;
1112 537 : q = RgX_recipspec_shallow(S+2,l+1,l+1); lQ = lgpol(q); q+=2;
1113 : /* We work on _spec_ FpX's, all the l[xzq] below are lgpol's */
1114 :
1115 : /* initialize */
1116 537 : gel(x,0) = Fq_inv(gel(q,0), T, p);
1117 537 : if (lQ>1) gel(q,1) = Fq_red(gel(q,1), T, p);
1118 537 : if (lQ>1 && signe(gel(q,1)))
1119 537 : {
1120 537 : GEN u = gel(q, 1);
1121 537 : if (!gequal1(gel(x,0))) u = Fq_mul(u, Fq_sqr(gel(x,0), T, p), T, p);
1122 537 : gel(x,1) = Fq_neg(u, T, p); lx = 2;
1123 : }
1124 : else
1125 0 : lx = 1;
1126 537 : nold = 1;
1127 4004 : for (; mask > 1; )
1128 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1129 3467 : long i, lnew, nnew = nold << 1;
1130 :
1131 3467 : if (mask & 1) nnew--;
1132 3467 : mask >>= 1;
1133 :
1134 3467 : lnew = nnew + 1;
1135 3467 : lq = ZXX_lgrenormalizespec(q, minss(lQ,lnew));
1136 3467 : z = FpXQX_mulspec(x, q, T, p, lx, lq); /* FIXME: high product */
1137 3467 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1138 3467 : z += 2;
1139 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1140 6934 : for (i = nold; i < lz; i++) if (signe(gel(z,i))) break;
1141 3467 : nold = nnew;
1142 3467 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1143 :
1144 : /* z + i represents (x*q - 1) / t^i */
1145 3467 : lz = ZXX_lgrenormalizespec (z+i, lz-i);
1146 3467 : z = FpXQX_mulspec(x, z+i, T, p, lx, lz); /* FIXME: low product */
1147 3467 : lz = lgpol(z); z += 2;
1148 3467 : if (lz > lnew-i) lz = ZXX_lgrenormalizespec(z, lnew-i);
1149 :
1150 3467 : lx = lz+ i;
1151 3467 : y = x + i; /* x -= z * t^i, in place */
1152 39529 : for (i = 0; i < lz; i++) gel(y,i) = Fq_neg(gel(z,i), T, p);
1153 : }
1154 537 : x -= 2; setlg(x, lx + 2); x[1] = S[1];
1155 537 : return gerepilecopy(av, x);
1156 : }
1157 :
1158 : GEN
1159 3476 : FpXQX_invBarrett(GEN S, GEN T, GEN p)
1160 : {
1161 3476 : pari_sp ltop = avma;
1162 3476 : long l = lg(S);
1163 : GEN r;
1164 3476 : if (l<5) return pol_0(varn(S));
1165 3455 : if (l<=FpXQX_INVBARRETT_LIMIT)
1166 : {
1167 2918 : GEN c = gel(S,l-1), ci=gen_1;
1168 2918 : if (!gequal1(c))
1169 : {
1170 1870 : ci = Fq_inv(c, T, p);
1171 1870 : S = FqX_Fq_mul(S, ci, T, p);
1172 1870 : r = FpXQX_invBarrett_basecase(S, T, p);
1173 1870 : r = FqX_Fq_mul(r, ci, T, p);
1174 : } else
1175 1048 : r = FpXQX_invBarrett_basecase(S, T, p);
1176 : }
1177 : else
1178 537 : r = FpXQX_invBarrett_Newton(S, T, p);
1179 3455 : return gerepileupto(ltop, r);
1180 : }
1181 :
1182 : GEN
1183 12685 : FpXQX_get_red(GEN S, GEN T, GEN p)
1184 : {
1185 12685 : if (typ(S)==t_POL && lg(S)>FpXQX_BARRETT_LIMIT)
1186 1174 : retmkvec2(FpXQX_invBarrett(S,T,p),S);
1187 11511 : return S;
1188 : }
1189 :
1190 : /* Compute x mod S where 2 <= degpol(S) <= l+1 <= 2*(degpol(S)-1)
1191 : * and mg is the Barrett inverse of S. */
1192 : static GEN
1193 153180 : FpXQX_divrem_Barrettspec(GEN x, long l, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
1194 : {
1195 : GEN q, r;
1196 153180 : long lt = degpol(S); /*We discard the leading term*/
1197 : long ld, lm, lT, lmg;
1198 153180 : ld = l-lt;
1199 153180 : lm = minss(ld, lgpol(mg));
1200 153180 : lT = ZXX_lgrenormalizespec(S+2,lt);
1201 153180 : lmg = ZXX_lgrenormalizespec(mg+2,lm);
1202 153180 : q = FpXX_recipspec(x+lt,ld,ld); /* q = rec(x) lq<=ld*/
1203 153180 : q = FpXQX_mulspec(q+2,mg+2,T,p,lgpol(q),lmg); /* q = rec(x) * mg lq<=ld+lm*/
1204 153180 : q = FpXX_recipspec(q+2,minss(ld,lgpol(q)),ld); /* q = rec (rec(x) * mg) lq<=ld*/
1205 153180 : if (!pr) return q;
1206 152221 : r = FpXQX_mulspec(q+2,S+2,T,p,lgpol(q),lT); /* r = q*pol lr<=ld+lt*/
1207 152221 : r = FpXX_subspec(x,r+2,p,lt,minss(lt,lgpol(r))); /* r = x - r lr<=lt */
1208 152221 : if (pr == ONLY_REM) return r;
1209 68338 : *pr = r; return q;
1210 : }
1211 :
1212 : static GEN
1213 85254 : FpXQX_divrem_Barrett(GEN x, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
1214 : {
1215 85254 : GEN q = NULL, r = FpXQX_red(x, T, p);
1216 85254 : long l = lgpol(r), lt = degpol(S), lm = 2*lt-1, v = varn(S);
1217 : long i;
1218 85254 : if (l <= lt)
1219 : {
1220 0 : if (pr == ONLY_REM) return r;
1221 0 : if (pr == ONLY_DIVIDES) return signe(r)? NULL: pol_0(v);
1222 0 : if (pr) *pr = r;
1223 0 : return pol_0(v);
1224 : }
1225 85254 : if (lt <= 1)
1226 21 : return FpXQX_divrem_basecase(r,S,T,p,pr);
1227 85233 : if (pr != ONLY_REM && l>lm)
1228 : {
1229 1350 : q = cgetg(l-lt+2, t_POL); q[1] = S[1];
1230 87536 : for (i=0;i<l-lt;i++) gel(q+2,i) = gen_0;
1231 : }
1232 153261 : while (l>lm)
1233 : {
1234 68028 : GEN zr, zq = FpXQX_divrem_Barrettspec(r+2+l-lm,lm,mg,S,T,p,&zr);
1235 68028 : long lz = lgpol(zr);
1236 68028 : if (pr != ONLY_REM)
1237 : {
1238 34904 : long lq = lgpol(zq);
1239 113973 : for(i=0; i<lq; i++) gel(q+2+l-lm,i) = gel(zq,2+i);
1240 : }
1241 284012 : for(i=0; i<lz; i++) gel(r+2+l-lm,i) = gel(zr,2+i);
1242 68028 : l = l-lm+lz;
1243 : }
1244 85233 : if (pr == ONLY_REM)
1245 : {
1246 83883 : if (l > lt)
1247 83883 : r = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,ONLY_REM);
1248 : else
1249 0 : r = FpXQX_renormalize(r, l+2);
1250 83883 : setvarn(r, v); return r;
1251 : }
1252 1350 : if (l > lt)
1253 : {
1254 1269 : GEN zq = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,pr ? &r: NULL);
1255 1269 : if (!q) q = zq;
1256 : else
1257 : {
1258 1269 : long lq = lgpol(zq);
1259 8050 : for(i=0; i<lq; i++) gel(q+2,i) = gel(zq,2+i);
1260 : }
1261 : }
1262 81 : else if (pr)
1263 81 : r = FpX_renormalize(r, l+2);
1264 1350 : setvarn(q, v); q = FpXQX_renormalize(q, lg(q));
1265 1350 : if (pr == ONLY_DIVIDES) return signe(r)? NULL: q;
1266 1350 : if (pr) { setvarn(r, v); *pr = r; }
1267 1350 : return q;
1268 : }
1269 :
1270 : GEN
1271 1069301 : FpXQX_divrem(GEN x, GEN S, GEN T, GEN p, GEN *pr)
1272 : {
1273 : GEN B, y;
1274 : long dy, dx, d;
1275 1069301 : if (pr==ONLY_REM) return FpXQX_rem(x, S, T, p);
1276 1069301 : y = get_FpXQX_red(S, &B);
1277 1069301 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1278 1069270 : if (lgefint(p) == 3)
1279 : {
1280 : GEN a, b, t, z;
1281 1050659 : pari_sp tetpil, av = avma;
1282 1050659 : ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1283 1050647 : z = FlxqX_divrem(a, b, t, pp, pr);
1284 1050651 : if (pr == ONLY_DIVIDES && !z) return gc_NULL(av);
1285 1050651 : tetpil=avma;
1286 1050651 : z = FlxX_to_ZXX(z);
1287 1050567 : if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
1288 1016590 : *pr = FlxX_to_ZXX(*pr);
1289 33977 : else return gerepile(av, tetpil, z);
1290 1016594 : gerepileallsp(av,tetpil,2, pr, &z);
1291 1016744 : return z;
1292 : }
1293 18611 : if (!B && d+3 < FpXQX_DIVREM_BARRETT_LIMIT)
1294 17240 : return FpXQX_divrem_basecase(x,y,T,p,pr);
1295 : else
1296 : {
1297 1371 : pari_sp av=avma;
1298 1371 : GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1299 1371 : GEN q = FpXQX_divrem_Barrett(x,mg,y,T,p,pr);
1300 1371 : if (!q) return gc_NULL(av);
1301 1371 : if (!pr || pr==ONLY_DIVIDES) return gerepilecopy(av, q);
1302 398 : return gc_all(av, 2, &q, pr);
1303 : }
1304 : }
1305 :
1306 : GEN
1307 263338 : FpXQX_rem(GEN x, GEN S, GEN T, GEN p)
1308 : {
1309 263338 : GEN B, y = get_FpXQX_red(S, &B);
1310 263338 : long dy = degpol(y), dx = degpol(x), d = dx-dy;
1311 263338 : if (d < 0) return FpXQX_red(x, T, p);
1312 242723 : if (lgefint(p) == 3)
1313 : {
1314 2736 : pari_sp av = avma;
1315 : GEN a, b, t, z;
1316 2736 : ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1317 2736 : z = FlxqX_rem(a, b, t, pp);
1318 2736 : z = FlxX_to_ZXX(z);
1319 2736 : return gerepileupto(av, z);
1320 : }
1321 239987 : if (!B && d+3 < FpXQX_REM_BARRETT_LIMIT)
1322 156104 : return FpXQX_divrem_basecase(x,y, T, p, ONLY_REM);
1323 : else
1324 : {
1325 83883 : pari_sp av=avma;
1326 83883 : GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1327 83883 : GEN r = FpXQX_divrem_Barrett(x, mg, y, T, p, ONLY_REM);
1328 83883 : return gerepileupto(av, r);
1329 : }
1330 : }
1331 :
1332 : /* x + y*z mod p */
1333 : INLINE GEN
1334 36478 : Fq_addmul(GEN x, GEN y, GEN z, GEN T, GEN p)
1335 : {
1336 : pari_sp av;
1337 36478 : if (!signe(y) || !signe(z)) return Fq_red(x, T, p);
1338 36478 : if (!signe(x)) return Fq_mul(z,y, T, p);
1339 36478 : av = avma;
1340 36478 : return gerepileupto(av, Fq_add(x, Fq_mul(y, z, T, p), T, p));
1341 : }
1342 :
1343 : GEN
1344 71939 : FpXQX_div_by_X_x(GEN a, GEN x, GEN T, GEN p, GEN *r)
1345 : {
1346 71939 : long l = lg(a), i;
1347 : GEN z;
1348 71939 : if (lgefint(p)==3)
1349 : {
1350 53700 : pari_sp av = avma;
1351 : GEN ap, xp, t, z;
1352 53700 : ulong pp = to_FlxqX(a, NULL, T, p, &ap, NULL, &t);
1353 53700 : xp = ZX_to_Flx(to_ZX(x, get_FpX_var(T)), pp);
1354 53700 : z = FlxX_to_ZXX(FlxqX_div_by_X_x(ap, xp, t, pp, r));
1355 53700 : if (!r) return gerepileupto(av, z);
1356 0 : *r = Flx_to_ZX(*r);
1357 0 : return gc_all(av, 2, &z, r);
1358 : }
1359 18239 : if (l <= 3)
1360 : {
1361 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1362 0 : return pol_0(varn(a));
1363 : }
1364 18239 : l--; z = cgetg(l, t_POL); z[1] = a[1];
1365 18239 : gel(z, l-1) = gel(a,l);
1366 54717 : for (i=l-2; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1367 36478 : gel(z, i) = Fq_addmul(gel(a,i+1), x, gel(z,i+1), T, p);
1368 18239 : if (r) *r = Fq_addmul(gel(a,2), x, gel(z,2), T, p);
1369 18239 : return z;
1370 : }
1371 :
1372 : struct _FpXQXQ {
1373 : GEN T, S;
1374 : GEN p;
1375 : };
1376 :
1377 118444 : static GEN _FpXQX_mul(void *data, GEN a,GEN b)
1378 : {
1379 118444 : struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1380 118444 : return FpXQX_mul(a,b,d->T,d->p);
1381 : }
1382 :
1383 1729 : static GEN _FpXQX_sqr(void *data, GEN a)
1384 : {
1385 1729 : struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1386 1729 : return FpXQX_sqr(a, d->T, d->p);
1387 : }
1388 :
1389 : GEN
1390 63 : FpXQX_powu(GEN x, ulong n, GEN T, GEN p)
1391 : {
1392 : struct _FpXQXQ D;
1393 63 : if (n==0) return pol_1(varn(x));
1394 63 : D.T = T; D.p = p;
1395 63 : return gen_powu(x, n, (void *)&D, _FpXQX_sqr, _FpXQX_mul);
1396 : }
1397 :
1398 : GEN
1399 16678 : FpXQXV_prod(GEN V, GEN T, GEN p)
1400 : {
1401 16678 : if (lgefint(p) == 3)
1402 : {
1403 0 : pari_sp av = avma;
1404 0 : ulong pp = p[2];
1405 0 : GEN Tl = ZXT_to_FlxT(T, pp);
1406 0 : GEN Vl = ZXXV_to_FlxXV(V, pp, get_FpX_var(T));
1407 0 : Tl = FlxqXV_prod(Vl, Tl, pp);
1408 0 : return gerepileupto(av, FlxX_to_ZXX(Tl));
1409 : }
1410 : else
1411 : {
1412 : struct _FpXQXQ d;
1413 16678 : d.T=T; d.p=p;
1414 16678 : return gen_product(V, (void*)&d, &_FpXQX_mul);
1415 : }
1416 : }
1417 :
1418 : static GEN
1419 9954 : _FpXQX_divrem(void * E, GEN x, GEN y, GEN *r)
1420 : {
1421 9954 : struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1422 9954 : return FpXQX_divrem(x, y, d->T, d->p, r);
1423 : }
1424 :
1425 : static GEN
1426 121791 : _FpXQX_add(void * E, GEN x, GEN y)
1427 : {
1428 121791 : struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1429 121791 : return FpXX_add(x, y, d->p);
1430 : }
1431 :
1432 : static GEN
1433 4638 : _FpXQX_sub(void * E, GEN x, GEN y) {
1434 4638 : struct _FpXQXQ *d = (struct _FpXQXQ*) E;
1435 4638 : return FpXX_sub(x,y, d->p);
1436 : }
1437 :
1438 : static struct bb_ring FpXQX_ring = { _FpXQX_add, _FpXQX_mul, _FpXQX_sqr };
1439 :
1440 : GEN
1441 623 : FpXQX_digits(GEN x, GEN B, GEN T, GEN p)
1442 : {
1443 623 : long d = degpol(B), n = (lgpol(x)+d-1)/d;
1444 : struct _FpXQXQ D;
1445 623 : D.T = T; D.p = p;
1446 623 : return gen_digits(x, B, n, (void *)&D, &FpXQX_ring, _FpXQX_divrem);
1447 : }
1448 :
1449 : GEN
1450 189 : FpXQXV_FpXQX_fromdigits(GEN x, GEN B, GEN T, GEN p)
1451 : {
1452 : struct _FpXQXQ D;
1453 189 : D.T = T; D.p = p;
1454 189 : return gen_fromdigits(x,B,(void *)&D, &FpXQX_ring);
1455 : }
1456 :
1457 : /* Q an FpXY (t_POL with FpX coeffs), evaluate at X = x */
1458 : GEN
1459 51552 : FpXY_evalx(GEN Q, GEN x, GEN p)
1460 : {
1461 51552 : long i, lb = lg(Q);
1462 : GEN z;
1463 51552 : z = cgetg(lb, t_POL); z[1] = Q[1];
1464 445045 : for (i=2; i<lb; i++)
1465 : {
1466 393493 : GEN q = gel(Q,i);
1467 393493 : gel(z,i) = typ(q) == t_INT? modii(q,p): FpX_eval(q, x, p);
1468 : }
1469 51552 : return FpX_renormalize(z, lb);
1470 : }
1471 : /* Q an FpXY, evaluate at Y = y */
1472 : GEN
1473 18799 : FpXY_evaly(GEN Q, GEN y, GEN p, long vx)
1474 : {
1475 18799 : pari_sp av = avma;
1476 18799 : long i, lb = lg(Q);
1477 : GEN z;
1478 18799 : if (!signe(Q)) return pol_0(vx);
1479 18771 : if (lb == 3 || !signe(y)) {
1480 84 : z = gel(Q, 2);
1481 84 : return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1482 : }
1483 18687 : z = gel(Q, lb-1);
1484 18687 : if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1485 242061 : for (i=lb-2; i>=2; i--) z = Fq_add(gel(Q,i), FpX_Fp_mul(z, y, p), NULL, p);
1486 18687 : return gerepileupto(av, z);
1487 : }
1488 : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
1489 : GEN
1490 13657 : FpXY_eval(GEN Q, GEN y, GEN x, GEN p)
1491 : {
1492 13657 : pari_sp av = avma;
1493 13657 : return gerepileuptoint(av, FpX_eval(FpXY_evalx(Q, x, p), y, p));
1494 : }
1495 :
1496 : GEN
1497 3852 : FpXY_FpXQV_evalx(GEN P, GEN x, GEN T, GEN p)
1498 : {
1499 3852 : long i, lP = lg(P);
1500 3852 : GEN res = cgetg(lP,t_POL);
1501 3852 : res[1] = P[1];
1502 61377 : for(i=2; i<lP; i++)
1503 115050 : gel(res,i) = typ(gel(P,i))==t_INT? icopy(gel(P,i)):
1504 57525 : FpX_FpXQV_eval(gel(P,i), x, T, p);
1505 3852 : return FlxX_renormalize(res, lP);
1506 : }
1507 :
1508 : GEN
1509 154 : FpXY_FpXQ_evalx(GEN P, GEN x, GEN T, GEN p)
1510 : {
1511 154 : pari_sp av = avma;
1512 154 : long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(P),1);
1513 154 : GEN xp = FpXQ_powers(x, n, T, p);
1514 154 : return gerepileupto(av, FpXY_FpXQV_evalx(P, xp, T, p));
1515 : }
1516 :
1517 : /*******************************************************************/
1518 : /* */
1519 : /* (Fp[X]/T(X))[Y] / S(Y) */
1520 : /* */
1521 : /*******************************************************************/
1522 :
1523 : /*Preliminary implementation to speed up FpX_ffisom*/
1524 : typedef struct {
1525 : GEN S, T, p;
1526 : } FpXYQQ_muldata;
1527 :
1528 : /* reduce x in Fp[X, Y] in the algebra Fp[X,Y]/ (S(X),T(Y)) */
1529 : static GEN
1530 476 : FpXYQQ_redswap(GEN x, GEN S, GEN T, GEN p)
1531 : {
1532 476 : pari_sp ltop=avma;
1533 476 : long n = get_FpX_degree(S);
1534 476 : long m = get_FpX_degree(T);
1535 476 : long v = get_FpX_var(T);
1536 476 : GEN V = RgXY_swap(x,m,v);
1537 476 : V = FpXQX_red(V,S,p);
1538 476 : V = RgXY_swap(V,n,v);
1539 476 : return gerepilecopy(ltop,V);
1540 : }
1541 : static GEN
1542 280 : FpXYQQ_sqr(void *data, GEN x)
1543 : {
1544 280 : FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1545 280 : return FpXYQQ_redswap(FpXQX_sqr(x, D->T, D->p),D->S,D->T,D->p);
1546 :
1547 : }
1548 : static GEN
1549 196 : FpXYQQ_mul(void *data, GEN x, GEN y)
1550 : {
1551 196 : FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1552 196 : return FpXYQQ_redswap(FpXQX_mul(x,y, D->T, D->p),D->S,D->T,D->p);
1553 : }
1554 :
1555 : /* x in Z[X,Y], S in Z[X] over Fq = Z[Y]/(p,T); compute lift(x^n mod (S,T,p)) */
1556 : GEN
1557 182 : FpXYQQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1558 : {
1559 182 : pari_sp av = avma;
1560 : FpXYQQ_muldata D;
1561 : GEN y;
1562 182 : if (lgefint(p) == 3)
1563 : {
1564 0 : ulong pp = to_FlxqX(x, NULL, T, p, &x, NULL, &T);
1565 0 : S = ZX_to_Flx(S, pp);
1566 0 : y = FlxX_to_ZXX( FlxYqq_pow(x, n, S, T, pp) );
1567 0 : y = gerepileupto(av, y);
1568 : }
1569 : else
1570 : {
1571 182 : D.S = S;
1572 182 : D.T = T;
1573 182 : D.p = p;
1574 182 : y = gen_pow(x, n, (void*)&D, &FpXYQQ_sqr, &FpXYQQ_mul);
1575 : }
1576 182 : return y;
1577 : }
1578 :
1579 : GEN
1580 49194 : FpXQXQ_mul(GEN x, GEN y, GEN S, GEN T, GEN p) {
1581 49194 : return FpXQX_rem(FpXQX_mul(x, y, T, p), S, T, p);
1582 : }
1583 :
1584 : GEN
1585 176117 : FpXQXQ_sqr(GEN x, GEN S, GEN T, GEN p) {
1586 176117 : return FpXQX_rem(FpXQX_sqr(x, T, p), S, T, p);
1587 : }
1588 :
1589 : /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
1590 : * return lift(1 / (x mod (p,pol))) */
1591 : GEN
1592 0 : FpXQXQ_invsafe(GEN x, GEN S, GEN T, GEN p)
1593 : {
1594 0 : GEN V, z = FpXQX_extgcd(get_FpXQX_mod(S), x, T, p, NULL, &V);
1595 0 : if (degpol(z)) return NULL;
1596 0 : z = gel(z,2);
1597 0 : z = typ(z)==t_INT ? Fp_invsafe(z,p) : FpXQ_invsafe(z,T,p);
1598 0 : if (!z) return NULL;
1599 0 : return typ(z)==t_INT ? FpXX_Fp_mul(V, z, p): FpXQX_FpXQ_mul(V, z, T, p);
1600 : }
1601 :
1602 : GEN
1603 0 : FpXQXQ_inv(GEN x, GEN S, GEN T,GEN p)
1604 : {
1605 0 : pari_sp av = avma;
1606 0 : GEN U = FpXQXQ_invsafe(x, S, T, p);
1607 0 : if (!U) pari_err_INV("FpXQXQ_inv",x);
1608 0 : return gerepileupto(av, U);
1609 : }
1610 :
1611 : GEN
1612 0 : FpXQXQ_div(GEN x,GEN y,GEN S, GEN T,GEN p)
1613 : {
1614 0 : pari_sp av = avma;
1615 0 : return gerepileupto(av, FpXQXQ_mul(x, FpXQXQ_inv(y,S,T,p),S,T,p));
1616 : }
1617 :
1618 : static GEN
1619 125668 : _FpXQXQ_cmul(void *data, GEN P, long a, GEN x) {
1620 125668 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1621 125668 : GEN y = gel(P,a+2);
1622 251313 : return typ(y)==t_INT ? FpXX_Fp_mul(x,y, d->p):
1623 125645 : FpXX_FpX_mul(x,y,d->p);
1624 : }
1625 : static GEN
1626 22956 : _FpXQXQ_red(void *data, GEN x) {
1627 22956 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1628 22956 : return FpXQX_red(x, d->T, d->p);
1629 : }
1630 : static GEN
1631 45926 : _FpXQXQ_mul(void *data, GEN x, GEN y) {
1632 45926 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1633 45926 : return FpXQXQ_mul(x,y, d->S,d->T, d->p);
1634 : }
1635 : static GEN
1636 176117 : _FpXQXQ_sqr(void *data, GEN x) {
1637 176117 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1638 176117 : return FpXQXQ_sqr(x, d->S,d->T, d->p);
1639 : }
1640 :
1641 : static GEN
1642 18963 : _FpXQXQ_one(void *data) {
1643 18963 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1644 18963 : return pol_1(get_FpXQX_var(d->S));
1645 : }
1646 :
1647 : static GEN
1648 132 : _FpXQXQ_zero(void *data) {
1649 132 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1650 132 : return pol_0(get_FpXQX_var(d->S));
1651 : }
1652 :
1653 : static struct bb_algebra FpXQXQ_algebra = { _FpXQXQ_red, _FpXQX_add,
1654 : _FpXQX_sub, _FpXQXQ_mul, _FpXQXQ_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1655 :
1656 : const struct bb_algebra *
1657 352 : get_FpXQXQ_algebra(void **E, GEN S, GEN T, GEN p)
1658 : {
1659 352 : GEN z = new_chunk(sizeof(struct _FpXQXQ));
1660 352 : struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1661 352 : e->T = FpX_get_red(T, p);
1662 352 : e->S = FpXQX_get_red(S, e->T, p);
1663 352 : e->p = p; *E = (void*)e;
1664 352 : return &FpXQXQ_algebra;
1665 : }
1666 :
1667 : static struct bb_algebra FpXQX_algebra = { _FpXQXQ_red, _FpXQX_add,
1668 : _FpXQX_sub, _FpXQX_mul, _FpXQX_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1669 :
1670 : const struct bb_algebra *
1671 0 : get_FpXQX_algebra(void **E, GEN T, GEN p, long v)
1672 : {
1673 0 : GEN z = new_chunk(sizeof(struct _FpXQXQ));
1674 0 : struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1675 0 : e->T = FpX_get_red(T, p);
1676 0 : e->S = pol_x(v);
1677 0 : e->p = p; *E = (void*)e;
1678 0 : return &FpXQX_algebra;
1679 : }
1680 :
1681 : /* x over Fq, return lift(x^n) mod S */
1682 : GEN
1683 1908 : FpXQXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1684 : {
1685 1908 : pari_sp ltop = avma;
1686 : GEN y;
1687 : struct _FpXQXQ D;
1688 1908 : long s = signe(n);
1689 1908 : if (!s) return pol_1(varn(x));
1690 1908 : if (is_pm1(n)) /* +/- 1 */
1691 0 : return (s < 0)? FpXQXQ_inv(x,S,T,p): ZXX_copy(x);
1692 1908 : if (lgefint(p) == 3)
1693 : {
1694 35 : ulong pp = to_FlxqX(x, S, T, p, &x, &S, &T);
1695 35 : GEN z = FlxqXQ_pow(x, n, S, T, pp);
1696 35 : y = FlxX_to_ZXX(z);
1697 35 : return gerepileupto(ltop, y);
1698 : }
1699 : else
1700 : {
1701 1873 : T = FpX_get_red(T, p);
1702 1873 : S = FpXQX_get_red(S, T, p);
1703 1873 : D.S = S; D.T = T; D.p = p;
1704 1873 : if (s < 0) x = FpXQXQ_inv(x,S,T,p);
1705 1873 : y = gen_pow_i(x, n, (void*)&D,&_FpXQXQ_sqr,&_FpXQXQ_mul);
1706 1873 : return gerepilecopy(ltop, y);
1707 : }
1708 : }
1709 :
1710 : /* generates the list of powers of x of degree 0,1,2,...,l*/
1711 : GEN
1712 1965 : FpXQXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
1713 : {
1714 : struct _FpXQXQ D;
1715 1965 : int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1716 1965 : T = FpX_get_red(T, p);
1717 1965 : S = FpXQX_get_red(S, T, p);
1718 1965 : D.S = S; D.T = T; D.p = p;
1719 1965 : return gen_powers(x, l, use_sqr, (void*)&D, &_FpXQXQ_sqr, &_FpXQXQ_mul,&_FpXQXQ_one);
1720 : }
1721 :
1722 : /* Let v a linear form, return the linear form z->v(tau*z)
1723 : that is, v*(M_tau) */
1724 :
1725 : INLINE GEN
1726 248 : FpXQX_recipspec(GEN x, long l, long n)
1727 : {
1728 248 : return RgX_recipspec_shallow(x, l, n);
1729 : }
1730 :
1731 : static GEN
1732 88 : FpXQXQ_transmul_init(GEN tau, GEN S, GEN T, GEN p)
1733 : {
1734 : GEN bht;
1735 88 : GEN h, Sp = get_FpXQX_red(S, &h);
1736 88 : long n = degpol(Sp), vT = varn(Sp);
1737 88 : GEN ft = FpXQX_recipspec(Sp+2, n+1, n+1);
1738 88 : GEN bt = FpXQX_recipspec(tau+2, lgpol(tau), n);
1739 88 : setvarn(ft, vT); setvarn(bt, vT);
1740 88 : if (h)
1741 16 : bht = FpXQXn_mul(bt, h, n-1, T, p);
1742 : else
1743 : {
1744 72 : GEN bh = FpXQX_div(FpXX_shift(tau, n-1), S, T, p);
1745 72 : bht = FpXQX_recipspec(bh+2, lgpol(bh), n-1);
1746 72 : setvarn(bht, vT);
1747 : }
1748 88 : return mkvec3(bt, bht, ft);
1749 : }
1750 :
1751 : static GEN
1752 200 : FpXQXQ_transmul(GEN tau, GEN a, long n, GEN T, GEN p)
1753 : {
1754 200 : pari_sp ltop = avma;
1755 : GEN t1, t2, t3, vec;
1756 200 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
1757 200 : if (signe(a)==0) return pol_0(varn(a));
1758 200 : t2 = FpXX_shift(FpXQX_mul(bt, a, T, p),1-n);
1759 200 : if (signe(bht)==0) return gerepilecopy(ltop, t2);
1760 114 : t1 = FpXX_shift(FpXQX_mul(ft, a, T, p),-n);
1761 114 : t3 = FpXQXn_mul(t1, bht, n-1, T, p);
1762 114 : vec = FpXX_sub(t2, FpXX_shift(t3, 1), p);
1763 114 : return gerepileupto(ltop, vec);
1764 : }
1765 :
1766 : static GEN
1767 44 : polxn_FpXX(long n, long v, long vT)
1768 : {
1769 44 : long i, a = n+2;
1770 44 : GEN p = cgetg(a+1, t_POL);
1771 44 : p[1] = evalsigne(1)|evalvarn(v);
1772 440 : for (i = 2; i < a; i++) gel(p,i) = pol_0(vT);
1773 44 : gel(p,a) = pol_1(vT); return p;
1774 : }
1775 :
1776 : GEN
1777 44 : FpXQXQ_minpoly(GEN x, GEN S, GEN T, GEN p)
1778 : {
1779 44 : pari_sp ltop = avma;
1780 : long vS, vT, n;
1781 : GEN v_x, g, tau;
1782 44 : vS = get_FpXQX_var(S);
1783 44 : vT = get_FpX_var(T);
1784 44 : n = get_FpXQX_degree(S);
1785 44 : g = pol_1(vS);
1786 44 : tau = pol_1(vS);
1787 44 : S = FpXQX_get_red(S, T, p);
1788 44 : v_x = FpXQXQ_powers(x, usqrt(2*n), S, T, p);
1789 88 : while(signe(tau) != 0)
1790 : {
1791 : long i, j, m, k1;
1792 : GEN M, v, tr;
1793 : GEN g_prime, c;
1794 44 : if (degpol(g) == n) { tau = pol_1(vS); g = pol_1(vS); }
1795 44 : v = random_FpXQX(n, vS, T, p);
1796 44 : tr = FpXQXQ_transmul_init(tau, S, T, p);
1797 44 : v = FpXQXQ_transmul(tr, v, n, T, p);
1798 44 : m = 2*(n-degpol(g));
1799 44 : k1 = usqrt(m);
1800 44 : tr = FpXQXQ_transmul_init(gel(v_x,k1+1), S, T, p);
1801 44 : c = cgetg(m+2,t_POL);
1802 44 : c[1] = evalsigne(1)|evalvarn(vS);
1803 200 : for (i=0; i<m; i+=k1)
1804 : {
1805 156 : long mj = minss(m-i, k1);
1806 552 : for (j=0; j<mj; j++)
1807 396 : gel(c,m+1-(i+j)) = FpXQX_dotproduct(v, gel(v_x,j+1), T, p);
1808 156 : v = FpXQXQ_transmul(tr, v, n, T, p);
1809 : }
1810 44 : c = FpXX_renormalize(c, m+2);
1811 : /* now c contains <v,x^i> , i = 0..m-1 */
1812 44 : M = FpXQX_halfgcd(polxn_FpXX(m, vS, vT), c, T, p);
1813 44 : g_prime = gmael(M, 2, 2);
1814 44 : if (degpol(g_prime) < 1) continue;
1815 44 : g = FpXQX_mul(g, g_prime, T, p);
1816 44 : tau = FpXQXQ_mul(tau, FpXQX_FpXQXQV_eval(g_prime, v_x, S, T, p), S, T, p);
1817 : }
1818 44 : g = FpXQX_normalize(g,T, p);
1819 44 : return gerepilecopy(ltop,g);
1820 : }
1821 :
1822 : GEN
1823 0 : FpXQXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
1824 : {
1825 0 : return RgXV_to_RgM(FpXQXQ_powers(y,m-1,S,T,p),n);
1826 : }
1827 :
1828 : GEN
1829 4191 : FpXQX_FpXQXQV_eval(GEN P, GEN V, GEN S, GEN T, GEN p)
1830 : {
1831 : struct _FpXQXQ D;
1832 4191 : T = FpX_get_red(T, p);
1833 4191 : S = FpXQX_get_red(S, T, p);
1834 4191 : D.S=S; D.T=T; D.p=p;
1835 4191 : return gen_bkeval_powers(P, degpol(P), V, (void*)&D, &FpXQXQ_algebra,
1836 : _FpXQXQ_cmul);
1837 : }
1838 :
1839 : GEN
1840 855 : FpXQX_FpXQXQ_eval(GEN Q, GEN x, GEN S, GEN T, GEN p)
1841 : {
1842 : struct _FpXQXQ D;
1843 855 : int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1844 855 : T = FpX_get_red(T, p);
1845 855 : S = FpXQX_get_red(S, T, p);
1846 855 : D.S=S; D.T=T; D.p=p;
1847 855 : return gen_bkeval(Q, degpol(Q), x, use_sqr, (void*)&D, &FpXQXQ_algebra,
1848 : _FpXQXQ_cmul);
1849 : }
1850 :
1851 : static GEN
1852 507 : FpXQXQ_autpow_sqr(void * E, GEN x)
1853 : {
1854 507 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1855 507 : GEN S = D->S, T = D->T, p = D->p;
1856 507 : GEN phi = gel(x,1), S1 = gel(x,2);
1857 507 : long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(S1)+1,1);
1858 507 : GEN V = FpXQ_powers(phi, n, T, p);
1859 507 : GEN phi2 = FpX_FpXQV_eval(phi, V, T, p);
1860 507 : GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1861 507 : GEN S2 = FpXQX_FpXQXQ_eval(Sphi, S1, S, T, p);
1862 507 : return mkvec2(phi2, S2);
1863 : }
1864 :
1865 : static GEN
1866 325 : FpXQXQ_autpow_mul(void * E, GEN x, GEN y)
1867 : {
1868 325 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1869 325 : GEN S = D->S, T = D->T, p = D->p;
1870 325 : GEN phi1 = gel(x,1), S1 = gel(x,2);
1871 325 : GEN phi2 = gel(y,1), S2 = gel(y,2);
1872 325 : long n = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+1, 1);
1873 325 : GEN V = FpXQ_powers(phi2, n, T, p);
1874 325 : GEN phi3 = FpX_FpXQV_eval(phi1, V, T, p);
1875 325 : GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1876 325 : GEN S3 = FpXQX_FpXQXQ_eval(Sphi, S2, S, T, p);
1877 325 : return mkvec2(phi3, S3);
1878 : }
1879 :
1880 : GEN
1881 493 : FpXQXQ_autpow(GEN aut, long n, GEN S, GEN T, GEN p)
1882 : {
1883 493 : pari_sp av = avma;
1884 : struct _FpXQXQ D;
1885 493 : T = FpX_get_red(T, p);
1886 493 : S = FpXQX_get_red(S, T, p);
1887 493 : D.S=S; D.T=T; D.p=p;
1888 493 : aut = gen_powu_i(aut,n,&D,FpXQXQ_autpow_sqr,FpXQXQ_autpow_mul);
1889 493 : return gerepilecopy(av, aut);
1890 : }
1891 :
1892 : static GEN
1893 1 : FpXQXQ_auttrace_mul(void *E, GEN x, GEN y)
1894 : {
1895 1 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1896 1 : GEN S = D->S, T = D->T;
1897 1 : GEN p = D->p;
1898 1 : GEN S1 = gel(x,1), a1 = gel(x,2);
1899 1 : GEN S2 = gel(y,1), a2 = gel(y,2);
1900 1 : long n = brent_kung_optpow(maxss(degpol(S1),degpol(a1)),2,1);
1901 1 : GEN V = FpXQXQ_powers(S2, n, S, T, p);
1902 1 : GEN S3 = FpXQX_FpXQXQV_eval(S1, V, S, T, p);
1903 1 : GEN aS = FpXQX_FpXQXQV_eval(a1, V, S, T, p);
1904 1 : GEN a3 = FpXX_add(aS, a2, p);
1905 1 : return mkvec2(S3, a3);
1906 : }
1907 :
1908 : static GEN
1909 1 : FpXQXQ_auttrace_sqr(void *E, GEN x)
1910 1 : { return FpXQXQ_auttrace_mul(E, x, x); }
1911 :
1912 : GEN
1913 8 : FpXQXQ_auttrace(GEN aut, long n, GEN S, GEN T, GEN p)
1914 : {
1915 8 : pari_sp av = avma;
1916 : struct _FpXQXQ D;
1917 8 : T = FpX_get_red(T, p);
1918 8 : S = FpXQX_get_red(S, T, p);
1919 8 : D.S=S; D.T=T; D.p=p;
1920 8 : aut = gen_powu_i(aut,n,&D,FpXQXQ_auttrace_sqr,FpXQXQ_auttrace_mul);
1921 8 : return gerepilecopy(av, aut);
1922 : }
1923 :
1924 : static GEN
1925 1433 : FpXQXQ_autsum_mul(void *E, GEN x, GEN y)
1926 : {
1927 1433 : struct _FpXQXQ *D = (struct _FpXQXQ *) E;
1928 1433 : GEN S = D->S, T = D->T, p = D->p;
1929 1433 : GEN phi1 = gel(x,1), S1 = gel(x,2), a1 = gel(x,3);
1930 1433 : GEN phi2 = gel(y,1), S2 = gel(y,2), a2 = gel(y,3);
1931 1433 : long n2 = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+lgpol(a1)+1, 1);
1932 1433 : GEN V2 = FpXQ_powers(phi2, n2, T, p);
1933 1433 : GEN phi3 = FpX_FpXQV_eval(phi1, V2, T, p);
1934 1433 : GEN Sphi = FpXY_FpXQV_evalx(S1, V2, T, p);
1935 1433 : GEN aphi = FpXY_FpXQV_evalx(a1, V2, T, p);
1936 1433 : long n = brent_kung_optpow(maxss(degpol(Sphi),degpol(aphi)),2,1);
1937 1433 : GEN V = FpXQXQ_powers(S2, n, S, T, p);
1938 1433 : GEN S3 = FpXQX_FpXQXQV_eval(Sphi, V, S, T, p);
1939 1433 : GEN aS = FpXQX_FpXQXQV_eval(aphi, V, S, T, p);
1940 1433 : GEN a3 = FpXQXQ_mul(aS, a2, S, T, p);
1941 1433 : return mkvec3(phi3, S3, a3);
1942 : }
1943 :
1944 : static GEN
1945 1275 : FpXQXQ_autsum_sqr(void * T, GEN x)
1946 1275 : { return FpXQXQ_autsum_mul(T,x,x); }
1947 :
1948 : GEN
1949 1233 : FpXQXQ_autsum(GEN aut, long n, GEN S, GEN T, GEN p)
1950 : {
1951 1233 : pari_sp av = avma;
1952 : struct _FpXQXQ D;
1953 1233 : T = FpX_get_red(T, p);
1954 1233 : S = FpXQX_get_red(S, T, p);
1955 1233 : D.S=S; D.T=T; D.p=p;
1956 1233 : aut = gen_powu_i(aut,n,&D,FpXQXQ_autsum_sqr,FpXQXQ_autsum_mul);
1957 1233 : return gerepilecopy(av, aut);
1958 : }
1959 :
1960 : GEN
1961 44573 : FpXQXn_mul(GEN x, GEN y, long n, GEN T, GEN p)
1962 : {
1963 44573 : pari_sp av = avma;
1964 : GEN z, kx, ky;
1965 : long d;
1966 44573 : if (ZXX_is_ZX(y) && ZXX_is_ZX(x))
1967 6937 : return FpXn_mul(x,y,n,p);
1968 37636 : d = get_FpX_degree(T);
1969 37636 : kx = RgXX_to_Kronecker(x, d);
1970 37636 : ky = RgXX_to_Kronecker(y, d);
1971 37636 : z = Kronecker_to_FpXQX(ZXn_mul(ky,kx,(2*d-1)*n), T, p);
1972 37636 : return gerepileupto(av, z);
1973 : }
1974 :
1975 : GEN
1976 0 : FpXQXn_sqr(GEN x, long n, GEN T, GEN p)
1977 : {
1978 0 : pari_sp av = avma;
1979 : GEN z, kx;
1980 : long d;
1981 0 : if (ZXX_is_ZX(x)) return ZXn_sqr(x, n);
1982 0 : d = get_FpX_degree(T);
1983 0 : kx= RgXX_to_Kronecker(x, d);
1984 0 : z = Kronecker_to_FpXQX(ZXn_sqr(kx, (2*d-1)*n), T, p);
1985 0 : return gerepileupto(av, z);
1986 : }
1987 :
1988 : /* (f*g) \/ x^n */
1989 : static GEN
1990 7399 : FpXQX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
1991 : {
1992 7399 : return FpXX_shift(FpXQX_mul(f,g,T, p),-n);
1993 : }
1994 :
1995 : static GEN
1996 4697 : FpXQXn_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
1997 : {
1998 4697 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
1999 4697 : return FpXX_add(FpXQX_mulhigh_i(fl, g, n2, T, p), FpXQXn_mul(fh, g, n - n2, T, p), p);
2000 : }
2001 :
2002 : /* Compute intformal(x^n*S)/x^(n+1) */
2003 : static GEN
2004 763 : FpXX_integXn(GEN x, long n, GEN p)
2005 : {
2006 763 : long i, lx = lg(x);
2007 : GEN y;
2008 763 : if (lx == 2) return ZXX_copy(x);
2009 763 : y = cgetg(lx, t_POL); y[1] = x[1];
2010 4317 : for (i=2; i<lx; i++)
2011 : {
2012 3554 : ulong j = n+i-1;
2013 3554 : GEN xi = gel(x,i);
2014 3554 : if (!signe(xi))
2015 0 : gel(y,i) = gen_0;
2016 : else
2017 3554 : gel(y,i) = typ(xi)==t_INT ? Fp_divu(xi, j, p)
2018 3554 : : FpX_divu(xi, j, p);
2019 : }
2020 763 : return ZXX_renormalize(y, lx);;
2021 : }
2022 :
2023 : /* Compute intformal(x^n*S)/x^(n+1) */
2024 : static GEN
2025 2702 : ZlXX_integXn(GEN x, long n, GEN p, ulong pp)
2026 : {
2027 2702 : long i, lx = lg(x);
2028 : GEN y;
2029 2702 : if (lx == 2) return ZXX_copy(x);
2030 2576 : if (!pp) return FpXX_integXn(x, n, p);
2031 1813 : y = cgetg(lx, t_POL); y[1] = x[1];
2032 6883 : for (i=2; i<lx; i++)
2033 : {
2034 5070 : GEN xi = gel(x,i);
2035 5070 : if (!signe(xi))
2036 14 : gel(y,i) = gen_0;
2037 : else
2038 : {
2039 : ulong j;
2040 5056 : long v = u_lvalrem(n+i-1, pp, &j);
2041 5056 : if (typ(xi)==t_INT)
2042 : {
2043 0 : if (v==0)
2044 0 : gel(y,i) = Fp_divu(xi, j, p);
2045 : else
2046 0 : gel(y,i) = Fp_divu(diviuexact(xi, upowuu(pp, v)), j, p);
2047 : } else
2048 : {
2049 5056 : if (v==0)
2050 5056 : gel(y,i) = FpX_divu(xi, j, p);
2051 : else
2052 0 : gel(y,i) = FpX_divu(ZX_divuexact(xi, upowuu(pp, v)), j, p);
2053 : }
2054 : }
2055 : }
2056 1813 : return ZXX_renormalize(y, lx);;
2057 : }
2058 :
2059 : GEN
2060 714 : ZlXQXn_expint(GEN h, long e, GEN T, GEN p, ulong pp)
2061 : {
2062 714 : pari_sp av = avma, av2;
2063 714 : long v = varn(h), n=1;
2064 714 : GEN f = pol_1(v), g = pol_1(v);
2065 714 : ulong mask = quadratic_prec_mask(e);
2066 714 : av2 = avma;
2067 2702 : for (;mask>1;)
2068 : {
2069 : GEN u, w;
2070 2702 : long n2 = n;
2071 2702 : n<<=1; if (mask & 1) n--;
2072 2702 : mask >>= 1;
2073 2702 : u = FpXQXn_mul(g, FpXQX_mulhigh_i(f, FpXXn_red(h, n2-1), n2-1, T, p), n-n2, T, p);
2074 2702 : u = FpXX_add(u, FpXX_shift(FpXXn_red(h, n-1), 1-n2), p);
2075 2702 : w = FpXQXn_mul(f, ZlXX_integXn(u, n2-1, p, pp), n-n2, T, p);
2076 2702 : f = FpXX_add(f, FpXX_shift(w, n2), p);
2077 2702 : if (mask<=1) break;
2078 1988 : u = FpXQXn_mul(g, FpXQXn_mulhigh(f, g, n2, n, T, p), n-n2, T, p);
2079 1988 : g = FpXX_sub(g, FpXX_shift(u, n2), p);
2080 1988 : if (gc_needed(av2,2))
2081 : {
2082 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_exp, e = %ld", n);
2083 0 : gerepileall(av2, 2, &f, &g);
2084 : }
2085 : }
2086 714 : return gerepileupto(av, f);
2087 : }
2088 :
2089 : GEN
2090 178 : FpXQXn_expint(GEN h, long e, GEN T, GEN p)
2091 178 : { return ZlXQXn_expint(h, e, T, p, 0); }
2092 :
2093 : GEN
2094 0 : FpXQXn_exp(GEN h, long e, GEN T, GEN p)
2095 : {
2096 0 : if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
2097 0 : pari_err_DOMAIN("FpXQXn_exp","valuation", "<", gen_1, h);
2098 0 : return FpXQXn_expint(FpXX_deriv(h, p), e, T, p);
2099 : }
2100 :
2101 : GEN
2102 721 : FpXQXn_div(GEN g, GEN f, long e, GEN T, GEN p)
2103 : {
2104 721 : pari_sp av = avma, av2;
2105 : ulong mask;
2106 : GEN W, a;
2107 721 : long v = varn(f), n = 1;
2108 :
2109 721 : if (!signe(f)) pari_err_INV("FpXXn_inv",f);
2110 721 : a = Fq_inv(gel(f,2), T, p);
2111 721 : if (e == 1 && !g) return scalarpol(a, v);
2112 721 : else if (e == 2 && !g)
2113 : {
2114 : GEN b;
2115 0 : if (degpol(f) <= 0) return scalarpol(a, v);
2116 0 : b = Fq_neg(gel(f,3),T,p);
2117 0 : if (signe(b)==0) return scalarpol(a, v);
2118 0 : if (!is_pm1(a)) b = Fq_mul(b, Fq_sqr(a, T, p), T, p);
2119 0 : W = deg1pol_shallow(b, a, v);
2120 0 : return gerepilecopy(av, W);
2121 : }
2122 721 : W = scalarpol_shallow(Fq_inv(gel(f,2), T, p),v);
2123 721 : mask = quadratic_prec_mask(e);
2124 721 : av2 = avma;
2125 3430 : for (;mask>1;)
2126 : {
2127 : GEN u, fr;
2128 2709 : long n2 = n;
2129 2709 : n<<=1; if (mask & 1) n--;
2130 2709 : mask >>= 1;
2131 2709 : fr = FpXXn_red(f, n);
2132 2709 : if (mask>1 || !g)
2133 : {
2134 2702 : u = FpXQXn_mul(W, FpXQXn_mulhigh(fr, W, n2, n, T, p), n-n2, T, p);
2135 2702 : W = FpXX_sub(W, FpXX_shift(u, n2), p);
2136 : }
2137 : else
2138 : {
2139 7 : GEN y = FpXQXn_mul(g, W, n, T, p), yt = FpXXn_red(y, n-n2);
2140 7 : u = FpXQXn_mul(yt, FpXQXn_mulhigh(fr, W, n2, n, T, p), n-n2, T, p);
2141 7 : W = FpXX_sub(y, FpXX_shift(u, n2), p);
2142 : }
2143 2709 : if (gc_needed(av2,2))
2144 : {
2145 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_inv, e = %ld", n);
2146 0 : W = gerepileupto(av2, W);
2147 : }
2148 : }
2149 721 : return gerepileupto(av, W);
2150 : }
2151 :
2152 : GEN
2153 714 : FpXQXn_inv(GEN f, long e, GEN T, GEN p)
2154 714 : { return FpXQXn_div(NULL, f, e, T, p); }
|
