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 1253143 : to_FlxqX(GEN P, GEN Q, GEN T, GEN p, GEN *pt_P, GEN *pt_Q, GEN *pt_T)
32 : {
33 1253143 : ulong pp = uel(p,2);
34 1253143 : long v = get_FpX_var(T);
35 1253129 : *pt_P = ZXX_to_FlxX(P, pp, v);
36 1253027 : if (pt_Q) *pt_Q = ZXX_to_FlxX(Q, pp, v);
37 1252947 : *pt_T = ZXT_to_FlxT(T, pp);
38 1253027 : return pp;
39 : }
40 :
41 : static GEN
42 126 : ZXX_copy(GEN a) { return gcopy(a); }
43 :
44 : GEN
45 39976 : FpXX_red(GEN z, GEN p)
46 : {
47 : GEN res;
48 39976 : long i, l = lg(z);
49 39976 : res = cgetg(l,t_POL); res[1] = z[1];
50 281829 : for (i=2; i<l; i++)
51 : {
52 241853 : GEN zi = gel(z,i), c;
53 241853 : if (typ(zi)==t_INT)
54 14924 : c = modii(zi,p);
55 : else
56 : {
57 226929 : pari_sp av = avma;
58 226929 : c = FpX_red(zi,p);
59 226929 : 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 241853 : gel(res,i) = c;
65 : }
66 39976 : 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 1334211 : get_FpXQX_red(GEN T, GEN *B)
220 : {
221 1334211 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
222 94923 : *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 1892748 : Kronecker_to_FpXQX(GEN Z, GEN T, GEN p)
238 : {
239 1892748 : long i,j,lx,l, N = (get_FpX_degree(T)<<1) + 1;
240 1892725 : GEN x, t = cgetg(N,t_POL), z = FpX_red(Z, p);
241 1891613 : t[1] = evalvarn(get_FpX_var(T));
242 1891701 : l = lg(z); lx = (l-2) / (N-2);
243 1891701 : x = cgetg(lx+3,t_POL);
244 1892362 : x[1] = z[1];
245 29941861 : for (i=2; i<lx+2; i++)
246 : {
247 238564910 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
248 28049606 : z += (N-2);
249 28049606 : gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
250 : }
251 1892255 : N = (l-2) % (N-2) + 2;
252 3386332 : for (j=2; j<N; j++) gel(t,j) = gel(z,j);
253 1892255 : gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
254 1892439 : return FpXQX_renormalize(x, i+1);
255 : }
256 :
257 : GEN
258 1939917 : FpXQX_red(GEN z, GEN T, GEN p)
259 : {
260 1939917 : long i, l = lg(z);
261 1939917 : GEN res = cgetg(l,t_POL); res[1] = z[1];
262 16308412 : for(i=2;i<l;i++)
263 14373058 : if (typ(gel(z,i)) == t_INT)
264 160149 : gel(res,i) = modii(gel(z,i),p);
265 : else
266 14212909 : gel(res,i) = FpXQ_red(gel(z,i),T,p);
267 1935354 : 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 3643190 : ZXX_is_ZX_spec(GEN a,long na)
372 : {
373 : long i;
374 3954668 : for(i=0;i<na;i++)
375 3894121 : if(typ(gel(a,i))!=t_INT) return 0;
376 60547 : 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 140819 : FpXX_FpX_mulspec(GEN P, GEN U, GEN p, long v, long lU)
384 : {
385 140819 : long i, lP =lg(P);
386 : GEN res;
387 140819 : res = cgetg(lP, t_POL); res[1] = P[1];
388 7696125 : for(i=2; i<lP; i++)
389 : {
390 7555306 : GEN Pi = gel(P,i);
391 7555306 : gel(res,i) = typ(Pi)==t_INT? FpX_Fp_mulspec(U, Pi, p, lU):
392 7540629 : FpX_mulspec(U, Pi+2, p, lU, lgpol(Pi));
393 7555306 : setvarn(gel(res,i),v);
394 : }
395 140819 : 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 15174 : FpXY_FpY_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
404 : {
405 15174 : pari_sp av = avma;
406 15174 : long v = get_FpX_var(T);
407 15174 : GEN z = RgXY_swapspec(x,get_FpX_degree(T)-1,v,lx);
408 15174 : z = FpXX_FpX_mulspec(z,y,p,v,ly);
409 15174 : z = RgXY_swapspec(z+2,lx+ly+3,v,lgpol(z));
410 15174 : return gerepilecopy(av,z);
411 : }
412 :
413 : static GEN
414 1698395 : FpXQX_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
415 : {
416 1698395 : pari_sp av = avma;
417 : GEN z, kx, ky;
418 : long n;
419 1698395 : if (ZXX_is_ZX_spec(y,ly))
420 : {
421 15687 : if (ZXX_is_ZX_spec(x,lx))
422 8271 : return FpX_mulspec(x,y,p,lx,ly);
423 : else
424 7416 : return FpXY_FpY_mulspec(x,y,T,p,lx,ly);
425 1682980 : } else if (ZXX_is_ZX_spec(x,lx))
426 7758 : return FpXY_FpY_mulspec(y,x,T,p,ly,lx);
427 1675427 : n = get_FpX_degree(T);
428 1675459 : kx = RgXX_to_Kronecker_spec(x, lx, n);
429 1675570 : ky = RgXX_to_Kronecker_spec(y, ly, n);
430 1675657 : z = Kronecker_to_FpXQX(ZX_mul(ky,kx), T, p);
431 1675273 : return gerepileupto(av, z);
432 : }
433 :
434 : GEN
435 1386105 : FpXQX_mul(GEN x, GEN y, GEN T, GEN p)
436 : {
437 1386105 : GEN z = FpXQX_mulspec(x+2,y+2,T,p,lgpol(x),lgpol(y));
438 1386710 : 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 553708 : FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p)
454 : {
455 : long i, lP;
456 : GEN res;
457 553708 : res = cgetg_copy(P, &lP); res[1] = P[1];
458 2067166 : for(i=2; i<lP; i++)
459 2749406 : gel(res,i) = typ(gel(P,i))==t_INT? FpX_Fp_mul(U, gel(P,i), p):
460 1235863 : FpXQ_mul(U, gel(P,i), T,p);
461 553482 : return FpXQX_renormalize(res,lP);
462 : }
463 :
464 : /* x and y in Z[Y][X]. Assume T irreducible mod p */
465 : static GEN
466 173392 : 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 173392 : if (!signe(y)) pari_err_INV("FpX_divrem",y);
473 173392 : vx=varn(x); dy=degpol(y); dx=degpol(x);
474 173392 : 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 173207 : lead = leading_coeff(y);
486 173207 : 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 171799 : av0 = avma; dz = dx-dy;
498 171799 : lead = gequal1(lead)? NULL: gclone(Fq_inv(lead,T,p));
499 171799 : set_avma(av0);
500 171799 : z = cgetg(dz+3,t_POL); z[1] = x[1];
501 171799 : x += 2; y += 2; z += 2;
502 178277 : for (dy1=dy-1; dy1>=0 && !signe(gel(y, dy1)); dy1--);
503 :
504 171799 : p1 = gel(x,dx); av = avma;
505 171799 : gel(z,dz) = lead? gerepileupto(av, Fq_mul(p1,lead, T, p)): gcopy(p1);
506 521176 : for (i=dx-1; i>=dy; i--)
507 : {
508 349377 : av=avma; p1=gel(x,i);
509 1164389 : for (j=i-dy1; j<=i && j<=dz; j++)
510 815012 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
511 349377 : if (lead) p1 = Fq_mul(p1, lead, NULL,p);
512 349377 : tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil,Fq_red(p1,T,p));
513 : }
514 171799 : if (!pr) { guncloneNULL(lead); return z-2; }
515 :
516 168640 : rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
517 180683 : for (sx=0; ; i--)
518 : {
519 180683 : p1 = gel(x,i);
520 692366 : for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
521 511683 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
522 180683 : 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 168640 : 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 168640 : lr=i+3; rem -= lr;
533 168640 : rem[0] = evaltyp(t_POL) | _evallg(lr);
534 168640 : rem[1] = z[-1];
535 168640 : p1 = gerepile((pari_sp)rem,tetpil,p1);
536 168640 : rem += 2; gel(rem,i) = p1;
537 1479870 : for (i--; i>=0; i--)
538 : {
539 1311230 : av=avma; p1 = gel(x,i);
540 4109454 : for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
541 2798224 : p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j), NULL,p), NULL,p);
542 1311230 : tetpil=avma; gel(rem,i) = gerepile(av,tetpil, Fq_red(p1, T, p));
543 : }
544 168640 : rem -= 2;
545 168640 : guncloneNULL(lead);
546 168640 : if (!sx) (void)FpXQX_renormalize(rem, lr);
547 168640 : if (pr == ONLY_REM) return gerepileupto(av0,rem);
548 13573 : *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 : if (res)
672 0 : gerepileall(av, 8, &a,&b,&u1,&v1,&u,&v,&res->res,&res->lc);
673 : else
674 0 : gerepileall(av, 6, &a,&b,&u1,&v1,&u,&v);
675 : }
676 : }
677 275 : M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
678 0 : return res ? gc_all(av, 5, &M, pa, pb, &res->res, &res->lc)
679 275 : : gc_all(av, 3, &M, pa, pb);
680 : }
681 :
682 : static GEN FpXQX_halfres_i(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res);
683 :
684 : static GEN
685 215 : FpXQX_halfres_split(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res)
686 : {
687 215 : pari_sp av = avma;
688 : GEN Q, R, S, V1, V2;
689 : GEN x1, y1, r, q;
690 215 : long l = lgpol(x), n = l>>1, k, vT = get_FpX_var(T);
691 215 : if (lgpol(y) <= n)
692 0 : { *a = RgX_copy(x); *b = RgX_copy(y); return matid2_FpXQXM(varn(x)); }
693 215 : if (res)
694 : {
695 0 : res->lc = to_ZX(leading_coeff(y), vT);
696 0 : res->deg0 -= n;
697 0 : res->deg1 -= n;
698 0 : res->off += n;
699 : }
700 215 : R = FpXQX_halfres_i(FpXX_shift(x,-n),FpXX_shift(y,-n), T, p, a, b, res);
701 215 : if (res)
702 : {
703 0 : res->off -= n;
704 0 : res->deg0 += n;
705 0 : res->deg1 += n;
706 : }
707 215 : V1 = FpXQXM_FpXQX_mul2(R, Flxn_red(x,n), Flxn_red(y,n), T, p);
708 215 : x1 = FpXX_add(FpXX_shift(*a,n), gel(V1,1), p);
709 215 : y1 = FpXX_add(FpXX_shift(*b,n), gel(V1,2), p);
710 215 : if (lgpol(y1) <= n)
711 : {
712 54 : *a = x1; *b = y1;
713 0 : return res ? gc_all(av, 5, &R, a, b, &res->res, &res->lc)
714 54 : : gc_all(av, 3, &R, a, b);
715 : }
716 161 : k = 2*n-degpol(y1);
717 161 : q = FpXQX_divrem(x1, y1, T, p, &r);
718 161 : if (res)
719 : {
720 0 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
721 0 : if (dy1 < degpol(y))
722 0 : FpXQX_halfres_update(res->deg0, res->deg1, dy1, T, p, res);
723 0 : res->lc = to_ZX(leading_coeff(y1), vT);
724 0 : res->deg0 = dx1;
725 0 : res->deg1 = dy1;
726 0 : if (dr >= n)
727 : {
728 0 : FpXQX_halfres_update(dx1, dy1, dr, T, p, res);
729 0 : res->deg0 = dy1;
730 0 : res->deg1 = dr;
731 : }
732 0 : res->deg0 -= k;
733 0 : res->deg1 -= k;
734 0 : res->off += k;
735 : }
736 161 : S = FpXQX_halfres_i(FpXX_shift(y1,-k), FpXX_shift(r,-k), T, p, a, b, res);
737 161 : if (res)
738 : {
739 0 : res->deg0 += k;
740 0 : res->deg1 += k;
741 0 : res->off -= k;
742 : }
743 161 : Q = FpXQXM_mul2(S,FpXQX_FpXQXM_qmul(q, R, T, p), T, p);
744 161 : V2 = FpXQXM_FpXQX_mul2(S, FpXXn_red(y1,k), FpXXn_red(r,k), T, p);
745 161 : *a = FpXX_add(FpXX_shift(*a,k), gel(V2,1), p);
746 161 : *b = FpXX_add(FpXX_shift(*b,k), gel(V2,2), p);
747 0 : return res ? gc_all(av, 5, &Q, a, b, &res->res, &res->lc)
748 161 : : gc_all(av, 3, &Q, a, b);
749 : }
750 :
751 : static GEN
752 490 : FpXQX_halfres_i(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, struct FpXQX_res *res)
753 : {
754 490 : if (lgpol(x) < FpXQX_HALFGCD_LIMIT)
755 275 : return FpXQX_halfres_basecase(x, y, T, p, a, b, res);
756 215 : return FpXQX_halfres_split(x, y, T, p, a, b, res);
757 : }
758 :
759 : static GEN
760 114 : FpXQX_halfgcd_all_i(GEN x, GEN y, GEN T, GEN p, GEN *pa, GEN *pb)
761 : {
762 : GEN a, b;
763 114 : GEN R = FpXQX_halfres_i(x, y, T, p, &a, &b, NULL);
764 114 : if (pa) *pa = a;
765 114 : if (pb) *pb = b;
766 114 : return R;
767 : }
768 :
769 : /* Return M in GL_2(Fp[X]/(T)[Y]) such that:
770 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
771 : */
772 :
773 : GEN
774 114 : FpXQX_halfgcd_all(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b)
775 : {
776 114 : pari_sp av = avma;
777 : GEN R,q,r;
778 114 : if (lgefint(p)==3)
779 : {
780 0 : ulong pp = to_FlxqX(x, y, T, p, &x, &y, &T);
781 0 : R = FlxXM_to_ZXXM(FlxqX_halfgcd(x, y, T, pp));
782 0 : if (a) *a = Flx_to_ZX(*a);
783 0 : if (b) *b = Flx_to_ZX(*b);
784 0 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
785 : }
786 114 : if (!signe(x))
787 : {
788 0 : if (a) *a = RgX_copy(y);
789 0 : if (b) *b = RgX_copy(x);
790 0 : return matJ2_FpXQXM(varn(x));
791 : }
792 114 : if (degpol(y)<degpol(x)) return FpXQX_halfgcd_all_i(x, y, T, p, a, b);
793 26 : q = FpXQX_divrem(y, x, T, p, &r);
794 26 : R = FpXQX_halfgcd_all_i(x, r, T, p, a, b);
795 26 : gcoeff(R,1,1) = FpXX_sub(gcoeff(R,1,1),
796 26 : FpXQX_mul(q, gcoeff(R,1,2), T, p), p);
797 26 : gcoeff(R,2,1) = FpXX_sub(gcoeff(R,2,1),
798 26 : FpXQX_mul(q, gcoeff(R,2,2), T, p), p);
799 26 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
800 : }
801 :
802 : GEN
803 44 : FpXQX_halfgcd(GEN x, GEN y, GEN T, GEN p)
804 44 : { return FpXQX_halfgcd_all(x, y, T, p, NULL, NULL); }
805 :
806 : static GEN
807 3896 : FpXQX_gcd_basecase(GEN a, GEN b, GEN T, GEN p)
808 : {
809 3896 : pari_sp av = avma, av0=avma;
810 39080 : while (signe(b))
811 : {
812 : GEN c;
813 35184 : if (gc_needed(av0,2))
814 : {
815 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_gcd (d = %ld)",degpol(b));
816 0 : gerepileall(av0,2, &a,&b);
817 : }
818 35184 : av = avma; c = FpXQX_rem(a, b, T, p); a=b; b=c;
819 : }
820 3896 : return gc_const(av, a);
821 : }
822 :
823 : GEN
824 14755 : FpXQX_gcd(GEN x, GEN y, GEN T, GEN p)
825 : {
826 14755 : pari_sp av = avma;
827 14755 : if (lgefint(p) == 3)
828 : {
829 : GEN Pl, Ql, Tl, U;
830 10750 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
831 10750 : U = FlxqX_gcd(Pl, Ql, Tl, pp);
832 10750 : return gerepileupto(av, FlxX_to_ZXX(U));
833 : }
834 4005 : x = FpXQX_red(x, T, p);
835 4005 : y = FpXQX_red(y, T, p);
836 4005 : if (!signe(x)) return gerepileupto(av, y);
837 3959 : while (lgpol(y)>=FpXQX_GCD_LIMIT)
838 : {
839 63 : if (lgpol(y)<=(lgpol(x)>>1))
840 : {
841 0 : GEN r = FpXQX_rem(x, y, T, p);
842 0 : x = y; y = r;
843 : }
844 63 : (void) FpXQX_halfgcd_all(x,y, T, p, &x, &y);
845 63 : if (gc_needed(av,2))
846 : {
847 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_gcd (y = %ld)",degpol(y));
848 0 : gerepileall(av,2,&x,&y);
849 : }
850 : }
851 3896 : return gerepileupto(av, FpXQX_gcd_basecase(x, y, T, p));
852 : }
853 :
854 : static GEN
855 0 : FpXQX_extgcd_basecase(GEN a, GEN b, GEN T, GEN p, GEN *ptu, GEN *ptv)
856 : {
857 0 : pari_sp av=avma;
858 : GEN u,v,d,d1,v1;
859 0 : long vx = varn(a);
860 0 : d = a; d1 = b;
861 0 : v = pol_0(vx); v1 = pol_1(vx);
862 0 : while (signe(d1))
863 : {
864 0 : GEN r, q = FpXQX_divrem(d, d1, T, p, &r);
865 0 : v = FpXX_sub(v,FpXQX_mul(q,v1,T, p),p);
866 0 : u=v; v=v1; v1=u;
867 0 : u=r; d=d1; d1=u;
868 0 : if (gc_needed(av,2))
869 : {
870 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_extgcd (d = %ld)",degpol(d));
871 0 : gerepileall(av,5, &d,&d1,&u,&v,&v1);
872 : }
873 : }
874 0 : if (ptu) *ptu = FpXQX_div(FpXX_sub(d,FpXQX_mul(b,v, T, p), p), a, T, p);
875 0 : *ptv = v; return d;
876 : }
877 :
878 : static GEN
879 0 : FpXQX_extgcd_halfgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
880 : {
881 : GEN u,v;
882 0 : GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
883 0 : long i, n = 0, vs = varn(x);
884 0 : while (lgpol(y) >= FpXQX_EXTGCD_LIMIT)
885 : {
886 0 : if (lgpol(y)<=(lgpol(x)>>1))
887 : {
888 0 : GEN r, q = FpXQX_divrem(x, y, T, p, &r);
889 0 : x = y; y = r;
890 0 : gel(V,++n) = mkmat22(pol_0(vs),pol_1(vs),pol_1(vs),FpXX_neg(q,p));
891 : } else
892 0 : gel(V,++n) = FpXQX_halfgcd_all(x, y, T, p, &x, &y);
893 : }
894 0 : y = FpXQX_extgcd_basecase(x,y, T, p, &u,&v);
895 0 : for (i = n; i>1; i--)
896 : {
897 0 : GEN R = gel(V,i);
898 0 : GEN u1 = FpXQX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), T, p);
899 0 : GEN v1 = FpXQX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), T, p);
900 0 : u = u1; v = v1;
901 : }
902 : {
903 0 : GEN R = gel(V,1);
904 0 : if (ptu)
905 0 : *ptu = FpXQX_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), T, p);
906 0 : *ptv = FpXQX_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), T, p);
907 : }
908 0 : return y;
909 : }
910 :
911 : /* x and y in Z[Y][X], return lift(gcd(x mod T,p, y mod T,p)). Set u and v st
912 : * ux + vy = gcd (mod T,p) */
913 : GEN
914 133795 : FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
915 : {
916 133795 : pari_sp av = avma;
917 : GEN d;
918 133795 : if (lgefint(p) == 3)
919 : {
920 : GEN Pl, Ql, Tl, Dl;
921 133795 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
922 133800 : Dl = FlxqX_extgcd(Pl, Ql, Tl, pp, ptu, ptv);
923 133804 : if (ptu) *ptu = FlxX_to_ZXX(*ptu);
924 133799 : *ptv = FlxX_to_ZXX(*ptv);
925 133799 : d = FlxX_to_ZXX(Dl);
926 : }
927 : else
928 : {
929 0 : x = FpXQX_red(x, T, p);
930 0 : y = FpXQX_red(y, T, p);
931 0 : if (lgpol(y)>=FpXQX_EXTGCD_LIMIT)
932 0 : d = FpXQX_extgcd_halfgcd(x, y, T, p, ptu, ptv);
933 : else
934 0 : d = FpXQX_extgcd_basecase(x, y, T, p, ptu, ptv);
935 : }
936 133798 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
937 : }
938 :
939 : static GEN
940 0 : FpXQX_halfres(GEN x, GEN y, GEN T, GEN p, GEN *a, GEN *b, GEN *r)
941 : {
942 : struct FpXQX_res res;
943 : GEN V;
944 0 : long dB, vT = get_FpX_var(T);
945 :
946 0 : res.res = *r;
947 0 : res.lc = to_ZX(leading_coeff(y),vT);
948 0 : res.deg0 = degpol(x);
949 0 : res.deg1 = degpol(y);
950 0 : res.off = 0;
951 0 : V = FpXQX_halfres_i(x, y, T, p, a, b, &res);
952 0 : dB = degpol(*b);
953 0 : if (dB < degpol(y))
954 0 : FpXQX_halfres_update(res.deg0, res.deg1, dB, T, p, &res);
955 0 : *r = res.res;
956 0 : return V;
957 : }
958 :
959 : /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab), with R=A%B, a=deg(A) ...*/
960 : static GEN
961 28 : FpXQX_resultant_basecase(GEN a, GEN b, GEN T, GEN p)
962 : {
963 28 : pari_sp av = avma;
964 28 : long vT = get_FpX_var(T), da,db,dc;
965 28 : GEN c,lb, res = pol_1(vT);
966 :
967 28 : if (!signe(a) || !signe(b)) return pol_0(vT);
968 :
969 28 : da = degpol(a);
970 28 : db = degpol(b);
971 28 : if (db > da)
972 : {
973 0 : swapspec(a,b, da,db);
974 0 : if (both_odd(da,db)) res = FpX_neg(res, p);
975 : }
976 28 : if (!da) return pol_1(vT); /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
977 98 : while (db)
978 : {
979 70 : lb = to_ZX(gel(b,db+2),vT);
980 70 : c = FpXQX_rem(a,b, T,p);
981 70 : a = b; b = c; dc = degpol(c);
982 70 : if (dc < 0) { set_avma(av); return pol_0(vT); }
983 :
984 70 : if (both_odd(da,db)) res = FpX_neg(res, p);
985 70 : if (!ZX_equal1(lb)) res = FpXQ_mul(res, FpXQ_powu(lb, da - dc, T, p), T, p);
986 70 : if (gc_needed(av,2))
987 : {
988 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_resultant (da = %ld)",da);
989 0 : gerepileall(av,3, &a,&b,&res);
990 : }
991 70 : da = db; /* = degpol(a) */
992 70 : db = dc; /* = degpol(b) */
993 : }
994 28 : res = FpXQ_mul(res, FpXQ_powu(gel(b,2), da, T, p), T, p);
995 28 : return gerepileupto(av, res);
996 : }
997 :
998 : GEN
999 63 : FpXQX_resultant(GEN x, GEN y, GEN T, GEN p)
1000 : {
1001 63 : pari_sp av = avma;
1002 63 : long dx, dy, vT = get_FpX_var(T);
1003 63 : GEN res = pol_1(vT);
1004 63 : if (!signe(x) || !signe(y)) return pol_0(vT);
1005 63 : if (lgefint(p) == 3)
1006 : {
1007 35 : pari_sp av = avma;
1008 : GEN Pl, Ql, Tl, R;
1009 35 : ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
1010 35 : R = FlxqX_resultant(Pl, Ql, Tl, pp);
1011 35 : return gerepileupto(av, Flx_to_ZX(R));
1012 : }
1013 :
1014 28 : dx = degpol(x); dy = degpol(y);
1015 28 : if (dx < dy)
1016 : {
1017 14 : swap(x,y);
1018 14 : if (both_odd(dx, dy))
1019 0 : res = Fp_neg(res, p);
1020 : }
1021 28 : while (lgpol(y) >= FpXQX_GCD_LIMIT)
1022 : {
1023 0 : if (lgpol(y)<=(lgpol(x)>>1))
1024 : {
1025 0 : GEN r = FpXQX_rem(x, y, T, p);
1026 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
1027 0 : GEN ly = FpX_red(gel(y,dy+2),p);
1028 0 : if (!ZX_equal1(ly)) res = FpXQ_mul(res, FpXQ_powu(ly, dx - dr, T, p), T, p);
1029 0 : if (both_odd(dx, dy))
1030 0 : res = Fp_neg(res, p);
1031 0 : x = y; y = r;
1032 : }
1033 0 : (void) FpXQX_halfres(x, y, T, p, &x, &y, &res);
1034 0 : if (gc_needed(av,2))
1035 : {
1036 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_resultant (y = %ld)",degpol(y));
1037 0 : gerepileall(av,3,&x,&y,&res);
1038 : }
1039 : }
1040 28 : return gerepileupto(av, FpXQ_mul(res, FpXQX_resultant_basecase(x, y, T, p), T, p));
1041 : }
1042 :
1043 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
1044 : GEN
1045 35 : FpXQX_disc(GEN P, GEN T, GEN p)
1046 : {
1047 35 : pari_sp av = avma;
1048 35 : GEN L, dP = FpXX_deriv(P, p), D = FpXQX_resultant(P, dP, T, p);
1049 : long dd;
1050 35 : if (!signe(D)) return pol_0(get_FpX_var(T));
1051 35 : dd = degpol(P) - 2 - degpol(dP); /* >= -1; > -1 iff p | deg(P) */
1052 35 : L = leading_coeff(P);
1053 35 : if (dd && !gequal1(L))
1054 0 : D = (dd == -1)? FpXQ_div(D,L,T,p): FpXQ_mul(D, FpXQ_powu(L, dd, T, p), T, p);
1055 35 : if (degpol(P) & 2) D = FpX_neg(D, p);
1056 35 : return gerepileupto(av, D);
1057 : }
1058 :
1059 : GEN
1060 396 : FpXQX_dotproduct(GEN x, GEN y, GEN T, GEN p)
1061 : {
1062 396 : long i, l = minss(lg(x), lg(y));
1063 : pari_sp av;
1064 : GEN c;
1065 396 : if (l == 2) return gen_0;
1066 396 : av = avma; c = gmul(gel(x,2),gel(y,2));
1067 1642 : for (i=3; i<l; i++) c = gadd(c, gmul(gel(x,i),gel(y,i)));
1068 396 : return gerepileupto(av, Fq_red(c,T,p));
1069 : }
1070 :
1071 : /***********************************************************************/
1072 : /** **/
1073 : /** Barrett reduction **/
1074 : /** **/
1075 : /***********************************************************************/
1076 :
1077 : /* Return new lgpol */
1078 : static long
1079 316761 : ZXX_lgrenormalizespec(GEN x, long lx)
1080 : {
1081 : long i;
1082 317271 : for (i = lx-1; i>=0; i--)
1083 317271 : if (signe(gel(x,i))) break;
1084 316761 : return i+1;
1085 : }
1086 :
1087 : static GEN
1088 2918 : FpXQX_invBarrett_basecase(GEN S, GEN T, GEN p)
1089 : {
1090 2918 : long i, l=lg(S)-1, lr = l-1, k;
1091 2918 : GEN r=cgetg(lr, t_POL); r[1]=S[1];
1092 2918 : gel(r,2) = gen_1;
1093 24894 : for (i=3; i<lr; i++)
1094 : {
1095 21976 : pari_sp av = avma;
1096 21976 : GEN u = gel(S,l-i+2);
1097 221103 : for (k=3; k<i; k++)
1098 199127 : u = Fq_add(u, Fq_mul(gel(S,l-i+k), gel(r,k), NULL, p), NULL, p);
1099 21976 : gel(r,i) = gerepileupto(av, Fq_red(Fq_neg(u, NULL, p), T, p));
1100 : }
1101 2918 : return FpXQX_renormalize(r,lr);
1102 : }
1103 :
1104 : INLINE GEN
1105 306360 : FpXX_recipspec(GEN x, long l, long n)
1106 : {
1107 306360 : return RgX_recipspec_shallow(x, l, n);
1108 : }
1109 :
1110 : static GEN
1111 537 : FpXQX_invBarrett_Newton(GEN S, GEN T, GEN p)
1112 : {
1113 537 : pari_sp av = avma;
1114 537 : long nold, lx, lz, lq, l = degpol(S), i, lQ;
1115 537 : GEN q, y, z, x = cgetg(l+2, t_POL) + 2;
1116 537 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1117 38210 : for (i=0;i<l;i++) gel(x,i) = gen_0;
1118 537 : q = RgX_recipspec_shallow(S+2,l+1,l+1); lQ = lgpol(q); q+=2;
1119 : /* We work on _spec_ FpX's, all the l[xzq] below are lgpol's */
1120 :
1121 : /* initialize */
1122 537 : gel(x,0) = Fq_inv(gel(q,0), T, p);
1123 537 : if (lQ>1) gel(q,1) = Fq_red(gel(q,1), T, p);
1124 537 : if (lQ>1 && signe(gel(q,1)))
1125 537 : {
1126 537 : GEN u = gel(q, 1);
1127 537 : if (!gequal1(gel(x,0))) u = Fq_mul(u, Fq_sqr(gel(x,0), T, p), T, p);
1128 537 : gel(x,1) = Fq_neg(u, T, p); lx = 2;
1129 : }
1130 : else
1131 0 : lx = 1;
1132 537 : nold = 1;
1133 4004 : for (; mask > 1; )
1134 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1135 3467 : long i, lnew, nnew = nold << 1;
1136 :
1137 3467 : if (mask & 1) nnew--;
1138 3467 : mask >>= 1;
1139 :
1140 3467 : lnew = nnew + 1;
1141 3467 : lq = ZXX_lgrenormalizespec(q, minss(lQ,lnew));
1142 3467 : z = FpXQX_mulspec(x, q, T, p, lx, lq); /* FIXME: high product */
1143 3467 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1144 3467 : z += 2;
1145 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1146 6934 : for (i = nold; i < lz; i++) if (signe(gel(z,i))) break;
1147 3467 : nold = nnew;
1148 3467 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1149 :
1150 : /* z + i represents (x*q - 1) / t^i */
1151 3467 : lz = ZXX_lgrenormalizespec (z+i, lz-i);
1152 3467 : z = FpXQX_mulspec(x, z+i, T, p, lx, lz); /* FIXME: low product */
1153 3467 : lz = lgpol(z); z += 2;
1154 3467 : if (lz > lnew-i) lz = ZXX_lgrenormalizespec(z, lnew-i);
1155 :
1156 3467 : lx = lz+ i;
1157 3467 : y = x + i; /* x -= z * t^i, in place */
1158 39529 : for (i = 0; i < lz; i++) gel(y,i) = Fq_neg(gel(z,i), T, p);
1159 : }
1160 537 : x -= 2; setlg(x, lx + 2); x[1] = S[1];
1161 537 : return gerepilecopy(av, x);
1162 : }
1163 :
1164 : GEN
1165 3476 : FpXQX_invBarrett(GEN S, GEN T, GEN p)
1166 : {
1167 3476 : pari_sp ltop = avma;
1168 3476 : long l = lg(S);
1169 : GEN r;
1170 3476 : if (l<5) return pol_0(varn(S));
1171 3455 : if (l<=FpXQX_INVBARRETT_LIMIT)
1172 : {
1173 2918 : GEN c = gel(S,l-1), ci=gen_1;
1174 2918 : if (!gequal1(c))
1175 : {
1176 1870 : ci = Fq_inv(c, T, p);
1177 1870 : S = FqX_Fq_mul(S, ci, T, p);
1178 1870 : r = FpXQX_invBarrett_basecase(S, T, p);
1179 1870 : r = FqX_Fq_mul(r, ci, T, p);
1180 : } else
1181 1048 : r = FpXQX_invBarrett_basecase(S, T, p);
1182 : }
1183 : else
1184 537 : r = FpXQX_invBarrett_Newton(S, T, p);
1185 3455 : return gerepileupto(ltop, r);
1186 : }
1187 :
1188 : GEN
1189 12685 : FpXQX_get_red(GEN S, GEN T, GEN p)
1190 : {
1191 12685 : if (typ(S)==t_POL && lg(S)>FpXQX_BARRETT_LIMIT)
1192 1174 : retmkvec2(FpXQX_invBarrett(S,T,p),S);
1193 11511 : return S;
1194 : }
1195 :
1196 : /* Compute x mod S where 2 <= degpol(S) <= l+1 <= 2*(degpol(S)-1)
1197 : * and mg is the Barrett inverse of S. */
1198 : static GEN
1199 153180 : FpXQX_divrem_Barrettspec(GEN x, long l, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
1200 : {
1201 : GEN q, r;
1202 153180 : long lt = degpol(S); /*We discard the leading term*/
1203 : long ld, lm, lT, lmg;
1204 153180 : ld = l-lt;
1205 153180 : lm = minss(ld, lgpol(mg));
1206 153180 : lT = ZXX_lgrenormalizespec(S+2,lt);
1207 153180 : lmg = ZXX_lgrenormalizespec(mg+2,lm);
1208 153180 : q = FpXX_recipspec(x+lt,ld,ld); /* q = rec(x) lq<=ld*/
1209 153180 : q = FpXQX_mulspec(q+2,mg+2,T,p,lgpol(q),lmg); /* q = rec(x) * mg lq<=ld+lm*/
1210 153180 : q = FpXX_recipspec(q+2,minss(ld,lgpol(q)),ld); /* q = rec (rec(x) * mg) lq<=ld*/
1211 153180 : if (!pr) return q;
1212 152221 : r = FpXQX_mulspec(q+2,S+2,T,p,lgpol(q),lT); /* r = q*pol lr<=ld+lt*/
1213 152221 : r = FpXX_subspec(x,r+2,p,lt,minss(lt,lgpol(r))); /* r = x - r lr<=lt */
1214 152221 : if (pr == ONLY_REM) return r;
1215 68338 : *pr = r; return q;
1216 : }
1217 :
1218 : static GEN
1219 85254 : FpXQX_divrem_Barrett(GEN x, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
1220 : {
1221 85254 : GEN q = NULL, r = FpXQX_red(x, T, p);
1222 85254 : long l = lgpol(r), lt = degpol(S), lm = 2*lt-1, v = varn(S);
1223 : long i;
1224 85254 : if (l <= lt)
1225 : {
1226 0 : if (pr == ONLY_REM) return r;
1227 0 : if (pr == ONLY_DIVIDES) return signe(r)? NULL: pol_0(v);
1228 0 : if (pr) *pr = r;
1229 0 : return pol_0(v);
1230 : }
1231 85254 : if (lt <= 1)
1232 21 : return FpXQX_divrem_basecase(r,S,T,p,pr);
1233 85233 : if (pr != ONLY_REM && l>lm)
1234 : {
1235 1350 : q = cgetg(l-lt+2, t_POL); q[1] = S[1];
1236 87536 : for (i=0;i<l-lt;i++) gel(q+2,i) = gen_0;
1237 : }
1238 153261 : while (l>lm)
1239 : {
1240 68028 : GEN zr, zq = FpXQX_divrem_Barrettspec(r+2+l-lm,lm,mg,S,T,p,&zr);
1241 68028 : long lz = lgpol(zr);
1242 68028 : if (pr != ONLY_REM)
1243 : {
1244 34904 : long lq = lgpol(zq);
1245 113973 : for(i=0; i<lq; i++) gel(q+2+l-lm,i) = gel(zq,2+i);
1246 : }
1247 284012 : for(i=0; i<lz; i++) gel(r+2+l-lm,i) = gel(zr,2+i);
1248 68028 : l = l-lm+lz;
1249 : }
1250 85233 : if (pr == ONLY_REM)
1251 : {
1252 83883 : if (l > lt)
1253 83883 : r = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,ONLY_REM);
1254 : else
1255 0 : r = FpXQX_renormalize(r, l+2);
1256 83883 : setvarn(r, v); return r;
1257 : }
1258 1350 : if (l > lt)
1259 : {
1260 1269 : GEN zq = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,pr ? &r: NULL);
1261 1269 : if (!q) q = zq;
1262 : else
1263 : {
1264 1269 : long lq = lgpol(zq);
1265 8050 : for(i=0; i<lq; i++) gel(q+2,i) = gel(zq,2+i);
1266 : }
1267 : }
1268 81 : else if (pr)
1269 81 : r = FpX_renormalize(r, l+2);
1270 1350 : setvarn(q, v); q = FpXQX_renormalize(q, lg(q));
1271 1350 : if (pr == ONLY_DIVIDES) return signe(r)? NULL: q;
1272 1350 : if (pr) { setvarn(r, v); *pr = r; }
1273 1350 : return q;
1274 : }
1275 :
1276 : GEN
1277 1070786 : FpXQX_divrem(GEN x, GEN S, GEN T, GEN p, GEN *pr)
1278 : {
1279 : GEN B, y;
1280 : long dy, dx, d;
1281 1070786 : if (pr==ONLY_REM) return FpXQX_rem(x, S, T, p);
1282 1070786 : y = get_FpXQX_red(S, &B);
1283 1070782 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1284 1070746 : if (lgefint(p) == 3)
1285 : {
1286 : GEN a, b, t, z;
1287 1052108 : pari_sp tetpil, av = avma;
1288 1052108 : ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1289 1052029 : z = FlxqX_divrem(a, b, t, pp, pr);
1290 1052061 : if (pr == ONLY_DIVIDES && !z) return gc_NULL(av);
1291 1052061 : tetpil=avma;
1292 1052061 : z = FlxX_to_ZXX(z);
1293 1051916 : if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
1294 1017797 : *pr = FlxX_to_ZXX(*pr);
1295 34119 : else return gerepile(av, tetpil, z);
1296 1017796 : gerepileallsp(av,tetpil,2, pr, &z);
1297 1018068 : return z;
1298 : }
1299 18638 : if (!B && d+3 < FpXQX_DIVREM_BARRETT_LIMIT)
1300 17267 : return FpXQX_divrem_basecase(x,y,T,p,pr);
1301 : else
1302 : {
1303 1371 : pari_sp av=avma;
1304 1371 : GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1305 1371 : GEN q = FpXQX_divrem_Barrett(x,mg,y,T,p,pr);
1306 1371 : if (!q) return gc_NULL(av);
1307 1371 : if (!pr || pr==ONLY_DIVIDES) return gerepilecopy(av, q);
1308 398 : return gc_all(av, 2, &q, pr);
1309 : }
1310 : }
1311 :
1312 : GEN
1313 263338 : FpXQX_rem(GEN x, GEN S, GEN T, GEN p)
1314 : {
1315 263338 : GEN B, y = get_FpXQX_red(S, &B);
1316 263338 : long dy = degpol(y), dx = degpol(x), d = dx-dy;
1317 263338 : if (d < 0) return FpXQX_red(x, T, p);
1318 242723 : if (lgefint(p) == 3)
1319 : {
1320 2736 : pari_sp av = avma;
1321 : GEN a, b, t, z;
1322 2736 : ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1323 2736 : z = FlxqX_rem(a, b, t, pp);
1324 2736 : z = FlxX_to_ZXX(z);
1325 2736 : return gerepileupto(av, z);
1326 : }
1327 239987 : if (!B && d+3 < FpXQX_REM_BARRETT_LIMIT)
1328 156104 : return FpXQX_divrem_basecase(x,y, T, p, ONLY_REM);
1329 : else
1330 : {
1331 83883 : pari_sp av=avma;
1332 83883 : GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1333 83883 : GEN r = FpXQX_divrem_Barrett(x, mg, y, T, p, ONLY_REM);
1334 83883 : return gerepileupto(av, r);
1335 : }
1336 : }
1337 :
1338 : /* x + y*z mod p */
1339 : INLINE GEN
1340 36478 : Fq_addmul(GEN x, GEN y, GEN z, GEN T, GEN p)
1341 : {
1342 : pari_sp av;
1343 36478 : if (!signe(y) || !signe(z)) return Fq_red(x, T, p);
1344 36478 : if (!signe(x)) return Fq_mul(z,y, T, p);
1345 36478 : av = avma;
1346 36478 : return gerepileupto(av, Fq_add(x, Fq_mul(y, z, T, p), T, p));
1347 : }
1348 :
1349 : GEN
1350 71939 : FpXQX_div_by_X_x(GEN a, GEN x, GEN T, GEN p, GEN *r)
1351 : {
1352 71939 : long l = lg(a), i;
1353 : GEN z;
1354 71939 : if (lgefint(p)==3)
1355 : {
1356 53700 : pari_sp av = avma;
1357 : GEN ap, xp, t, z;
1358 53700 : ulong pp = to_FlxqX(a, NULL, T, p, &ap, NULL, &t);
1359 53700 : xp = ZX_to_Flx(to_ZX(x, get_FpX_var(T)), pp);
1360 53700 : z = FlxX_to_ZXX(FlxqX_div_by_X_x(ap, xp, t, pp, r));
1361 53700 : if (!r) return gerepileupto(av, z);
1362 0 : *r = Flx_to_ZX(*r);
1363 0 : return gc_all(av, 2, &z, r);
1364 : }
1365 18239 : if (l <= 3)
1366 : {
1367 0 : if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1368 0 : return pol_0(varn(a));
1369 : }
1370 18239 : l--; z = cgetg(l, t_POL); z[1] = a[1];
1371 18239 : gel(z, l-1) = gel(a,l);
1372 54717 : for (i=l-2; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1373 36478 : gel(z, i) = Fq_addmul(gel(a,i+1), x, gel(z,i+1), T, p);
1374 18239 : if (r) *r = Fq_addmul(gel(a,2), x, gel(z,2), T, p);
1375 18239 : return z;
1376 : }
1377 :
1378 : struct _FpXQXQ {
1379 : GEN T, S;
1380 : GEN p;
1381 : };
1382 :
1383 118643 : static GEN _FpXQX_mul(void *data, GEN a,GEN b)
1384 : {
1385 118643 : struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1386 118643 : return FpXQX_mul(a,b,d->T,d->p);
1387 : }
1388 :
1389 1729 : static GEN _FpXQX_sqr(void *data, GEN a)
1390 : {
1391 1729 : struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1392 1729 : return FpXQX_sqr(a, d->T, d->p);
1393 : }
1394 :
1395 : GEN
1396 63 : FpXQX_powu(GEN x, ulong n, GEN T, GEN p)
1397 : {
1398 : struct _FpXQXQ D;
1399 63 : if (n==0) return pol_1(varn(x));
1400 63 : D.T = T; D.p = p;
1401 63 : return gen_powu(x, n, (void *)&D, _FpXQX_sqr, _FpXQX_mul);
1402 : }
1403 :
1404 : GEN
1405 16717 : FpXQXV_prod(GEN V, GEN T, GEN p)
1406 : {
1407 16717 : if (lgefint(p) == 3)
1408 : {
1409 0 : pari_sp av = avma;
1410 0 : ulong pp = p[2];
1411 0 : GEN Tl = ZXT_to_FlxT(T, pp);
1412 0 : GEN Vl = ZXXV_to_FlxXV(V, pp, get_FpX_var(T));
1413 0 : Tl = FlxqXV_prod(Vl, Tl, pp);
1414 0 : return gerepileupto(av, FlxX_to_ZXX(Tl));
1415 : }
1416 : else
1417 : {
1418 : struct _FpXQXQ d;
1419 16717 : d.T=T; d.p=p;
1420 16717 : return gen_product(V, (void*)&d, &_FpXQX_mul);
1421 : }
1422 : }
1423 :
1424 : static GEN
1425 9954 : _FpXQX_divrem(void * E, GEN x, GEN y, GEN *r)
1426 : {
1427 9954 : struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1428 9954 : return FpXQX_divrem(x, y, d->T, d->p, r);
1429 : }
1430 :
1431 : static GEN
1432 121791 : _FpXQX_add(void * E, GEN x, GEN y)
1433 : {
1434 121791 : struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1435 121791 : return FpXX_add(x, y, d->p);
1436 : }
1437 :
1438 : static GEN
1439 4638 : _FpXQX_sub(void * E, GEN x, GEN y) {
1440 4638 : struct _FpXQXQ *d = (struct _FpXQXQ*) E;
1441 4638 : return FpXX_sub(x,y, d->p);
1442 : }
1443 :
1444 : static struct bb_ring FpXQX_ring = { _FpXQX_add, _FpXQX_mul, _FpXQX_sqr };
1445 :
1446 : GEN
1447 623 : FpXQX_digits(GEN x, GEN B, GEN T, GEN p)
1448 : {
1449 623 : long d = degpol(B), n = (lgpol(x)+d-1)/d;
1450 : struct _FpXQXQ D;
1451 623 : D.T = T; D.p = p;
1452 623 : return gen_digits(x, B, n, (void *)&D, &FpXQX_ring, _FpXQX_divrem);
1453 : }
1454 :
1455 : GEN
1456 189 : FpXQXV_FpXQX_fromdigits(GEN x, GEN B, GEN T, GEN p)
1457 : {
1458 : struct _FpXQXQ D;
1459 189 : D.T = T; D.p = p;
1460 189 : return gen_fromdigits(x,B,(void *)&D, &FpXQX_ring);
1461 : }
1462 :
1463 : /* Q an FpXY (t_POL with FpX coeffs), evaluate at X = x */
1464 : GEN
1465 51552 : FpXY_evalx(GEN Q, GEN x, GEN p)
1466 : {
1467 51552 : long i, lb = lg(Q);
1468 : GEN z;
1469 51552 : z = cgetg(lb, t_POL); z[1] = Q[1];
1470 445045 : for (i=2; i<lb; i++)
1471 : {
1472 393493 : GEN q = gel(Q,i);
1473 393493 : gel(z,i) = typ(q) == t_INT? modii(q,p): FpX_eval(q, x, p);
1474 : }
1475 51552 : return FpX_renormalize(z, lb);
1476 : }
1477 : /* Q an FpXY, evaluate at Y = y */
1478 : GEN
1479 18799 : FpXY_evaly(GEN Q, GEN y, GEN p, long vx)
1480 : {
1481 18799 : pari_sp av = avma;
1482 18799 : long i, lb = lg(Q);
1483 : GEN z;
1484 18799 : if (!signe(Q)) return pol_0(vx);
1485 18771 : if (lb == 3 || !signe(y)) {
1486 84 : z = gel(Q, 2);
1487 84 : return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1488 : }
1489 18687 : z = gel(Q, lb-1);
1490 18687 : if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1491 242061 : for (i=lb-2; i>=2; i--) z = Fq_add(gel(Q,i), FpX_Fp_mul(z, y, p), NULL, p);
1492 18687 : return gerepileupto(av, z);
1493 : }
1494 : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
1495 : GEN
1496 13657 : FpXY_eval(GEN Q, GEN y, GEN x, GEN p)
1497 : {
1498 13657 : pari_sp av = avma;
1499 13657 : return gerepileuptoint(av, FpX_eval(FpXY_evalx(Q, x, p), y, p));
1500 : }
1501 :
1502 : GEN
1503 3852 : FpXY_FpXQV_evalx(GEN P, GEN x, GEN T, GEN p)
1504 : {
1505 3852 : long i, lP = lg(P);
1506 3852 : GEN res = cgetg(lP,t_POL);
1507 3852 : res[1] = P[1];
1508 61377 : for(i=2; i<lP; i++)
1509 115050 : gel(res,i) = typ(gel(P,i))==t_INT? icopy(gel(P,i)):
1510 57525 : FpX_FpXQV_eval(gel(P,i), x, T, p);
1511 3852 : return FlxX_renormalize(res, lP);
1512 : }
1513 :
1514 : GEN
1515 154 : FpXY_FpXQ_evalx(GEN P, GEN x, GEN T, GEN p)
1516 : {
1517 154 : pari_sp av = avma;
1518 154 : long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(P),1);
1519 154 : GEN xp = FpXQ_powers(x, n, T, p);
1520 154 : return gerepileupto(av, FpXY_FpXQV_evalx(P, xp, T, p));
1521 : }
1522 :
1523 : /*******************************************************************/
1524 : /* */
1525 : /* (Fp[X]/T(X))[Y] / S(Y) */
1526 : /* */
1527 : /*******************************************************************/
1528 :
1529 : /*Preliminary implementation to speed up FpX_ffisom*/
1530 : typedef struct {
1531 : GEN S, T, p;
1532 : } FpXYQQ_muldata;
1533 :
1534 : /* reduce x in Fp[X, Y] in the algebra Fp[X,Y]/ (S(X),T(Y)) */
1535 : static GEN
1536 476 : FpXYQQ_redswap(GEN x, GEN S, GEN T, GEN p)
1537 : {
1538 476 : pari_sp ltop=avma;
1539 476 : long n = get_FpX_degree(S);
1540 476 : long m = get_FpX_degree(T);
1541 476 : long v = get_FpX_var(T);
1542 476 : GEN V = RgXY_swap(x,m,v);
1543 476 : V = FpXQX_red(V,S,p);
1544 476 : V = RgXY_swap(V,n,v);
1545 476 : return gerepilecopy(ltop,V);
1546 : }
1547 : static GEN
1548 280 : FpXYQQ_sqr(void *data, GEN x)
1549 : {
1550 280 : FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1551 280 : return FpXYQQ_redswap(FpXQX_sqr(x, D->T, D->p),D->S,D->T,D->p);
1552 :
1553 : }
1554 : static GEN
1555 196 : FpXYQQ_mul(void *data, GEN x, GEN y)
1556 : {
1557 196 : FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1558 196 : return FpXYQQ_redswap(FpXQX_mul(x,y, D->T, D->p),D->S,D->T,D->p);
1559 : }
1560 :
1561 : /* x in Z[X,Y], S in Z[X] over Fq = Z[Y]/(p,T); compute lift(x^n mod (S,T,p)) */
1562 : GEN
1563 182 : FpXYQQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1564 : {
1565 182 : pari_sp av = avma;
1566 : FpXYQQ_muldata D;
1567 : GEN y;
1568 182 : if (lgefint(p) == 3)
1569 : {
1570 0 : ulong pp = to_FlxqX(x, NULL, T, p, &x, NULL, &T);
1571 0 : S = ZX_to_Flx(S, pp);
1572 0 : y = FlxX_to_ZXX( FlxYqq_pow(x, n, S, T, pp) );
1573 0 : y = gerepileupto(av, y);
1574 : }
1575 : else
1576 : {
1577 182 : D.S = S;
1578 182 : D.T = T;
1579 182 : D.p = p;
1580 182 : y = gen_pow(x, n, (void*)&D, &FpXYQQ_sqr, &FpXYQQ_mul);
1581 : }
1582 182 : return y;
1583 : }
1584 :
1585 : GEN
1586 49194 : FpXQXQ_mul(GEN x, GEN y, GEN S, GEN T, GEN p) {
1587 49194 : return FpXQX_rem(FpXQX_mul(x, y, T, p), S, T, p);
1588 : }
1589 :
1590 : GEN
1591 176117 : FpXQXQ_sqr(GEN x, GEN S, GEN T, GEN p) {
1592 176117 : return FpXQX_rem(FpXQX_sqr(x, T, p), S, T, p);
1593 : }
1594 :
1595 : /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
1596 : * return lift(1 / (x mod (p,pol))) */
1597 : GEN
1598 0 : FpXQXQ_invsafe(GEN x, GEN S, GEN T, GEN p)
1599 : {
1600 0 : GEN V, z = FpXQX_extgcd(get_FpXQX_mod(S), x, T, p, NULL, &V);
1601 0 : if (degpol(z)) return NULL;
1602 0 : z = gel(z,2);
1603 0 : z = typ(z)==t_INT ? Fp_invsafe(z,p) : FpXQ_invsafe(z,T,p);
1604 0 : if (!z) return NULL;
1605 0 : return typ(z)==t_INT ? FpXX_Fp_mul(V, z, p): FpXQX_FpXQ_mul(V, z, T, p);
1606 : }
1607 :
1608 : GEN
1609 0 : FpXQXQ_inv(GEN x, GEN S, GEN T,GEN p)
1610 : {
1611 0 : pari_sp av = avma;
1612 0 : GEN U = FpXQXQ_invsafe(x, S, T, p);
1613 0 : if (!U) pari_err_INV("FpXQXQ_inv",x);
1614 0 : return gerepileupto(av, U);
1615 : }
1616 :
1617 : GEN
1618 0 : FpXQXQ_div(GEN x,GEN y,GEN S, GEN T,GEN p)
1619 : {
1620 0 : pari_sp av = avma;
1621 0 : return gerepileupto(av, FpXQXQ_mul(x, FpXQXQ_inv(y,S,T,p),S,T,p));
1622 : }
1623 :
1624 : static GEN
1625 125668 : _FpXQXQ_cmul(void *data, GEN P, long a, GEN x) {
1626 125668 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1627 125668 : GEN y = gel(P,a+2);
1628 251313 : return typ(y)==t_INT ? FpXX_Fp_mul(x,y, d->p):
1629 125645 : FpXX_FpX_mul(x,y,d->p);
1630 : }
1631 : static GEN
1632 22956 : _FpXQXQ_red(void *data, GEN x) {
1633 22956 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1634 22956 : return FpXQX_red(x, d->T, d->p);
1635 : }
1636 : static GEN
1637 45926 : _FpXQXQ_mul(void *data, GEN x, GEN y) {
1638 45926 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1639 45926 : return FpXQXQ_mul(x,y, d->S,d->T, d->p);
1640 : }
1641 : static GEN
1642 176117 : _FpXQXQ_sqr(void *data, GEN x) {
1643 176117 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1644 176117 : return FpXQXQ_sqr(x, d->S,d->T, d->p);
1645 : }
1646 :
1647 : static GEN
1648 18963 : _FpXQXQ_one(void *data) {
1649 18963 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1650 18963 : return pol_1(get_FpXQX_var(d->S));
1651 : }
1652 :
1653 : static GEN
1654 132 : _FpXQXQ_zero(void *data) {
1655 132 : struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1656 132 : return pol_0(get_FpXQX_var(d->S));
1657 : }
1658 :
1659 : static struct bb_algebra FpXQXQ_algebra = { _FpXQXQ_red, _FpXQX_add,
1660 : _FpXQX_sub, _FpXQXQ_mul, _FpXQXQ_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1661 :
1662 : const struct bb_algebra *
1663 352 : get_FpXQXQ_algebra(void **E, GEN S, GEN T, GEN p)
1664 : {
1665 352 : GEN z = new_chunk(sizeof(struct _FpXQXQ));
1666 352 : struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1667 352 : e->T = FpX_get_red(T, p);
1668 352 : e->S = FpXQX_get_red(S, e->T, p);
1669 352 : e->p = p; *E = (void*)e;
1670 352 : return &FpXQXQ_algebra;
1671 : }
1672 :
1673 : static struct bb_algebra FpXQX_algebra = { _FpXQXQ_red, _FpXQX_add,
1674 : _FpXQX_sub, _FpXQX_mul, _FpXQX_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1675 :
1676 : const struct bb_algebra *
1677 0 : get_FpXQX_algebra(void **E, GEN T, GEN p, long v)
1678 : {
1679 0 : GEN z = new_chunk(sizeof(struct _FpXQXQ));
1680 0 : struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1681 0 : e->T = FpX_get_red(T, p);
1682 0 : e->S = pol_x(v);
1683 0 : e->p = p; *E = (void*)e;
1684 0 : return &FpXQX_algebra;
1685 : }
1686 :
1687 : /* x over Fq, return lift(x^n) mod S */
1688 : GEN
1689 1908 : FpXQXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1690 : {
1691 1908 : pari_sp ltop = avma;
1692 : GEN y;
1693 : struct _FpXQXQ D;
1694 1908 : long s = signe(n);
1695 1908 : if (!s) return pol_1(varn(x));
1696 1908 : if (is_pm1(n)) /* +/- 1 */
1697 0 : return (s < 0)? FpXQXQ_inv(x,S,T,p): ZXX_copy(x);
1698 1908 : if (lgefint(p) == 3)
1699 : {
1700 35 : ulong pp = to_FlxqX(x, S, T, p, &x, &S, &T);
1701 35 : GEN z = FlxqXQ_pow(x, n, S, T, pp);
1702 35 : y = FlxX_to_ZXX(z);
1703 35 : return gerepileupto(ltop, y);
1704 : }
1705 : else
1706 : {
1707 1873 : T = FpX_get_red(T, p);
1708 1873 : S = FpXQX_get_red(S, T, p);
1709 1873 : D.S = S; D.T = T; D.p = p;
1710 1873 : if (s < 0) x = FpXQXQ_inv(x,S,T,p);
1711 1873 : y = gen_pow_i(x, n, (void*)&D,&_FpXQXQ_sqr,&_FpXQXQ_mul);
1712 1873 : return gerepilecopy(ltop, y);
1713 : }
1714 : }
1715 :
1716 : /* generates the list of powers of x of degree 0,1,2,...,l*/
1717 : GEN
1718 1965 : FpXQXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
1719 : {
1720 : struct _FpXQXQ D;
1721 1965 : int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1722 1965 : T = FpX_get_red(T, p);
1723 1965 : S = FpXQX_get_red(S, T, p);
1724 1965 : D.S = S; D.T = T; D.p = p;
1725 1965 : return gen_powers(x, l, use_sqr, (void*)&D, &_FpXQXQ_sqr, &_FpXQXQ_mul,&_FpXQXQ_one);
1726 : }
1727 :
1728 : /* Let v a linear form, return the linear form z->v(tau*z)
1729 : that is, v*(M_tau) */
1730 :
1731 : INLINE GEN
1732 248 : FpXQX_recipspec(GEN x, long l, long n)
1733 : {
1734 248 : return RgX_recipspec_shallow(x, l, n);
1735 : }
1736 :
1737 : static GEN
1738 88 : FpXQXQ_transmul_init(GEN tau, GEN S, GEN T, GEN p)
1739 : {
1740 : GEN bht;
1741 88 : GEN h, Sp = get_FpXQX_red(S, &h);
1742 88 : long n = degpol(Sp), vT = varn(Sp);
1743 88 : GEN ft = FpXQX_recipspec(Sp+2, n+1, n+1);
1744 88 : GEN bt = FpXQX_recipspec(tau+2, lgpol(tau), n);
1745 88 : setvarn(ft, vT); setvarn(bt, vT);
1746 88 : if (h)
1747 16 : bht = FpXQXn_mul(bt, h, n-1, T, p);
1748 : else
1749 : {
1750 72 : GEN bh = FpXQX_div(FpXX_shift(tau, n-1), S, T, p);
1751 72 : bht = FpXQX_recipspec(bh+2, lgpol(bh), n-1);
1752 72 : setvarn(bht, vT);
1753 : }
1754 88 : return mkvec3(bt, bht, ft);
1755 : }
1756 :
1757 : static GEN
1758 200 : FpXQXQ_transmul(GEN tau, GEN a, long n, GEN T, GEN p)
1759 : {
1760 200 : pari_sp ltop = avma;
1761 : GEN t1, t2, t3, vec;
1762 200 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
1763 200 : if (signe(a)==0) return pol_0(varn(a));
1764 200 : t2 = FpXX_shift(FpXQX_mul(bt, a, T, p),1-n);
1765 200 : if (signe(bht)==0) return gerepilecopy(ltop, t2);
1766 114 : t1 = FpXX_shift(FpXQX_mul(ft, a, T, p),-n);
1767 114 : t3 = FpXQXn_mul(t1, bht, n-1, T, p);
1768 114 : vec = FpXX_sub(t2, FpXX_shift(t3, 1), p);
1769 114 : return gerepileupto(ltop, vec);
1770 : }
1771 :
1772 : static GEN
1773 44 : polxn_FpXX(long n, long v, long vT)
1774 : {
1775 44 : long i, a = n+2;
1776 44 : GEN p = cgetg(a+1, t_POL);
1777 44 : p[1] = evalsigne(1)|evalvarn(v);
1778 440 : for (i = 2; i < a; i++) gel(p,i) = pol_0(vT);
1779 44 : gel(p,a) = pol_1(vT); return p;
1780 : }
1781 :
1782 : GEN
1783 44 : FpXQXQ_minpoly(GEN x, GEN S, GEN T, GEN p)
1784 : {
1785 44 : pari_sp ltop = avma;
1786 : long vS, vT, n;
1787 : GEN v_x, g, tau;
1788 44 : vS = get_FpXQX_var(S);
1789 44 : vT = get_FpX_var(T);
1790 44 : n = get_FpXQX_degree(S);
1791 44 : g = pol_1(vS);
1792 44 : tau = pol_1(vS);
1793 44 : S = FpXQX_get_red(S, T, p);
1794 44 : v_x = FpXQXQ_powers(x, usqrt(2*n), S, T, p);
1795 88 : while(signe(tau) != 0)
1796 : {
1797 : long i, j, m, k1;
1798 : GEN M, v, tr;
1799 : GEN g_prime, c;
1800 44 : if (degpol(g) == n) { tau = pol_1(vS); g = pol_1(vS); }
1801 44 : v = random_FpXQX(n, vS, T, p);
1802 44 : tr = FpXQXQ_transmul_init(tau, S, T, p);
1803 44 : v = FpXQXQ_transmul(tr, v, n, T, p);
1804 44 : m = 2*(n-degpol(g));
1805 44 : k1 = usqrt(m);
1806 44 : tr = FpXQXQ_transmul_init(gel(v_x,k1+1), S, T, p);
1807 44 : c = cgetg(m+2,t_POL);
1808 44 : c[1] = evalsigne(1)|evalvarn(vS);
1809 200 : for (i=0; i<m; i+=k1)
1810 : {
1811 156 : long mj = minss(m-i, k1);
1812 552 : for (j=0; j<mj; j++)
1813 396 : gel(c,m+1-(i+j)) = FpXQX_dotproduct(v, gel(v_x,j+1), T, p);
1814 156 : v = FpXQXQ_transmul(tr, v, n, T, p);
1815 : }
1816 44 : c = FpXX_renormalize(c, m+2);
1817 : /* now c contains <v,x^i> , i = 0..m-1 */
1818 44 : M = FpXQX_halfgcd(polxn_FpXX(m, vS, vT), c, T, p);
1819 44 : g_prime = gmael(M, 2, 2);
1820 44 : if (degpol(g_prime) < 1) continue;
1821 44 : g = FpXQX_mul(g, g_prime, T, p);
1822 44 : tau = FpXQXQ_mul(tau, FpXQX_FpXQXQV_eval(g_prime, v_x, S, T, p), S, T, p);
1823 : }
1824 44 : g = FpXQX_normalize(g,T, p);
1825 44 : return gerepilecopy(ltop,g);
1826 : }
1827 :
1828 : GEN
1829 0 : FpXQXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
1830 : {
1831 0 : return RgXV_to_RgM(FpXQXQ_powers(y,m-1,S,T,p),n);
1832 : }
1833 :
1834 : GEN
1835 4191 : FpXQX_FpXQXQV_eval(GEN P, GEN V, GEN S, GEN T, GEN p)
1836 : {
1837 : struct _FpXQXQ D;
1838 4191 : T = FpX_get_red(T, p);
1839 4191 : S = FpXQX_get_red(S, T, p);
1840 4191 : D.S=S; D.T=T; D.p=p;
1841 4191 : return gen_bkeval_powers(P, degpol(P), V, (void*)&D, &FpXQXQ_algebra,
1842 : _FpXQXQ_cmul);
1843 : }
1844 :
1845 : GEN
1846 855 : FpXQX_FpXQXQ_eval(GEN Q, GEN x, GEN S, GEN T, GEN p)
1847 : {
1848 : struct _FpXQXQ D;
1849 855 : int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1850 855 : T = FpX_get_red(T, p);
1851 855 : S = FpXQX_get_red(S, T, p);
1852 855 : D.S=S; D.T=T; D.p=p;
1853 855 : return gen_bkeval(Q, degpol(Q), x, use_sqr, (void*)&D, &FpXQXQ_algebra,
1854 : _FpXQXQ_cmul);
1855 : }
1856 :
1857 : static GEN
1858 507 : FpXQXQ_autpow_sqr(void * E, GEN x)
1859 : {
1860 507 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1861 507 : GEN S = D->S, T = D->T, p = D->p;
1862 507 : GEN phi = gel(x,1), S1 = gel(x,2);
1863 507 : long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(S1)+1,1);
1864 507 : GEN V = FpXQ_powers(phi, n, T, p);
1865 507 : GEN phi2 = FpX_FpXQV_eval(phi, V, T, p);
1866 507 : GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1867 507 : GEN S2 = FpXQX_FpXQXQ_eval(Sphi, S1, S, T, p);
1868 507 : return mkvec2(phi2, S2);
1869 : }
1870 :
1871 : static GEN
1872 325 : FpXQXQ_autpow_mul(void * E, GEN x, GEN y)
1873 : {
1874 325 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1875 325 : GEN S = D->S, T = D->T, p = D->p;
1876 325 : GEN phi1 = gel(x,1), S1 = gel(x,2);
1877 325 : GEN phi2 = gel(y,1), S2 = gel(y,2);
1878 325 : long n = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+1, 1);
1879 325 : GEN V = FpXQ_powers(phi2, n, T, p);
1880 325 : GEN phi3 = FpX_FpXQV_eval(phi1, V, T, p);
1881 325 : GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1882 325 : GEN S3 = FpXQX_FpXQXQ_eval(Sphi, S2, S, T, p);
1883 325 : return mkvec2(phi3, S3);
1884 : }
1885 :
1886 : GEN
1887 493 : FpXQXQ_autpow(GEN aut, long n, GEN S, GEN T, GEN p)
1888 : {
1889 493 : pari_sp av = avma;
1890 : struct _FpXQXQ D;
1891 493 : T = FpX_get_red(T, p);
1892 493 : S = FpXQX_get_red(S, T, p);
1893 493 : D.S=S; D.T=T; D.p=p;
1894 493 : aut = gen_powu_i(aut,n,&D,FpXQXQ_autpow_sqr,FpXQXQ_autpow_mul);
1895 493 : return gerepilecopy(av, aut);
1896 : }
1897 :
1898 : static GEN
1899 1 : FpXQXQ_auttrace_mul(void *E, GEN x, GEN y)
1900 : {
1901 1 : struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1902 1 : GEN S = D->S, T = D->T;
1903 1 : GEN p = D->p;
1904 1 : GEN S1 = gel(x,1), a1 = gel(x,2);
1905 1 : GEN S2 = gel(y,1), a2 = gel(y,2);
1906 1 : long n = brent_kung_optpow(maxss(degpol(S1),degpol(a1)),2,1);
1907 1 : GEN V = FpXQXQ_powers(S2, n, S, T, p);
1908 1 : GEN S3 = FpXQX_FpXQXQV_eval(S1, V, S, T, p);
1909 1 : GEN aS = FpXQX_FpXQXQV_eval(a1, V, S, T, p);
1910 1 : GEN a3 = FpXX_add(aS, a2, p);
1911 1 : return mkvec2(S3, a3);
1912 : }
1913 :
1914 : static GEN
1915 1 : FpXQXQ_auttrace_sqr(void *E, GEN x)
1916 1 : { return FpXQXQ_auttrace_mul(E, x, x); }
1917 :
1918 : GEN
1919 8 : FpXQXQ_auttrace(GEN aut, long n, GEN S, GEN T, GEN p)
1920 : {
1921 8 : pari_sp av = avma;
1922 : struct _FpXQXQ D;
1923 8 : T = FpX_get_red(T, p);
1924 8 : S = FpXQX_get_red(S, T, p);
1925 8 : D.S=S; D.T=T; D.p=p;
1926 8 : aut = gen_powu_i(aut,n,&D,FpXQXQ_auttrace_sqr,FpXQXQ_auttrace_mul);
1927 8 : return gerepilecopy(av, aut);
1928 : }
1929 :
1930 : static GEN
1931 1433 : FpXQXQ_autsum_mul(void *E, GEN x, GEN y)
1932 : {
1933 1433 : struct _FpXQXQ *D = (struct _FpXQXQ *) E;
1934 1433 : GEN S = D->S, T = D->T, p = D->p;
1935 1433 : GEN phi1 = gel(x,1), S1 = gel(x,2), a1 = gel(x,3);
1936 1433 : GEN phi2 = gel(y,1), S2 = gel(y,2), a2 = gel(y,3);
1937 1433 : long n2 = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+lgpol(a1)+1, 1);
1938 1433 : GEN V2 = FpXQ_powers(phi2, n2, T, p);
1939 1433 : GEN phi3 = FpX_FpXQV_eval(phi1, V2, T, p);
1940 1433 : GEN Sphi = FpXY_FpXQV_evalx(S1, V2, T, p);
1941 1433 : GEN aphi = FpXY_FpXQV_evalx(a1, V2, T, p);
1942 1433 : long n = brent_kung_optpow(maxss(degpol(Sphi),degpol(aphi)),2,1);
1943 1433 : GEN V = FpXQXQ_powers(S2, n, S, T, p);
1944 1433 : GEN S3 = FpXQX_FpXQXQV_eval(Sphi, V, S, T, p);
1945 1433 : GEN aS = FpXQX_FpXQXQV_eval(aphi, V, S, T, p);
1946 1433 : GEN a3 = FpXQXQ_mul(aS, a2, S, T, p);
1947 1433 : return mkvec3(phi3, S3, a3);
1948 : }
1949 :
1950 : static GEN
1951 1275 : FpXQXQ_autsum_sqr(void * T, GEN x)
1952 1275 : { return FpXQXQ_autsum_mul(T,x,x); }
1953 :
1954 : GEN
1955 1233 : FpXQXQ_autsum(GEN aut, long n, GEN S, GEN T, GEN p)
1956 : {
1957 1233 : pari_sp av = avma;
1958 : struct _FpXQXQ D;
1959 1233 : T = FpX_get_red(T, p);
1960 1233 : S = FpXQX_get_red(S, T, p);
1961 1233 : D.S=S; D.T=T; D.p=p;
1962 1233 : aut = gen_powu_i(aut,n,&D,FpXQXQ_autsum_sqr,FpXQXQ_autsum_mul);
1963 1233 : return gerepilecopy(av, aut);
1964 : }
1965 :
1966 : GEN
1967 44573 : FpXQXn_mul(GEN x, GEN y, long n, GEN T, GEN p)
1968 : {
1969 44573 : pari_sp av = avma;
1970 : GEN z, kx, ky;
1971 : long d;
1972 44573 : if (ZXX_is_ZX(y) && ZXX_is_ZX(x))
1973 6937 : return FpXn_mul(x,y,n,p);
1974 37636 : d = get_FpX_degree(T);
1975 37636 : kx = RgXX_to_Kronecker(x, d);
1976 37636 : ky = RgXX_to_Kronecker(y, d);
1977 37636 : z = Kronecker_to_FpXQX(ZXn_mul(ky,kx,(2*d-1)*n), T, p);
1978 37636 : return gerepileupto(av, z);
1979 : }
1980 :
1981 : GEN
1982 0 : FpXQXn_sqr(GEN x, long n, GEN T, GEN p)
1983 : {
1984 0 : pari_sp av = avma;
1985 : GEN z, kx;
1986 : long d;
1987 0 : if (ZXX_is_ZX(x)) return ZXn_sqr(x, n);
1988 0 : d = get_FpX_degree(T);
1989 0 : kx= RgXX_to_Kronecker(x, d);
1990 0 : z = Kronecker_to_FpXQX(ZXn_sqr(kx, (2*d-1)*n), T, p);
1991 0 : return gerepileupto(av, z);
1992 : }
1993 :
1994 : /* (f*g) \/ x^n */
1995 : static GEN
1996 7399 : FpXQX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
1997 : {
1998 7399 : return FpXX_shift(FpXQX_mul(f,g,T, p),-n);
1999 : }
2000 :
2001 : static GEN
2002 4697 : FpXQXn_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
2003 : {
2004 4697 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2005 4697 : return FpXX_add(FpXQX_mulhigh_i(fl, g, n2, T, p), FpXQXn_mul(fh, g, n - n2, T, p), p);
2006 : }
2007 :
2008 : /* Compute intformal(x^n*S)/x^(n+1) */
2009 : static GEN
2010 763 : FpXX_integXn(GEN x, long n, GEN p)
2011 : {
2012 763 : long i, lx = lg(x);
2013 : GEN y;
2014 763 : if (lx == 2) return ZXX_copy(x);
2015 763 : y = cgetg(lx, t_POL); y[1] = x[1];
2016 4317 : for (i=2; i<lx; i++)
2017 : {
2018 3554 : ulong j = n+i-1;
2019 3554 : GEN xi = gel(x,i);
2020 3554 : if (!signe(xi))
2021 0 : gel(y,i) = gen_0;
2022 : else
2023 3554 : gel(y,i) = typ(xi)==t_INT ? Fp_divu(xi, j, p)
2024 3554 : : FpX_divu(xi, j, p);
2025 : }
2026 763 : return ZXX_renormalize(y, lx);;
2027 : }
2028 :
2029 : /* Compute intformal(x^n*S)/x^(n+1) */
2030 : static GEN
2031 2702 : ZlXX_integXn(GEN x, long n, GEN p, ulong pp)
2032 : {
2033 2702 : long i, lx = lg(x);
2034 : GEN y;
2035 2702 : if (lx == 2) return ZXX_copy(x);
2036 2576 : if (!pp) return FpXX_integXn(x, n, p);
2037 1813 : y = cgetg(lx, t_POL); y[1] = x[1];
2038 6883 : for (i=2; i<lx; i++)
2039 : {
2040 5070 : GEN xi = gel(x,i);
2041 5070 : if (!signe(xi))
2042 14 : gel(y,i) = gen_0;
2043 : else
2044 : {
2045 : ulong j;
2046 5056 : long v = u_lvalrem(n+i-1, pp, &j);
2047 5056 : if (typ(xi)==t_INT)
2048 : {
2049 0 : if (v==0)
2050 0 : gel(y,i) = Fp_divu(xi, j, p);
2051 : else
2052 0 : gel(y,i) = Fp_divu(diviuexact(xi, upowuu(pp, v)), j, p);
2053 : } else
2054 : {
2055 5056 : if (v==0)
2056 5056 : gel(y,i) = FpX_divu(xi, j, p);
2057 : else
2058 0 : gel(y,i) = FpX_divu(ZX_divuexact(xi, upowuu(pp, v)), j, p);
2059 : }
2060 : }
2061 : }
2062 1813 : return ZXX_renormalize(y, lx);;
2063 : }
2064 :
2065 : GEN
2066 714 : ZlXQXn_expint(GEN h, long e, GEN T, GEN p, ulong pp)
2067 : {
2068 714 : pari_sp av = avma, av2;
2069 714 : long v = varn(h), n=1;
2070 714 : GEN f = pol_1(v), g = pol_1(v);
2071 714 : ulong mask = quadratic_prec_mask(e);
2072 714 : av2 = avma;
2073 2702 : for (;mask>1;)
2074 : {
2075 : GEN u, w;
2076 2702 : long n2 = n;
2077 2702 : n<<=1; if (mask & 1) n--;
2078 2702 : mask >>= 1;
2079 2702 : u = FpXQXn_mul(g, FpXQX_mulhigh_i(f, FpXXn_red(h, n2-1), n2-1, T, p), n-n2, T, p);
2080 2702 : u = FpXX_add(u, FpXX_shift(FpXXn_red(h, n-1), 1-n2), p);
2081 2702 : w = FpXQXn_mul(f, ZlXX_integXn(u, n2-1, p, pp), n-n2, T, p);
2082 2702 : f = FpXX_add(f, FpXX_shift(w, n2), p);
2083 2702 : if (mask<=1) break;
2084 1988 : u = FpXQXn_mul(g, FpXQXn_mulhigh(f, g, n2, n, T, p), n-n2, T, p);
2085 1988 : g = FpXX_sub(g, FpXX_shift(u, n2), p);
2086 1988 : if (gc_needed(av2,2))
2087 : {
2088 0 : if (DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_exp, e = %ld", n);
2089 0 : gerepileall(av2, 2, &f, &g);
2090 : }
2091 : }
2092 714 : return gerepileupto(av, f);
2093 : }
2094 :
2095 : GEN
2096 178 : FpXQXn_expint(GEN h, long e, GEN T, GEN p)
2097 178 : { return ZlXQXn_expint(h, e, T, p, 0); }
2098 :
2099 : GEN
2100 0 : FpXQXn_exp(GEN h, long e, GEN T, GEN p)
2101 : {
2102 0 : if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
2103 0 : pari_err_DOMAIN("FpXQXn_exp","valuation", "<", gen_1, h);
2104 0 : return FpXQXn_expint(FpXX_deriv(h, p), e, T, p);
2105 : }
2106 :
2107 : GEN
2108 721 : FpXQXn_div(GEN g, GEN f, long e, GEN T, GEN p)
2109 : {
2110 721 : pari_sp av = avma, av2;
2111 : ulong mask;
2112 : GEN W, a;
2113 721 : long v = varn(f), n = 1;
2114 :
2115 721 : if (!signe(f)) pari_err_INV("FpXXn_inv",f);
2116 721 : a = Fq_inv(gel(f,2), T, p);
2117 721 : if (e == 1 && !g) return scalarpol(a, v);
2118 721 : else if (e == 2 && !g)
2119 : {
2120 : GEN b;
2121 0 : if (degpol(f) <= 0) return scalarpol(a, v);
2122 0 : b = Fq_neg(gel(f,3),T,p);
2123 0 : if (signe(b)==0) return scalarpol(a, v);
2124 0 : if (!is_pm1(a)) b = Fq_mul(b, Fq_sqr(a, T, p), T, p);
2125 0 : W = deg1pol_shallow(b, a, v);
2126 0 : return gerepilecopy(av, W);
2127 : }
2128 721 : W = scalarpol_shallow(Fq_inv(gel(f,2), T, p),v);
2129 721 : mask = quadratic_prec_mask(e);
2130 721 : av2 = avma;
2131 3430 : for (;mask>1;)
2132 : {
2133 : GEN u, fr;
2134 2709 : long n2 = n;
2135 2709 : n<<=1; if (mask & 1) n--;
2136 2709 : mask >>= 1;
2137 2709 : fr = FpXXn_red(f, n);
2138 2709 : if (mask>1 || !g)
2139 : {
2140 2702 : u = FpXQXn_mul(W, FpXQXn_mulhigh(fr, W, n2, n, T, p), n-n2, T, p);
2141 2702 : W = FpXX_sub(W, FpXX_shift(u, n2), p);
2142 : }
2143 : else
2144 : {
2145 7 : GEN y = FpXQXn_mul(g, W, n, T, p), yt = FpXXn_red(y, n-n2);
2146 7 : u = FpXQXn_mul(yt, FpXQXn_mulhigh(fr, W, n2, n, T, p), n-n2, T, p);
2147 7 : W = FpXX_sub(y, FpXX_shift(u, n2), p);
2148 : }
2149 2709 : if (gc_needed(av2,2))
2150 : {
2151 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_inv, e = %ld", n);
2152 0 : W = gerepileupto(av2, W);
2153 : }
2154 : }
2155 721 : return gerepileupto(av, W);
2156 : }
2157 :
2158 : GEN
2159 714 : FpXQXn_inv(GEN f, long e, GEN T, GEN p)
2160 714 : { return FpXQXn_div(NULL, f, e, T, p); }
|
