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