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