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