1 /* Copyright (C) 2007 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 Fp */
19
20 static GEN
get_FpX_red(GEN T,GEN * B)21 get_FpX_red(GEN T, GEN *B)
22 {
23 if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 *B = gel(T,1); return gel(T,2);
25 }
26
27 /***********************************************************************/
28 /** **/
29 /** FpX **/
30 /** **/
31 /***********************************************************************/
32
33 /* FpX are polynomials over Z/pZ represented as t_POL with
34 * t_INT coefficients.
35 * 1) Coefficients should belong to {0,...,p-1}, though nonreduced
36 * coefficients should work but be slower.
37 *
38 * 2) p is not assumed to be prime, but it is assumed that impossible divisions
39 * will not happen.
40 * 3) Theses functions let some garbage on the stack, but are gerepileupto
41 * compatible.
42 */
43
44 static ulong
to_Flx(GEN * P,GEN * Q,GEN p)45 to_Flx(GEN *P, GEN *Q, GEN p)
46 {
47 ulong pp = uel(p,2);
48 *P = ZX_to_Flx(*P, pp);
49 if(Q) *Q = ZX_to_Flx(*Q, pp);
50 return pp;
51 }
52
53 static ulong
to_Flxq(GEN * P,GEN * T,GEN p)54 to_Flxq(GEN *P, GEN *T, GEN p)
55 {
56 ulong pp = uel(p,2);
57 if (P) *P = ZX_to_Flx(*P, pp);
58 *T = ZXT_to_FlxT(*T, pp); return pp;
59 }
60
61 GEN
Z_to_FpX(GEN a,GEN p,long v)62 Z_to_FpX(GEN a, GEN p, long v)
63 {
64 pari_sp av = avma;
65 GEN z = cgetg(3, t_POL);
66 GEN x = modii(a, p);
67 if (!signe(x)) { set_avma(av); return pol_0(v); }
68 z[1] = evalsigne(1) | evalvarn(v);
69 gel(z,2) = x; return z;
70 }
71
72 /* z in Z[X], return lift(z * Mod(1,p)), normalized*/
73 GEN
FpX_red(GEN z,GEN p)74 FpX_red(GEN z, GEN p)
75 {
76 long i, l = lg(z);
77 GEN x = cgetg(l, t_POL);
78 for (i=2; i<l; i++) gel(x,i) = modii(gel(z,i),p);
79 x[1] = z[1]; return FpX_renormalize(x,l);
80 }
81
82 GEN
FpXV_red(GEN x,GEN p)83 FpXV_red(GEN x, GEN p)
84 { pari_APPLY_type(t_VEC, FpX_red(gel(x,i), p)) }
85
86 GEN
FpXT_red(GEN x,GEN p)87 FpXT_red(GEN x, GEN p)
88 {
89 if (typ(x) == t_POL)
90 return FpX_red(x, p);
91 else
92 pari_APPLY_type(t_VEC, FpXT_red(gel(x,i), p))
93 }
94
95 GEN
FpX_normalize(GEN z,GEN p)96 FpX_normalize(GEN z, GEN p)
97 {
98 GEN p1 = leading_coeff(z);
99 if (lg(z) == 2 || equali1(p1)) return z;
100 return FpX_Fp_mul_to_monic(z, Fp_inv(p1,p), p);
101 }
102
103 GEN
FpX_center(GEN T,GEN p,GEN pov2)104 FpX_center(GEN T, GEN p, GEN pov2)
105 {
106 long i, l = lg(T);
107 GEN P = cgetg(l,t_POL);
108 for(i=2; i<l; i++) gel(P,i) = Fp_center(gel(T,i), p, pov2);
109 P[1] = T[1]; return P;
110 }
111 GEN
FpX_center_i(GEN T,GEN p,GEN pov2)112 FpX_center_i(GEN T, GEN p, GEN pov2)
113 {
114 long i, l = lg(T);
115 GEN P = cgetg(l,t_POL);
116 for(i=2; i<l; i++) gel(P,i) = Fp_center_i(gel(T,i), p, pov2);
117 P[1] = T[1]; return P;
118 }
119
120 GEN
FpX_add(GEN x,GEN y,GEN p)121 FpX_add(GEN x,GEN y,GEN p)
122 {
123 long lx = lg(x), ly = lg(y), i;
124 GEN z;
125 if (lx < ly) swapspec(x,y, lx,ly);
126 z = cgetg(lx,t_POL); z[1] = x[1];
127 for (i=2; i<ly; i++) gel(z,i) = Fp_add(gel(x,i),gel(y,i), p);
128 for ( ; i<lx; i++) gel(z,i) = modii(gel(x,i), p);
129 z = ZX_renormalize(z, lx);
130 if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
131 return z;
132 }
133
134 static GEN
Fp_red_FpX(GEN x,GEN p,long v)135 Fp_red_FpX(GEN x, GEN p, long v)
136 {
137 GEN z;
138 if (!signe(x)) return pol_0(v);
139 z = cgetg(3, t_POL);
140 gel(z,2) = Fp_red(x,p);
141 z[1] = evalvarn(v);
142 return FpX_renormalize(z, 3);
143 }
144
145 static GEN
Fp_neg_FpX(GEN x,GEN p,long v)146 Fp_neg_FpX(GEN x, GEN p, long v)
147 {
148 GEN z;
149 if (!signe(x)) return pol_0(v);
150 z = cgetg(3, t_POL);
151 gel(z,2) = Fp_neg(x,p);
152 z[1] = evalvarn(v);
153 return FpX_renormalize(z, 3);
154 }
155
156 GEN
FpX_Fp_add(GEN y,GEN x,GEN p)157 FpX_Fp_add(GEN y,GEN x,GEN p)
158 {
159 long i, lz = lg(y);
160 GEN z;
161 if (lz == 2) return Fp_red_FpX(x,p,varn(y));
162 z = cgetg(lz,t_POL); z[1] = y[1];
163 gel(z,2) = Fp_add(gel(y,2),x, p);
164 if (lz == 3) z = FpX_renormalize(z,lz);
165 else
166 for(i=3;i<lz;i++) gel(z,i) = icopy(gel(y,i));
167 return z;
168 }
169 GEN
FpX_Fp_add_shallow(GEN y,GEN x,GEN p)170 FpX_Fp_add_shallow(GEN y,GEN x,GEN p)
171 {
172 long i, lz = lg(y);
173 GEN z;
174 if (lz == 2) return scalar_ZX_shallow(x,varn(y));
175 z = cgetg(lz,t_POL); z[1] = y[1];
176 gel(z,2) = Fp_add(gel(y,2),x, p);
177 if (lz == 3) z = FpX_renormalize(z,lz);
178 else
179 for(i=3;i<lz;i++) gel(z,i) = gel(y,i);
180 return z;
181 }
182 GEN
FpX_Fp_sub(GEN y,GEN x,GEN p)183 FpX_Fp_sub(GEN y,GEN x,GEN p)
184 {
185 long i, lz = lg(y);
186 GEN z;
187 if (lz == 2) return Fp_neg_FpX(x,p,varn(y));
188 z = cgetg(lz,t_POL); z[1] = y[1];
189 gel(z,2) = Fp_sub(gel(y,2),x, p);
190 if (lz == 3) z = FpX_renormalize(z,lz);
191 else
192 for(i=3;i<lz;i++) gel(z,i) = icopy(gel(y,i));
193 return z;
194 }
195 GEN
FpX_Fp_sub_shallow(GEN y,GEN x,GEN p)196 FpX_Fp_sub_shallow(GEN y,GEN x,GEN p)
197 {
198 long i, lz = lg(y);
199 GEN z;
200 if (lz == 2) return Fp_neg_FpX(x,p,varn(y));
201 z = cgetg(lz,t_POL); z[1] = y[1];
202 gel(z,2) = Fp_sub(gel(y,2),x, p);
203 if (lz == 3) z = FpX_renormalize(z,lz);
204 else
205 for(i=3;i<lz;i++) gel(z,i) = gel(y,i);
206 return z;
207 }
208
209 GEN
FpX_neg(GEN x,GEN p)210 FpX_neg(GEN x,GEN p)
211 {
212 long i, lx = lg(x);
213 GEN y = cgetg(lx,t_POL);
214 y[1] = x[1];
215 for(i=2; i<lx; i++) gel(y,i) = Fp_neg(gel(x,i), p);
216 return ZX_renormalize(y, lx);
217 }
218
219 static GEN
FpX_subspec(GEN x,GEN y,GEN p,long nx,long ny)220 FpX_subspec(GEN x,GEN y,GEN p, long nx, long ny)
221 {
222 long i, lz;
223 GEN z;
224 if (nx >= ny)
225 {
226 lz = nx+2;
227 z = cgetg(lz,t_POL); z[1] = 0; z += 2;
228 for (i=0; i<ny; i++) gel(z,i) = Fp_sub(gel(x,i),gel(y,i), p);
229 for ( ; i<nx; i++) gel(z,i) = modii(gel(x,i), p);
230 }
231 else
232 {
233 lz = ny+2;
234 z = cgetg(lz,t_POL); z[1] = 0; z += 2;
235 for (i=0; i<nx; i++) gel(z,i) = Fp_sub(gel(x,i),gel(y,i), p);
236 for ( ; i<ny; i++) gel(z,i) = Fp_neg(gel(y,i), p);
237 }
238 z = FpX_renormalize(z-2, lz);
239 if (!lgpol(z)) { set_avma((pari_sp)(z + lz)); return pol_0(0); }
240 return z;
241 }
242
243 GEN
FpX_sub(GEN x,GEN y,GEN p)244 FpX_sub(GEN x,GEN y,GEN p)
245 {
246 GEN z = FpX_subspec(x+2,y+2,p,lgpol(x),lgpol(y));
247 setvarn(z, varn(x));
248 return z;
249 }
250
251 GEN
Fp_FpX_sub(GEN x,GEN y,GEN p)252 Fp_FpX_sub(GEN x, GEN y, GEN p)
253 {
254 long ly = lg(y), i;
255 GEN z;
256 if (ly <= 3) {
257 z = cgetg(3, t_POL);
258 x = (ly == 3)? Fp_sub(x, gel(y,2), p): modii(x, p);
259 if (!signe(x)) { set_avma((pari_sp)(z + 3)); return pol_0(varn(y)); }
260 z[1] = evalsigne(1)|y[1]; gel(z,2) = x; return z;
261 }
262 z = cgetg(ly,t_POL);
263 gel(z,2) = Fp_sub(x, gel(y,2), p);
264 for (i = 3; i < ly; i++) gel(z,i) = Fp_neg(gel(y,i), p);
265 z = ZX_renormalize(z, ly);
266 if (!lgpol(z)) { set_avma((pari_sp)(z + ly)); return pol_0(varn(x)); }
267 z[1] = y[1]; return z;
268 }
269
270 GEN
FpX_convol(GEN x,GEN y,GEN p)271 FpX_convol(GEN x, GEN y, GEN p)
272 {
273 long lx = lg(x), ly = lg(y), i;
274 GEN z;
275 if (lx < ly) swapspec(x,y, lx,ly);
276 z = cgetg(lx,t_POL); z[1] = x[1];
277 for (i=2; i<ly; i++) gel(z,i) = Fp_mul(gel(x,i),gel(y,i), p);
278 for ( ; i<lx; i++) gel(z,i) = modii(gel(x,i), p);
279 z = ZX_renormalize(z, lx);
280 if (!lgpol(z)) { set_avma((pari_sp)(z + lx)); return pol_0(varn(x)); }
281 return z;
282 }
283
284 GEN
FpX_mul(GEN x,GEN y,GEN p)285 FpX_mul(GEN x,GEN y,GEN p)
286 {
287 if (lgefint(p) == 3)
288 {
289 ulong pp = to_Flx(&x, &y, p);
290 return Flx_to_ZX(Flx_mul(x, y, pp));
291 }
292 return FpX_red(ZX_mul(x, y), p);
293 }
294
295 GEN
FpX_mulspec(GEN a,GEN b,GEN p,long na,long nb)296 FpX_mulspec(GEN a, GEN b, GEN p, long na, long nb)
297 { return FpX_red(ZX_mulspec(a, b, na, nb), p); }
298
299 GEN
FpX_sqr(GEN x,GEN p)300 FpX_sqr(GEN x,GEN p)
301 {
302 if (lgefint(p) == 3)
303 {
304 ulong pp = to_Flx(&x, NULL, p);
305 return Flx_to_ZX(Flx_sqr(x, pp));
306 }
307 return FpX_red(ZX_sqr(x), p);
308 }
309
310 GEN
FpX_mulu(GEN y,ulong x,GEN p)311 FpX_mulu(GEN y, ulong x,GEN p)
312 {
313 GEN z;
314 long i, l;
315 x = umodui(x, p);
316 if (!x) return zeropol(varn(y));
317 z = cgetg_copy(y, &l); z[1] = y[1];
318 for(i=2; i<l; i++) gel(z,i) = Fp_mulu(gel(y,i), x, p);
319 return z;
320 }
321
322 GEN
FpX_divu(GEN y,ulong x,GEN p)323 FpX_divu(GEN y, ulong x, GEN p)
324 {
325 return FpX_Fp_div(y, utoi(umodui(x, p)), p);
326 }
327
328 GEN
FpX_Fp_mulspec(GEN y,GEN x,GEN p,long ly)329 FpX_Fp_mulspec(GEN y,GEN x,GEN p,long ly)
330 {
331 GEN z;
332 long i;
333 if (!signe(x)) return pol_0(0);
334 z = cgetg(ly+2,t_POL); z[1] = evalsigne(1);
335 for(i=0; i<ly; i++) gel(z,i+2) = Fp_mul(gel(y,i), x, p);
336 return ZX_renormalize(z, ly+2);
337 }
338
339 GEN
FpX_Fp_mul(GEN y,GEN x,GEN p)340 FpX_Fp_mul(GEN y,GEN x,GEN p)
341 {
342 GEN z = FpX_Fp_mulspec(y+2,x,p,lgpol(y));
343 setvarn(z, varn(y)); return z;
344 }
345
346 GEN
FpX_Fp_div(GEN y,GEN x,GEN p)347 FpX_Fp_div(GEN y, GEN x, GEN p)
348 {
349 return FpX_Fp_mul(y, Fp_inv(x, p), p);
350 }
351
352 GEN
FpX_Fp_mul_to_monic(GEN y,GEN x,GEN p)353 FpX_Fp_mul_to_monic(GEN y,GEN x,GEN p)
354 {
355 GEN z;
356 long i, l;
357 z = cgetg_copy(y, &l); z[1] = y[1];
358 for(i=2; i<l-1; i++) gel(z,i) = Fp_mul(gel(y,i), x, p);
359 gel(z,l-1) = gen_1; return z;
360 }
361
362 struct _FpXQ {
363 GEN T, p, aut;
364 };
365
366 static GEN
_FpX_sqr(void * data,GEN x)367 _FpX_sqr(void *data, GEN x)
368 {
369 struct _FpXQ *D = (struct _FpXQ*)data;
370 return FpX_sqr(x, D->p);
371 }
372 static GEN
_FpX_mul(void * data,GEN x,GEN y)373 _FpX_mul(void *data, GEN x, GEN y)
374 {
375 struct _FpXQ *D = (struct _FpXQ*)data;
376 return FpX_mul(x,y, D->p);
377 }
378
379 GEN
FpX_powu(GEN x,ulong n,GEN p)380 FpX_powu(GEN x, ulong n, GEN p)
381 {
382 struct _FpXQ D;
383 if (n==0) return pol_1(varn(x));
384 D.p = p;
385 return gen_powu(x, n, (void *)&D, _FpX_sqr, _FpX_mul);
386 }
387
388 GEN
FpX_halve(GEN y,GEN p)389 FpX_halve(GEN y, GEN p)
390 {
391 GEN z;
392 long i, l;
393 z = cgetg_copy(y, &l); z[1] = y[1];
394 for(i=2; i<l; i++) gel(z,i) = Fp_halve(gel(y,i), p);
395 return z;
396 }
397
398 static GEN
FpX_divrem_basecase(GEN x,GEN y,GEN p,GEN * pr)399 FpX_divrem_basecase(GEN x, GEN y, GEN p, GEN *pr)
400 {
401 long vx, dx, dy, dy1, dz, i, j, sx, lr;
402 pari_sp av0, av;
403 GEN z,p1,rem,lead;
404
405 if (!signe(y)) pari_err_INV("FpX_divrem",y);
406 vx = varn(x);
407 dy = degpol(y);
408 dx = degpol(x);
409 if (dx < dy)
410 {
411 if (pr)
412 {
413 av0 = avma; x = FpX_red(x, p);
414 if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: pol_0(vx); }
415 if (pr == ONLY_REM) return x;
416 *pr = x;
417 }
418 return pol_0(vx);
419 }
420 lead = leading_coeff(y);
421 if (!dy) /* y is constant */
422 {
423 if (pr && pr != ONLY_DIVIDES)
424 {
425 if (pr == ONLY_REM) return pol_0(vx);
426 *pr = pol_0(vx);
427 }
428 av0 = avma;
429 if (equali1(lead)) return FpX_red(x, p);
430 else return gerepileupto(av0, FpX_Fp_div(x, lead, p));
431 }
432 av0 = avma; dz = dx-dy;
433 if (lgefint(p) == 3)
434 { /* assume ab != 0 mod p */
435 ulong pp = to_Flx(&x, &y, p);
436 z = Flx_divrem(x, y, pp, pr);
437 set_avma(av0); /* HACK: assume pr last on stack, then z */
438 if (!z) return NULL;
439 z = leafcopy(z);
440 if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
441 {
442 *pr = leafcopy(*pr);
443 *pr = Flx_to_ZX_inplace(*pr);
444 }
445 return Flx_to_ZX_inplace(z);
446 }
447 lead = equali1(lead)? NULL: gclone(Fp_inv(lead,p));
448 set_avma(av0);
449 z=cgetg(dz+3,t_POL); z[1] = x[1];
450 x += 2; y += 2; z += 2;
451 for (dy1=dy-1; dy1>=0 && !signe(gel(y, dy1)); dy1--);
452
453 p1 = gel(x,dx); av = avma;
454 gel(z,dz) = lead? gerepileuptoint(av, Fp_mul(p1,lead, p)): icopy(p1);
455 for (i=dx-1; i>=dy; i--)
456 {
457 av=avma; p1=gel(x,i);
458 for (j=i-dy1; j<=i && j<=dz; j++)
459 p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
460 if (lead) p1 = mulii(p1,lead);
461 gel(z,i-dy) = gerepileuptoint(av,modii(p1, p));
462 }
463 if (!pr) { guncloneNULL(lead); return z-2; }
464
465 rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
466 for (sx=0; ; i--)
467 {
468 p1 = gel(x,i);
469 for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
470 p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
471 p1 = modii(p1,p); if (signe(p1)) { sx = 1; break; }
472 if (!i) break;
473 set_avma(av);
474 }
475 if (pr == ONLY_DIVIDES)
476 {
477 guncloneNULL(lead);
478 if (sx) return gc_NULL(av0);
479 return gc_const((pari_sp)rem, z-2);
480 }
481 lr=i+3; rem -= lr;
482 rem[0] = evaltyp(t_POL) | evallg(lr);
483 rem[1] = z[-1];
484 p1 = gerepileuptoint((pari_sp)rem, p1);
485 rem += 2; gel(rem,i) = p1;
486 for (i--; i>=0; i--)
487 {
488 av=avma; p1 = gel(x,i);
489 for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
490 p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
491 gel(rem,i) = gerepileuptoint(av, modii(p1,p));
492 }
493 rem -= 2;
494 guncloneNULL(lead);
495 if (!sx) (void)FpX_renormalize(rem, lr);
496 if (pr == ONLY_REM) return gerepileupto(av0,rem);
497 *pr = rem; return z-2;
498 }
499
500 GEN
FpX_div_by_X_x(GEN a,GEN x,GEN p,GEN * r)501 FpX_div_by_X_x(GEN a, GEN x, GEN p, GEN *r)
502 {
503 long l = lg(a), i;
504 GEN z;
505 if (l <= 3)
506 {
507 if (r) *r = l == 2? gen_0: icopy(gel(a,2));
508 return pol_0(0);
509 }
510 l--; z = cgetg(l, t_POL); z[1] = evalsigne(1) | evalvarn(0);
511 gel(z, l-1) = gel(a,l);
512 for (i = l-2; i > 1; i--) /* z[i] = a[i+1] + x*z[i+1] */
513 gel(z,i) = Fp_addmul(gel(a,i+1), x, gel(z,i+1), p);
514 if (r) *r = Fp_addmul(gel(a,2), x, gel(z,2), p);
515 return z;
516 }
517
518 static GEN
_FpX_divrem(void * E,GEN x,GEN y,GEN * r)519 _FpX_divrem(void * E, GEN x, GEN y, GEN *r)
520 {
521 struct _FpXQ *D = (struct _FpXQ*) E;
522 return FpX_divrem(x, y, D->p, r);
523 }
524 static GEN
_FpX_add(void * E,GEN x,GEN y)525 _FpX_add(void * E, GEN x, GEN y) {
526 struct _FpXQ *D = (struct _FpXQ*) E;
527 return FpX_add(x, y, D->p);
528 }
529
530 static struct bb_ring FpX_ring = { _FpX_add,_FpX_mul,_FpX_sqr };
531
532 GEN
FpX_digits(GEN x,GEN T,GEN p)533 FpX_digits(GEN x, GEN T, GEN p)
534 {
535 pari_sp av = avma;
536 struct _FpXQ D;
537 long d = degpol(T), n = (lgpol(x)+d-1)/d;
538 GEN z;
539 D.p = p;
540 z = gen_digits(x,T,n,(void *)&D, &FpX_ring, _FpX_divrem);
541 return gerepileupto(av, z);
542 }
543
544 GEN
FpXV_FpX_fromdigits(GEN x,GEN T,GEN p)545 FpXV_FpX_fromdigits(GEN x, GEN T, GEN p)
546 {
547 pari_sp av = avma;
548 struct _FpXQ D;
549 GEN z;
550 D.p = p;
551 z = gen_fromdigits(x,T,(void *)&D, &FpX_ring);
552 return gerepileupto(av, z);
553 }
554
555 long
FpX_valrem(GEN x,GEN t,GEN p,GEN * py)556 FpX_valrem(GEN x, GEN t, GEN p, GEN *py)
557 {
558 pari_sp av=avma;
559 long k;
560 GEN r, y;
561
562 for (k=0; ; k++)
563 {
564 y = FpX_divrem(x, t, p, &r);
565 if (signe(r)) break;
566 x = y;
567 }
568 *py = gerepilecopy(av,x);
569 return k;
570 }
571
572 static GEN
FpX_halfgcd_basecase(GEN a,GEN b,GEN p)573 FpX_halfgcd_basecase(GEN a, GEN b, GEN p)
574 {
575 pari_sp av=avma;
576 GEN u,u1,v,v1;
577 long vx = varn(a);
578 long n = lgpol(a)>>1;
579 u1 = v = pol_0(vx);
580 u = v1 = pol_1(vx);
581 while (lgpol(b)>n)
582 {
583 GEN r, q = FpX_divrem(a,b,p, &r);
584 a = b; b = r; swap(u,u1); swap(v,v1);
585 u1 = FpX_sub(u1, FpX_mul(u, q, p), p);
586 v1 = FpX_sub(v1, FpX_mul(v, q ,p), p);
587 if (gc_needed(av,2))
588 {
589 if (DEBUGMEM>1) pari_warn(warnmem,"FpX_halfgcd (d = %ld)",degpol(b));
590 gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
591 }
592 }
593 return gerepilecopy(av, mkmat2(mkcol2(u,u1), mkcol2(v,v1)));
594 }
595 static GEN
FpX_addmulmul(GEN u,GEN v,GEN x,GEN y,GEN p)596 FpX_addmulmul(GEN u, GEN v, GEN x, GEN y, GEN p)
597 {
598 return FpX_add(FpX_mul(u, x, p),FpX_mul(v, y, p), p);
599 }
600
601 static GEN
FpXM_FpX_mul2(GEN M,GEN x,GEN y,GEN p)602 FpXM_FpX_mul2(GEN M, GEN x, GEN y, GEN p)
603 {
604 GEN res = cgetg(3, t_COL);
605 gel(res, 1) = FpX_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p);
606 gel(res, 2) = FpX_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p);
607 return res;
608 }
609
610 static GEN
FpXM_mul2(GEN A,GEN B,GEN p)611 FpXM_mul2(GEN A, GEN B, GEN p)
612 {
613 GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
614 GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
615 GEN M1 = FpX_mul(FpX_add(A11,A22, p), FpX_add(B11,B22, p), p);
616 GEN M2 = FpX_mul(FpX_add(A21,A22, p), B11, p);
617 GEN M3 = FpX_mul(A11, FpX_sub(B12,B22, p), p);
618 GEN M4 = FpX_mul(A22, FpX_sub(B21,B11, p), p);
619 GEN M5 = FpX_mul(FpX_add(A11,A12, p), B22, p);
620 GEN M6 = FpX_mul(FpX_sub(A21,A11, p), FpX_add(B11,B12, p), p);
621 GEN M7 = FpX_mul(FpX_sub(A12,A22, p), FpX_add(B21,B22, p), p);
622 GEN T1 = FpX_add(M1,M4, p), T2 = FpX_sub(M7,M5, p);
623 GEN T3 = FpX_sub(M1,M2, p), T4 = FpX_add(M3,M6, p);
624 retmkmat2(mkcol2(FpX_add(T1,T2, p), FpX_add(M2,M4, p)),
625 mkcol2(FpX_add(M3,M5, p), FpX_add(T3,T4, p)));
626 }
627
628 /* Return [0,1;1,-q]*M */
629 static GEN
FpX_FpXM_qmul(GEN q,GEN M,GEN p)630 FpX_FpXM_qmul(GEN q, GEN M, GEN p)
631 {
632 GEN u, v, res = cgetg(3, t_MAT);
633 u = FpX_sub(gcoeff(M,1,1), FpX_mul(gcoeff(M,2,1), q, p), p);
634 gel(res,1) = mkcol2(gcoeff(M,2,1), u);
635 v = FpX_sub(gcoeff(M,1,2), FpX_mul(gcoeff(M,2,2), q, p), p);
636 gel(res,2) = mkcol2(gcoeff(M,2,2), v);
637 return res;
638 }
639
640 static GEN
matid2_FpXM(long v)641 matid2_FpXM(long v)
642 {
643 retmkmat2(mkcol2(pol_1(v),pol_0(v)),
644 mkcol2(pol_0(v),pol_1(v)));
645 }
646
647 static GEN
FpX_shift(GEN a,long n)648 FpX_shift(GEN a, long n) { return RgX_shift_shallow(a, n); }
649
650 static GEN
FpX_halfgcd_split(GEN x,GEN y,GEN p)651 FpX_halfgcd_split(GEN x, GEN y, GEN p)
652 {
653 pari_sp av=avma;
654 GEN R, S, V;
655 GEN y1, r, q;
656 long l = lgpol(x), n = l>>1, k;
657 if (lgpol(y)<=n) return matid2_FpXM(varn(x));
658 R = FpX_halfgcd(FpX_shift(x,-n), FpX_shift(y,-n), p);
659 V = FpXM_FpX_mul2(R,x,y,p); y1 = gel(V,2);
660 if (lgpol(y1)<=n) return gerepilecopy(av, R);
661 q = FpX_divrem(gel(V,1), y1, p, &r);
662 k = 2*n-degpol(y1);
663 S = FpX_halfgcd(FpX_shift(y1,-k), FpX_shift(r,-k), p);
664 return gerepileupto(av, FpXM_mul2(S,FpX_FpXM_qmul(q,R,p),p));
665 }
666
667 /* Return M in GL_2(Fp[X]) such that:
668 if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
669 */
670
671 static GEN
FpX_halfgcd_i(GEN x,GEN y,GEN p)672 FpX_halfgcd_i(GEN x, GEN y, GEN p)
673 {
674 if (lg(x)<=FpX_HALFGCD_LIMIT) return FpX_halfgcd_basecase(x,y,p);
675 return FpX_halfgcd_split(x,y,p);
676 }
677
678 GEN
FpX_halfgcd(GEN x,GEN y,GEN p)679 FpX_halfgcd(GEN x, GEN y, GEN p)
680 {
681 pari_sp av = avma;
682 GEN M,q,r;
683 if (lgefint(p)==3)
684 {
685 ulong pp = to_Flx(&x, &y, p);
686 M = FlxM_to_ZXM(Flx_halfgcd(x, y, pp));
687 }
688 else
689 {
690 if (!signe(x))
691 {
692 long v = varn(x);
693 retmkmat2(mkcol2(pol_0(v),pol_1(v)),
694 mkcol2(pol_1(v),pol_0(v)));
695 }
696 if (degpol(y)<degpol(x)) return FpX_halfgcd_i(x,y,p);
697 q = FpX_divrem(y,x,p,&r);
698 M = FpX_halfgcd_i(x,r,p);
699 gcoeff(M,1,1) = FpX_sub(gcoeff(M,1,1), FpX_mul(q, gcoeff(M,1,2), p), p);
700 gcoeff(M,2,1) = FpX_sub(gcoeff(M,2,1), FpX_mul(q, gcoeff(M,2,2), p), p);
701 }
702 return gerepilecopy(av, M);
703 }
704
705 static GEN
FpX_gcd_basecase(GEN a,GEN b,GEN p)706 FpX_gcd_basecase(GEN a, GEN b, GEN p)
707 {
708 pari_sp av = avma, av0=avma;
709 while (signe(b))
710 {
711 GEN c;
712 if (gc_needed(av0,2))
713 {
714 if (DEBUGMEM>1) pari_warn(warnmem,"FpX_gcd (d = %ld)",degpol(b));
715 gerepileall(av0,2, &a,&b);
716 }
717 av = avma; c = FpX_rem(a,b,p); a=b; b=c;
718 }
719 return gc_const(av, a);
720 }
721
722 GEN
FpX_gcd(GEN x,GEN y,GEN p)723 FpX_gcd(GEN x, GEN y, GEN p)
724 {
725 pari_sp av = avma;
726 if (lgefint(p)==3)
727 {
728 ulong pp;
729 (void)new_chunk((lg(x) + lg(y)) << 2); /* scratch space */
730 pp = to_Flx(&x, &y, p);
731 x = Flx_gcd(x, y, pp);
732 set_avma(av); return Flx_to_ZX(x);
733 }
734 x = FpX_red(x, p);
735 y = FpX_red(y, p);
736 if (!signe(x)) return gerepileupto(av, y);
737 while (lg(y)>FpX_GCD_LIMIT)
738 {
739 GEN c;
740 if (lgpol(y)<=(lgpol(x)>>1))
741 {
742 GEN r = FpX_rem(x, y, p);
743 x = y; y = r;
744 }
745 c = FpXM_FpX_mul2(FpX_halfgcd(x,y, p), x, y, p);
746 x = gel(c,1); y = gel(c,2);
747 gerepileall(av,2,&x,&y);
748 }
749 return gerepileupto(av, FpX_gcd_basecase(x,y,p));
750 }
751
752 /* Return NULL if gcd can be computed else return a factor of p */
753 GEN
FpX_gcd_check(GEN x,GEN y,GEN p)754 FpX_gcd_check(GEN x, GEN y, GEN p)
755 {
756 pari_sp av = avma;
757 GEN a,b,c;
758
759 a = FpX_red(x, p);
760 b = FpX_red(y, p);
761 while (signe(b))
762 {
763 GEN g;
764 if (!invmod(leading_coeff(b), p, &g)) return gerepileuptoint(av,g);
765 b = FpX_Fp_mul_to_monic(b, g, p);
766 c = FpX_rem(a, b, p); a = b; b = c;
767 if (gc_needed(av,1))
768 {
769 if (DEBUGMEM>1) pari_warn(warnmem,"FpX_gcd_check (d = %ld)",degpol(b));
770 gerepileall(av,2,&a,&b);
771 }
772 }
773 return gc_NULL(av);
774 }
775
776 static GEN
FpX_extgcd_basecase(GEN a,GEN b,GEN p,GEN * ptu,GEN * ptv)777 FpX_extgcd_basecase(GEN a, GEN b, GEN p, GEN *ptu, GEN *ptv)
778 {
779 pari_sp av=avma;
780 GEN u,v,d,d1,v1;
781 long vx = varn(a);
782 d = a; d1 = b;
783 v = pol_0(vx); v1 = pol_1(vx);
784 while (signe(d1))
785 {
786 GEN r, q = FpX_divrem(d,d1,p, &r);
787 v = FpX_sub(v,FpX_mul(q,v1,p),p);
788 u=v; v=v1; v1=u;
789 u=r; d=d1; d1=u;
790 if (gc_needed(av,2))
791 {
792 if (DEBUGMEM>1) pari_warn(warnmem,"FpX_extgcd (d = %ld)",degpol(d));
793 gerepileall(av,5, &d,&d1,&u,&v,&v1);
794 }
795 }
796 if (ptu) *ptu = FpX_div(FpX_sub(d,FpX_mul(b,v,p),p),a,p);
797 *ptv = v; return d;
798 }
799
800 static GEN
FpX_extgcd_halfgcd(GEN x,GEN y,GEN p,GEN * ptu,GEN * ptv)801 FpX_extgcd_halfgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
802 {
803 pari_sp av=avma;
804 GEN u,v,R = matid2_FpXM(varn(x));
805 while (lg(y)>FpX_EXTGCD_LIMIT)
806 {
807 GEN M, c;
808 if (lgpol(y)<=(lgpol(x)>>1))
809 {
810 GEN r, q = FpX_divrem(x, y, p, &r);
811 x = y; y = r;
812 R = FpX_FpXM_qmul(q, R, p);
813 }
814 M = FpX_halfgcd(x,y, p);
815 c = FpXM_FpX_mul2(M, x,y, p);
816 R = FpXM_mul2(M, R, p);
817 x = gel(c,1); y = gel(c,2);
818 gerepileall(av,3,&x,&y,&R);
819 }
820 y = FpX_extgcd_basecase(x,y,p,&u,&v);
821 if (ptu) *ptu = FpX_addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1),p);
822 *ptv = FpX_addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2),p);
823 return y;
824 }
825
826 /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
827 * ux + vy = gcd (mod p) */
828 GEN
FpX_extgcd(GEN x,GEN y,GEN p,GEN * ptu,GEN * ptv)829 FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv)
830 {
831 GEN d;
832 pari_sp ltop=avma;
833 if (lgefint(p)==3)
834 {
835 ulong pp = to_Flx(&x, &y, p);
836 d = Flx_extgcd(x,y, pp, ptu,ptv);
837 d = Flx_to_ZX(d);
838 if (ptu) *ptu=Flx_to_ZX(*ptu);
839 *ptv=Flx_to_ZX(*ptv);
840 }
841 else
842 {
843 x = FpX_red(x, p);
844 y = FpX_red(y, p);
845 if (lg(y)>FpX_EXTGCD_LIMIT)
846 d = FpX_extgcd_halfgcd(x, y, p, ptu, ptv);
847 else
848 d = FpX_extgcd_basecase(x, y, p, ptu, ptv);
849 }
850 gerepileall(ltop,ptu?3:2,&d,ptv,ptu);
851 return d;
852 }
853
854 GEN
FpX_rescale(GEN P,GEN h,GEN p)855 FpX_rescale(GEN P, GEN h, GEN p)
856 {
857 long i, l = lg(P);
858 GEN Q = cgetg(l,t_POL), hi = h;
859 gel(Q,l-1) = gel(P,l-1);
860 for (i=l-2; i>=2; i--)
861 {
862 gel(Q,i) = Fp_mul(gel(P,i), hi, p);
863 if (i == 2) break;
864 hi = Fp_mul(hi,h, p);
865 }
866 Q[1] = P[1]; return Q;
867 }
868
869 GEN
FpX_deriv(GEN x,GEN p)870 FpX_deriv(GEN x, GEN p) { return FpX_red(ZX_deriv(x), p); }
871
872 /* Compute intformal(x^n*S)/x^(n+1) */
873 static GEN
FpX_integXn(GEN x,long n,GEN p)874 FpX_integXn(GEN x, long n, GEN p)
875 {
876 long i, lx = lg(x);
877 GEN y;
878 if (lx == 2) return ZX_copy(x);
879 y = cgetg(lx, t_POL); y[1] = x[1];
880 for (i=2; i<lx; i++)
881 {
882 GEN xi = gel(x,i);
883 if (!signe(xi))
884 gel(y,i) = gen_0;
885 else
886 {
887 ulong j = n+i-1;
888 ulong d = ugcd(j, umodiu(xi, j));
889 if (d==1)
890 gel(y,i) = Fp_divu(xi, j, p);
891 else
892 gel(y,i) = Fp_divu(diviuexact(xi, d), j/d, p);
893 }
894 }
895 return ZX_renormalize(y, lx);;
896 }
897
898 GEN
FpX_integ(GEN x,GEN p)899 FpX_integ(GEN x, GEN p)
900 {
901 long i, lx = lg(x);
902 GEN y;
903 if (lx == 2) return ZX_copy(x);
904 y = cgetg(lx+1, t_POL); y[1] = x[1];
905 gel(y,2) = gen_0;
906 for (i=3; i<=lx; i++)
907 gel(y,i) = signe(gel(x,i-1))? Fp_divu(gel(x,i-1), i-2, p): gen_0;
908 return ZX_renormalize(y, lx+1);;
909 }
910
911 INLINE GEN
FpXn_recip(GEN P,long n)912 FpXn_recip(GEN P, long n)
913 { return RgXn_recip_shallow(P, n); }
914
915 GEN
FpX_Newton(GEN P,long n,GEN p)916 FpX_Newton(GEN P, long n, GEN p)
917 {
918 pari_sp av = avma;
919 GEN dP = FpX_deriv(P, p);
920 GEN Q = FpXn_recip(FpX_div(FpX_shift(dP,n), P, p), n);
921 return gerepilecopy(av, Q);
922 }
923
924 GEN
FpX_fromNewton(GEN P,GEN p)925 FpX_fromNewton(GEN P, GEN p)
926 {
927 pari_sp av = avma;
928 if (lgefint(p)==3)
929 {
930 ulong pp = p[2];
931 GEN Q = Flx_fromNewton(ZX_to_Flx(P, pp), pp);
932 return gerepileupto(av, Flx_to_ZX(Q));
933 } else
934 {
935 long n = itos(modii(constant_coeff(P), p))+1;
936 GEN z = FpX_neg(FpX_shift(P,-1),p);
937 GEN Q = FpXn_recip(FpXn_expint(z, n, p), n);
938 return gerepilecopy(av, Q);
939 }
940 }
941
942 GEN
FpX_invLaplace(GEN x,GEN p)943 FpX_invLaplace(GEN x, GEN p)
944 {
945 pari_sp av = avma;
946 long i, d = degpol(x);
947 GEN t, y;
948 if (d <= 1) return gcopy(x);
949 t = Fp_inv(factorial_Fp(d, p), p);
950 y = cgetg(d+3, t_POL);
951 y[1] = x[1];
952 for (i=d; i>=2; i--)
953 {
954 gel(y,i+2) = Fp_mul(gel(x,i+2), t, p);
955 t = Fp_mulu(t, i, p);
956 }
957 gel(y,3) = gel(x,3);
958 gel(y,2) = gel(x,2);
959 return gerepilecopy(av, y);
960 }
961
962 GEN
FpX_Laplace(GEN x,GEN p)963 FpX_Laplace(GEN x, GEN p)
964 {
965 pari_sp av = avma;
966 long i, d = degpol(x);
967 GEN t = gen_1;
968 GEN y;
969 if (d <= 1) return gcopy(x);
970 y = cgetg(d+3, t_POL);
971 y[1] = x[1];
972 gel(y,2) = gel(x,2);
973 gel(y,3) = gel(x,3);
974 for (i=2; i<=d; i++)
975 {
976 t = Fp_mulu(t, i, p);
977 gel(y,i+2) = Fp_mul(gel(x,i+2), t, p);
978 }
979 return gerepilecopy(av, y);
980 }
981
982 int
FpX_is_squarefree(GEN f,GEN p)983 FpX_is_squarefree(GEN f, GEN p)
984 {
985 pari_sp av = avma;
986 GEN z = FpX_gcd(f,FpX_deriv(f,p),p);
987 set_avma(av);
988 return degpol(z)==0;
989 }
990
991 GEN
random_FpX(long d1,long v,GEN p)992 random_FpX(long d1, long v, GEN p)
993 {
994 long i, d = d1+2;
995 GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
996 for (i=2; i<d; i++) gel(y,i) = randomi(p);
997 return FpX_renormalize(y,d);
998 }
999
1000 GEN
FpX_dotproduct(GEN x,GEN y,GEN p)1001 FpX_dotproduct(GEN x, GEN y, GEN p)
1002 {
1003 long i, l = minss(lg(x), lg(y));
1004 pari_sp av;
1005 GEN c;
1006 if (l == 2) return gen_0;
1007 av = avma; c = mulii(gel(x,2),gel(y,2));
1008 for (i=3; i<l; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
1009 return gerepileuptoint(av, modii(c,p));
1010 }
1011
1012 /* Evaluation in Fp
1013 * x a ZX and y an Fp, return x(y) mod p
1014 *
1015 * If p is very large (several longs) and x has small coefficients(<<p),
1016 * then Brent & Kung algorithm is faster. */
1017 GEN
FpX_eval(GEN x,GEN y,GEN p)1018 FpX_eval(GEN x,GEN y,GEN p)
1019 {
1020 pari_sp av;
1021 GEN p1,r,res;
1022 long j, i=lg(x)-1;
1023 if (i<=2 || !signe(y))
1024 return (i==1)? gen_0: modii(gel(x,2),p);
1025 res=cgeti(lgefint(p));
1026 av=avma; p1=gel(x,i);
1027 /* specific attention to sparse polynomials (see poleval)*/
1028 /*You've guessed it! It's a copy-paste(tm)*/
1029 for (i--; i>=2; i=j-1)
1030 {
1031 for (j=i; !signe(gel(x,j)); j--)
1032 if (j==2)
1033 {
1034 if (i!=j) y = Fp_powu(y,i-j+1,p);
1035 p1=mulii(p1,y);
1036 goto fppoleval;/*sorry break(2) no implemented*/
1037 }
1038 r = (i==j)? y: Fp_powu(y,i-j+1,p);
1039 p1 = Fp_addmul(gel(x,j), p1, r, p);
1040 if ((i & 7) == 0) { affii(p1, res); p1 = res; set_avma(av); }
1041 }
1042 fppoleval:
1043 modiiz(p1,p,res); return gc_const(av, res);
1044 }
1045
1046 /* Tz=Tx*Ty where Tx and Ty coprime
1047 * return lift(chinese(Mod(x*Mod(1,p),Tx*Mod(1,p)),Mod(y*Mod(1,p),Ty*Mod(1,p))))
1048 * if Tz is NULL it is computed
1049 * As we do not return it, and the caller will frequently need it,
1050 * it must compute it and pass it.
1051 */
1052 GEN
FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p)1053 FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p)
1054 {
1055 pari_sp av = avma;
1056 GEN ax,p1;
1057 ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p);
1058 p1 = FpX_mul(ax, FpX_sub(y,x,p),p);
1059 p1 = FpX_add(x,p1,p);
1060 if (!Tz) Tz=FpX_mul(Tx,Ty,p);
1061 p1 = FpX_rem(p1,Tz,p);
1062 return gerepileupto(av,p1);
1063 }
1064
1065 /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab), with R=A%B, a=deg(A) ...*/
1066 GEN
FpX_resultant(GEN a,GEN b,GEN p)1067 FpX_resultant(GEN a, GEN b, GEN p)
1068 {
1069 long da,db,dc;
1070 pari_sp av;
1071 GEN c,lb, res = gen_1;
1072
1073 if (!signe(a) || !signe(b)) return gen_0;
1074 if (lgefint(p) == 3)
1075 {
1076 pari_sp av = avma;
1077 ulong pp = to_Flx(&a, &b, p);
1078 long r = Flx_resultant(a, b, pp);
1079 set_avma(av);
1080 return utoi(r);
1081 }
1082
1083 da = degpol(a);
1084 db = degpol(b);
1085 if (db > da)
1086 {
1087 swapspec(a,b, da,db);
1088 if (both_odd(da,db)) res = subii(p, res);
1089 }
1090 if (!da) return gen_1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
1091 av = avma;
1092 while (db)
1093 {
1094 lb = gel(b,db+2);
1095 c = FpX_rem(a,b, p);
1096 a = b; b = c; dc = degpol(c);
1097 if (dc < 0) return gc_const(av, gen_0);
1098
1099 if (both_odd(da,db)) res = subii(p, res);
1100 if (!equali1(lb)) res = Fp_mul(res, Fp_powu(lb, da - dc, p), p);
1101 if (gc_needed(av,2))
1102 {
1103 if (DEBUGMEM>1) pari_warn(warnmem,"FpX_resultant (da = %ld)",da);
1104 gerepileall(av,3, &a,&b,&res);
1105 }
1106 da = db; /* = degpol(a) */
1107 db = dc; /* = degpol(b) */
1108 }
1109 res = Fp_mul(res, Fp_powu(gel(b,2), da, p), p);
1110 return gerepileuptoint(av, res);
1111 }
1112
1113 /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
1114 GEN
FpX_disc(GEN P,GEN p)1115 FpX_disc(GEN P, GEN p)
1116 {
1117 pari_sp av = avma;
1118 GEN L, dP = FpX_deriv(P,p), D = FpX_resultant(P, dP, p);
1119 long dd;
1120 if (!signe(D)) return gen_0;
1121 dd = degpol(P) - 2 - degpol(dP); /* >= -1; > -1 iff p | deg(P) */
1122 L = leading_coeff(P);
1123 if (dd && !equali1(L))
1124 D = (dd == -1)? Fp_div(D,L,p): Fp_mul(D, Fp_powu(L, dd, p), p);
1125 if (degpol(P) & 2) D = Fp_neg(D ,p);
1126 return gerepileuptoint(av, D);
1127 }
1128
1129 GEN
FpXV_prod(GEN V,GEN p)1130 FpXV_prod(GEN V, GEN p)
1131 {
1132 struct _FpXQ D;
1133 D.p = p;
1134 return gen_product(V, (void *)&D, &_FpX_mul);
1135 }
1136
1137 GEN
FpV_roots_to_pol(GEN V,GEN p,long v)1138 FpV_roots_to_pol(GEN V, GEN p, long v)
1139 {
1140 pari_sp ltop=avma;
1141 long i;
1142 GEN g=cgetg(lg(V),t_VEC);
1143 for(i=1;i<lg(V);i++)
1144 gel(g,i) = deg1pol_shallow(gen_1,modii(negi(gel(V,i)),p),v);
1145 return gerepileupto(ltop,FpXV_prod(g,p));
1146 }
1147
1148 /* invert all elements of x mod p using Montgomery's multi-inverse trick.
1149 * Not stack-clean. */
1150 GEN
FpV_inv(GEN x,GEN p)1151 FpV_inv(GEN x, GEN p)
1152 {
1153 long i, lx = lg(x);
1154 GEN u, y = cgetg(lx, t_VEC);
1155
1156 gel(y,1) = gel(x,1);
1157 for (i=2; i<lx; i++) gel(y,i) = Fp_mul(gel(y,i-1), gel(x,i), p);
1158
1159 u = Fp_inv(gel(y,--i), p);
1160 for ( ; i > 1; i--)
1161 {
1162 gel(y,i) = Fp_mul(u, gel(y,i-1), p);
1163 u = Fp_mul(u, gel(x,i), p); /* u = 1 / (x[1] ... x[i-1]) */
1164 }
1165 gel(y,1) = u; return y;
1166 }
1167 GEN
FqV_inv(GEN x,GEN T,GEN p)1168 FqV_inv(GEN x, GEN T, GEN p)
1169 {
1170 long i, lx = lg(x);
1171 GEN u, y = cgetg(lx, t_VEC);
1172
1173 gel(y,1) = gel(x,1);
1174 for (i=2; i<lx; i++) gel(y,i) = Fq_mul(gel(y,i-1), gel(x,i), T,p);
1175
1176 u = Fq_inv(gel(y,--i), T,p);
1177 for ( ; i > 1; i--)
1178 {
1179 gel(y,i) = Fq_mul(u, gel(y,i-1), T,p);
1180 u = Fq_mul(u, gel(x,i), T,p); /* u = 1 / (x[1] ... x[i-1]) */
1181 }
1182 gel(y,1) = u; return y;
1183 }
1184
1185 /***********************************************************************/
1186 /** **/
1187 /** Barrett reduction **/
1188 /** **/
1189 /***********************************************************************/
1190
1191 static GEN
FpX_invBarrett_basecase(GEN T,GEN p)1192 FpX_invBarrett_basecase(GEN T, GEN p)
1193 {
1194 long i, l=lg(T)-1, lr = l-1, k;
1195 GEN r=cgetg(lr, t_POL); r[1]=T[1];
1196 gel(r,2) = gen_1;
1197 for (i=3; i<lr; i++)
1198 {
1199 pari_sp av = avma;
1200 GEN u = gel(T,l-i+2);
1201 for (k=3; k<i; k++)
1202 u = addii(u, mulii(gel(T,l-i+k), gel(r,k)));
1203 gel(r,i) = gerepileupto(av, modii(negi(u), p));
1204 }
1205 return FpX_renormalize(r,lr);
1206 }
1207
1208 /* Return new lgpol */
1209 static long
ZX_lgrenormalizespec(GEN x,long lx)1210 ZX_lgrenormalizespec(GEN x, long lx)
1211 {
1212 long i;
1213 for (i = lx-1; i>=0; i--)
1214 if (signe(gel(x,i))) break;
1215 return i+1;
1216 }
1217
1218 INLINE GEN
FpX_recipspec(GEN x,long l,long n)1219 FpX_recipspec(GEN x, long l, long n)
1220 {
1221 return RgX_recipspec_shallow(x, l, n);
1222 }
1223
1224 static GEN
FpX_invBarrett_Newton(GEN T,GEN p)1225 FpX_invBarrett_Newton(GEN T, GEN p)
1226 {
1227 pari_sp av = avma;
1228 long nold, lx, lz, lq, l = degpol(T), i, lQ;
1229 GEN q, y, z, x = cgetg(l+2, t_POL) + 2;
1230 ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1231 for (i=0;i<l;i++) gel(x,i) = gen_0;
1232 q = FpX_recipspec(T+2,l+1,l+1); lQ = lgpol(q); q+=2;
1233 /* We work on _spec_ FpX's, all the l[xzq] below are lgpol's */
1234
1235 /* initialize */
1236 gel(x,0) = Fp_inv(gel(q,0), p);
1237 if (lQ>1) gel(q,1) = Fp_red(gel(q,1), p);
1238 if (lQ>1 && signe(gel(q,1)))
1239 {
1240 GEN u = gel(q, 1);
1241 if (!equali1(gel(x,0))) u = Fp_mul(u, Fp_sqr(gel(x,0), p), p);
1242 gel(x,1) = Fp_neg(u, p); lx = 2;
1243 }
1244 else
1245 lx = 1;
1246 nold = 1;
1247 for (; mask > 1; )
1248 { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1249 long i, lnew, nnew = nold << 1;
1250
1251 if (mask & 1) nnew--;
1252 mask >>= 1;
1253
1254 lnew = nnew + 1;
1255 lq = ZX_lgrenormalizespec(q, minss(lQ,lnew));
1256 z = FpX_mulspec(x, q, p, lx, lq); /* FIXME: high product */
1257 lz = lgpol(z); if (lz > lnew) lz = lnew;
1258 z += 2;
1259 /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1260 for (i = nold; i < lz; i++) if (signe(gel(z,i))) break;
1261 nold = nnew;
1262 if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1263
1264 /* z + i represents (x*q - 1) / t^i */
1265 lz = ZX_lgrenormalizespec (z+i, lz-i);
1266 z = FpX_mulspec(x, z+i, p, lx, lz); /* FIXME: low product */
1267 lz = lgpol(z); z += 2;
1268 if (lz > lnew-i) lz = ZX_lgrenormalizespec(z, lnew-i);
1269
1270 lx = lz+ i;
1271 y = x + i; /* x -= z * t^i, in place */
1272 for (i = 0; i < lz; i++) gel(y,i) = Fp_neg(gel(z,i), p);
1273 }
1274 x -= 2; setlg(x, lx + 2); x[1] = T[1];
1275 return gerepilecopy(av, x);
1276 }
1277
1278 /* 1/polrecip(T)+O(x^(deg(T)-1)) */
1279 GEN
FpX_invBarrett(GEN T,GEN p)1280 FpX_invBarrett(GEN T, GEN p)
1281 {
1282 pari_sp ltop = avma;
1283 long l = lg(T);
1284 GEN r;
1285 if (l<5) return pol_0(varn(T));
1286 if (l<=FpX_INVBARRETT_LIMIT)
1287 {
1288 GEN c = gel(T,l-1), ci=gen_1;
1289 if (!equali1(c))
1290 {
1291 ci = Fp_inv(c, p);
1292 T = FpX_Fp_mul(T, ci, p);
1293 r = FpX_invBarrett_basecase(T, p);
1294 r = FpX_Fp_mul(r, ci, p);
1295 } else
1296 r = FpX_invBarrett_basecase(T, p);
1297 }
1298 else
1299 r = FpX_invBarrett_Newton(T, p);
1300 return gerepileupto(ltop, r);
1301 }
1302
1303 GEN
FpX_get_red(GEN T,GEN p)1304 FpX_get_red(GEN T, GEN p)
1305 {
1306 if (typ(T)==t_POL && lg(T)>FpX_BARRETT_LIMIT)
1307 retmkvec2(FpX_invBarrett(T,p),T);
1308 return T;
1309 }
1310
1311 /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1312 * and mg is the Barrett inverse of T. */
1313 static GEN
FpX_divrem_Barrettspec(GEN x,long l,GEN mg,GEN T,GEN p,GEN * pr)1314 FpX_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, GEN p, GEN *pr)
1315 {
1316 GEN q, r;
1317 long lt = degpol(T); /*We discard the leading term*/
1318 long ld, lm, lT, lmg;
1319 ld = l-lt;
1320 lm = minss(ld, lgpol(mg));
1321 lT = ZX_lgrenormalizespec(T+2,lt);
1322 lmg = ZX_lgrenormalizespec(mg+2,lm);
1323 q = FpX_recipspec(x+lt,ld,ld); /* q = rec(x) lq<=ld*/
1324 q = FpX_mulspec(q+2,mg+2,p,lgpol(q),lmg); /* q = rec(x) * mg lq<=ld+lm*/
1325 q = FpX_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lq<=ld*/
1326 if (!pr) return q;
1327 r = FpX_mulspec(q+2,T+2,p,lgpol(q),lT); /* r = q*pol lr<=ld+lt*/
1328 r = FpX_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - r lr<=lt */
1329 if (pr == ONLY_REM) return r;
1330 *pr = r; return q;
1331 }
1332
1333 static GEN
FpX_divrem_Barrett(GEN x,GEN mg,GEN T,GEN p,GEN * pr)1334 FpX_divrem_Barrett(GEN x, GEN mg, GEN T, GEN p, GEN *pr)
1335 {
1336 GEN q = NULL, r = FpX_red(x, p);
1337 long l = lgpol(r), lt = degpol(T), lm = 2*lt-1, v = varn(T);
1338 long i;
1339 if (l <= lt)
1340 {
1341 if (pr == ONLY_REM) return r;
1342 if (pr == ONLY_DIVIDES) return signe(r)? NULL: pol_0(v);
1343 if (pr) *pr = r;
1344 return pol_0(v);
1345 }
1346 if (lt <= 1)
1347 return FpX_divrem_basecase(r,T,p,pr);
1348 if (pr != ONLY_REM && l>lm)
1349 {
1350 q = cgetg(l-lt+2, t_POL); q[1] = T[1];
1351 for (i=0;i<l-lt;i++) gel(q+2,i) = gen_0;
1352 }
1353 while (l>lm)
1354 {
1355 GEN zr, zq = FpX_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,&zr);
1356 long lz = lgpol(zr);
1357 if (pr != ONLY_REM)
1358 {
1359 long lq = lgpol(zq);
1360 for(i=0; i<lq; i++) gel(q+2+l-lm,i) = gel(zq,2+i);
1361 }
1362 for(i=0; i<lz; i++) gel(r+2+l-lm,i) = gel(zr,2+i);
1363 l = l-lm+lz;
1364 }
1365 if (pr == ONLY_REM)
1366 {
1367 if (l > lt)
1368 r = FpX_divrem_Barrettspec(r+2, l, mg, T, p, ONLY_REM);
1369 else
1370 r = FpX_renormalize(r, l+2);
1371 setvarn(r, v); return r;
1372 }
1373 if (l > lt)
1374 {
1375 GEN zq = FpX_divrem_Barrettspec(r+2,l,mg,T,p, pr? &r: NULL);
1376 if (!q) q = zq;
1377 else
1378 {
1379 long lq = lgpol(zq);
1380 for(i=0; i<lq; i++) gel(q+2,i) = gel(zq,2+i);
1381 }
1382 }
1383 else if (pr)
1384 r = FpX_renormalize(r, l+2);
1385 setvarn(q, v); q = FpX_renormalize(q, lg(q));
1386 if (pr == ONLY_DIVIDES) return signe(r)? NULL: q;
1387 if (pr) { setvarn(r, v); *pr = r; }
1388 return q;
1389 }
1390
1391 GEN
FpX_divrem(GEN x,GEN T,GEN p,GEN * pr)1392 FpX_divrem(GEN x, GEN T, GEN p, GEN *pr)
1393 {
1394 GEN B, y;
1395 long dy, dx, d;
1396 if (pr==ONLY_REM) return FpX_rem(x, T, p);
1397 y = get_FpX_red(T, &B);
1398 dy = degpol(y); dx = degpol(x); d = dx-dy;
1399 if (!B && d+3 < FpX_DIVREM_BARRETT_LIMIT)
1400 return FpX_divrem_basecase(x,y,p,pr);
1401 else if (lgefint(p)==3)
1402 {
1403 pari_sp av = avma;
1404 ulong pp = to_Flxq(&x, &T, p);
1405 GEN z = Flx_divrem(x, T, pp, pr);
1406 if (!z) return NULL;
1407 if (!pr || pr == ONLY_DIVIDES)
1408 return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
1409 z = Flx_to_ZX(z);
1410 *pr = Flx_to_ZX(*pr);
1411 gerepileall(av, 2, &z, pr);
1412 return z;
1413 } else
1414 {
1415 pari_sp av=avma;
1416 GEN mg = B? B: FpX_invBarrett(y, p);
1417 GEN q1 = FpX_divrem_Barrett(x,mg,y,p,pr);
1418 if (!q1) return gc_NULL(av);
1419 if (!pr || pr==ONLY_DIVIDES) return gerepilecopy(av, q1);
1420 gerepileall(av,2,&q1,pr);
1421 return q1;
1422 }
1423 }
1424
1425 GEN
FpX_rem(GEN x,GEN T,GEN p)1426 FpX_rem(GEN x, GEN T, GEN p)
1427 {
1428 GEN B, y = get_FpX_red(T, &B);
1429 long dy = degpol(y), dx = degpol(x), d = dx-dy;
1430 if (d < 0) return FpX_red(x,p);
1431 if (!B && d+3 < FpX_REM_BARRETT_LIMIT)
1432 return FpX_divrem_basecase(x,y,p,ONLY_REM);
1433 else if (lgefint(p)==3)
1434 {
1435 pari_sp av = avma;
1436 ulong pp = to_Flxq(&x, &T, p);
1437 return Flx_to_ZX_inplace(gerepileuptoleaf(av, Flx_rem(x, T, pp)));
1438 } else
1439 {
1440 pari_sp av = avma;
1441 GEN mg = B? B: FpX_invBarrett(y, p);
1442 return gerepileupto(av, FpX_divrem_Barrett(x, mg, y, p, ONLY_REM));
1443 }
1444 }
1445
1446 static GEN
FpXV_producttree_dbl(GEN t,long n,GEN p)1447 FpXV_producttree_dbl(GEN t, long n, GEN p)
1448 {
1449 long i, j, k, m = n==1 ? 1: expu(n-1)+1;
1450 GEN T = cgetg(m+1, t_VEC);
1451 gel(T,1) = t;
1452 for (i=2; i<=m; i++)
1453 {
1454 GEN u = gel(T, i-1);
1455 long n = lg(u)-1;
1456 GEN t = cgetg(((n+1)>>1)+1, t_VEC);
1457 for (j=1, k=1; k<n; j++, k+=2)
1458 gel(t, j) = FpX_mul(gel(u, k), gel(u, k+1), p);
1459 gel(T, i) = t;
1460 }
1461 return T;
1462 }
1463
1464 static GEN
FpV_producttree(GEN xa,GEN s,GEN p,long vs)1465 FpV_producttree(GEN xa, GEN s, GEN p, long vs)
1466 {
1467 long n = lg(xa)-1;
1468 long j, k, ls = lg(s);
1469 GEN t = cgetg(ls, t_VEC);
1470 for (j=1, k=1; j<ls; k+=s[j++])
1471 gel(t, j) = s[j] == 1 ?
1472 deg1pol(gen_1, Fp_neg(gel(xa,k), p), vs):
1473 deg2pol_shallow(gen_1,
1474 Fp_neg(Fp_add(gel(xa,k), gel(xa,k+1), p), p),
1475 Fp_mul(gel(xa,k), gel(xa,k+1), p), vs);
1476 return FpXV_producttree_dbl(t, n, p);
1477 }
1478
1479 static GEN
FpX_FpXV_multirem_dbl_tree(GEN P,GEN T,GEN p)1480 FpX_FpXV_multirem_dbl_tree(GEN P, GEN T, GEN p)
1481 {
1482 long i,j,k;
1483 long m = lg(T)-1;
1484 GEN t;
1485 GEN Tp = cgetg(m+1, t_VEC);
1486 gel(Tp, m) = mkvec(P);
1487 for (i=m-1; i>=1; i--)
1488 {
1489 GEN u = gel(T, i);
1490 GEN v = gel(Tp, i+1);
1491 long n = lg(u)-1;
1492 t = cgetg(n+1, t_VEC);
1493 for (j=1, k=1; k<n; j++, k+=2)
1494 {
1495 gel(t, k) = FpX_rem(gel(v, j), gel(u, k), p);
1496 gel(t, k+1) = FpX_rem(gel(v, j), gel(u, k+1), p);
1497 }
1498 gel(Tp, i) = t;
1499 }
1500 return Tp;
1501 }
1502
1503 static GEN
FpX_FpV_multieval_tree(GEN P,GEN xa,GEN T,GEN p)1504 FpX_FpV_multieval_tree(GEN P, GEN xa, GEN T, GEN p)
1505 {
1506 pari_sp av = avma;
1507 long j,k;
1508 GEN Tp = FpX_FpXV_multirem_dbl_tree(P, T, p);
1509 GEN R = cgetg(lg(xa), t_VEC);
1510 GEN u = gel(T, 1);
1511 GEN v = gel(Tp, 1);
1512 long n = lg(u)-1;
1513 for (j=1, k=1; j<=n; j++)
1514 {
1515 long c, d = degpol(gel(u,j));
1516 for (c=1; c<=d; c++, k++)
1517 gel(R,k) = FpX_eval(gel(v, j), gel(xa,k), p);
1518 }
1519 return gerepileupto(av, R);
1520 }
1521
1522 static GEN
FpVV_polint_tree(GEN T,GEN R,GEN s,GEN xa,GEN ya,GEN p,long vs)1523 FpVV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, GEN p, long vs)
1524 {
1525 pari_sp av = avma;
1526 long m = lg(T)-1;
1527 long i, j, k, ls = lg(s);
1528 GEN Tp = cgetg(m+1, t_VEC);
1529 GEN t = cgetg(ls, t_VEC);
1530 for (j=1, k=1; j<ls; k+=s[j++])
1531 if (s[j]==2)
1532 {
1533 GEN a = Fp_mul(gel(ya,k), gel(R,k), p);
1534 GEN b = Fp_mul(gel(ya,k+1), gel(R,k+1), p);
1535 gel(t, j) = deg1pol(Fp_add(a, b, p),
1536 Fp_neg(Fp_add(Fp_mul(gel(xa,k), b, p ),
1537 Fp_mul(gel(xa,k+1), a, p), p), p), vs);
1538 }
1539 else
1540 gel(t, j) = scalarpol(Fp_mul(gel(ya,k), gel(R,k), p), vs);
1541 gel(Tp, 1) = t;
1542 for (i=2; i<=m; i++)
1543 {
1544 GEN u = gel(T, i-1);
1545 GEN t = cgetg(lg(gel(T,i)), t_VEC);
1546 GEN v = gel(Tp, i-1);
1547 long n = lg(v)-1;
1548 for (j=1, k=1; k<n; j++, k+=2)
1549 gel(t, j) = FpX_add(ZX_mul(gel(u, k), gel(v, k+1)),
1550 ZX_mul(gel(u, k+1), gel(v, k)), p);
1551 gel(Tp, i) = t;
1552 }
1553 return gerepilecopy(av, gmael(Tp,m,1));
1554 }
1555
1556 GEN
FpX_FpV_multieval(GEN P,GEN xa,GEN p)1557 FpX_FpV_multieval(GEN P, GEN xa, GEN p)
1558 {
1559 pari_sp av = avma;
1560 GEN s = producttree_scheme(lg(xa)-1);
1561 GEN T = FpV_producttree(xa, s, p, varn(P));
1562 return gerepileupto(av, FpX_FpV_multieval_tree(P, xa, T, p));
1563 }
1564
1565 GEN
FpV_polint(GEN xa,GEN ya,GEN p,long vs)1566 FpV_polint(GEN xa, GEN ya, GEN p, long vs)
1567 {
1568 pari_sp av = avma;
1569 GEN s, T, P, R;
1570 long m;
1571 if (lgefint(p) == 3)
1572 {
1573 ulong pp = p[2];
1574 P = Flv_polint(ZV_to_Flv(xa, pp), ZV_to_Flv(ya, pp), pp, evalvarn(vs));
1575 return gerepileupto(av, Flx_to_ZX(P));
1576 }
1577 s = producttree_scheme(lg(xa)-1);
1578 T = FpV_producttree(xa, s, p, vs);
1579 m = lg(T)-1;
1580 P = FpX_deriv(gmael(T, m, 1), p);
1581 R = FpV_inv(FpX_FpV_multieval_tree(P, xa, T, p), p);
1582 return gerepileupto(av, FpVV_polint_tree(T, R, s, xa, ya, p, vs));
1583 }
1584
1585 GEN
FpV_FpM_polint(GEN xa,GEN ya,GEN p,long vs)1586 FpV_FpM_polint(GEN xa, GEN ya, GEN p, long vs)
1587 {
1588 pari_sp av = avma;
1589 GEN s = producttree_scheme(lg(xa)-1);
1590 GEN T = FpV_producttree(xa, s, p, vs);
1591 long i, m = lg(T)-1, l = lg(ya)-1;
1592 GEN P = FpX_deriv(gmael(T, m, 1), p);
1593 GEN R = FpV_inv(FpX_FpV_multieval_tree(P, xa, T, p), p);
1594 GEN M = cgetg(l+1, t_VEC);
1595 for (i=1; i<=l; i++)
1596 gel(M,i) = FpVV_polint_tree(T, R, s, xa, gel(ya,i), p, vs);
1597 return gerepileupto(av, M);
1598 }
1599
1600 GEN
FpV_invVandermonde(GEN L,GEN den,GEN p)1601 FpV_invVandermonde(GEN L, GEN den, GEN p)
1602 {
1603 pari_sp av = avma;
1604 long i, n = lg(L);
1605 GEN M, R;
1606 GEN s = producttree_scheme(n-1);
1607 GEN tree = FpV_producttree(L, s, p, 0);
1608 long m = lg(tree)-1;
1609 GEN T = gmael(tree, m, 1);
1610 R = FpV_inv(FpX_FpV_multieval_tree(FpX_deriv(T, p), L, tree, p), p);
1611 if (den) R = FpC_Fp_mul(R, den, p);
1612 M = cgetg(n, t_MAT);
1613 for (i = 1; i < n; i++)
1614 {
1615 GEN P = FpX_Fp_mul(FpX_div_by_X_x(T, gel(L,i), p, NULL), gel(R,i), p);
1616 gel(M,i) = RgX_to_RgC(P, n-1);
1617 }
1618 return gerepilecopy(av, M);
1619 }
1620
1621 static GEN
FpXV_producttree(GEN xa,GEN s,GEN p)1622 FpXV_producttree(GEN xa, GEN s, GEN p)
1623 {
1624 long n = lg(xa)-1;
1625 long j, k, ls = lg(s);
1626 GEN t = cgetg(ls, t_VEC);
1627 for (j=1, k=1; j<ls; k+=s[j++])
1628 gel(t, j) = s[j] == 1 ?
1629 gel(xa,k): FpX_mul(gel(xa,k),gel(xa,k+1),p);
1630 return FpXV_producttree_dbl(t, n, p);
1631 }
1632
1633 static GEN
FpX_FpXV_multirem_tree(GEN P,GEN xa,GEN T,GEN s,GEN p)1634 FpX_FpXV_multirem_tree(GEN P, GEN xa, GEN T, GEN s, GEN p)
1635 {
1636 pari_sp av = avma;
1637 long j, k, ls = lg(s);
1638 GEN Tp = FpX_FpXV_multirem_dbl_tree(P, T, p);
1639 GEN R = cgetg(lg(xa), t_VEC);
1640 GEN v = gel(Tp, 1);
1641 for (j=1, k=1; j<ls; k+=s[j++])
1642 {
1643 gel(R,k) = FpX_rem(gel(v, j), gel(xa,k), p);
1644 if (s[j] == 2)
1645 gel(R,k+1) = FpX_rem(gel(v, j), gel(xa,k+1), p);
1646 }
1647 return gerepileupto(av, R);
1648 }
1649
1650 GEN
FpX_FpXV_multirem(GEN P,GEN xa,GEN p)1651 FpX_FpXV_multirem(GEN P, GEN xa, GEN p)
1652 {
1653 pari_sp av = avma;
1654 GEN s = producttree_scheme(lg(xa)-1);
1655 GEN T = FpXV_producttree(xa, s, p);
1656 return gerepileupto(av, FpX_FpXV_multirem_tree(P, xa, T, s, p));
1657 }
1658
1659 /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
1660 static GEN
FpXV_chinese_tree(GEN A,GEN P,GEN T,GEN R,GEN s,GEN p)1661 FpXV_chinese_tree(GEN A, GEN P, GEN T, GEN R, GEN s, GEN p)
1662 {
1663 long m = lg(T)-1, ls = lg(s);
1664 long i,j,k;
1665 GEN Tp = cgetg(m+1, t_VEC);
1666 GEN M = gel(T, 1);
1667 GEN t = cgetg(lg(M), t_VEC);
1668 for (j=1, k=1; j<ls; k+=s[j++])
1669 if (s[j] == 2)
1670 {
1671 pari_sp av = avma;
1672 GEN a = FpX_mul(gel(A,k), gel(R,k), p), b = FpX_mul(gel(A,k+1), gel(R,k+1), p);
1673 GEN tj = FpX_rem(FpX_add(FpX_mul(gel(P,k), b, p),
1674 FpX_mul(gel(P,k+1), a, p), p), gel(M,j), p);
1675 gel(t, j) = gerepileupto(av, tj);
1676 }
1677 else
1678 gel(t, j) = FpX_rem(FpX_mul(gel(A,k), gel(R,k), p), gel(M, j), p);
1679 gel(Tp, 1) = t;
1680 for (i=2; i<=m; i++)
1681 {
1682 GEN u = gel(T, i-1), M = gel(T, i);
1683 GEN t = cgetg(lg(M), t_VEC);
1684 GEN v = gel(Tp, i-1);
1685 long n = lg(v)-1;
1686 for (j=1, k=1; k<n; j++, k+=2)
1687 {
1688 pari_sp av = avma;
1689 gel(t, j) = gerepileupto(av, FpX_rem(FpX_add(FpX_mul(gel(u, k), gel(v, k+1), p),
1690 FpX_mul(gel(u, k+1), gel(v, k), p), p), gel(M, j), p));
1691 }
1692 if (k==n) gel(t, j) = gel(v, k);
1693 gel(Tp, i) = t;
1694 }
1695 return gmael(Tp,m,1);
1696 }
1697
1698 static GEN
FpXV_sqr(GEN x,GEN p)1699 FpXV_sqr(GEN x, GEN p)
1700 { pari_APPLY_type(t_VEC, FpX_sqr(gel(x,i), p)) }
1701
1702 static GEN
FpXT_sqr(GEN x,GEN p)1703 FpXT_sqr(GEN x, GEN p)
1704 {
1705 if (typ(x) == t_POL)
1706 return FpX_sqr(x, p);
1707 pari_APPLY_type(t_VEC, FpXT_sqr(gel(x,i), p))
1708 }
1709
1710 static GEN
FpXV_invdivexact(GEN x,GEN y,GEN p)1711 FpXV_invdivexact(GEN x, GEN y, GEN p)
1712 { pari_APPLY_type(t_VEC, FpXQ_inv(FpX_div(gel(x,i), gel(y,i),p), gel(y,i),p)) }
1713
1714 static GEN
FpXV_chinesetree(GEN P,GEN T,GEN s,GEN p)1715 FpXV_chinesetree(GEN P, GEN T, GEN s, GEN p)
1716 {
1717 GEN T2 = FpXT_sqr(T, p), P2 = FpXV_sqr(P, p);
1718 GEN mod = gmael(T,lg(T)-1,1);
1719 return FpXV_invdivexact(FpX_FpXV_multirem_tree(mod, P2, T2, s, p), P, p);
1720 }
1721
1722 static GEN
gc_chinese(pari_sp av,GEN T,GEN a,GEN * pt_mod)1723 gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
1724 {
1725 if (!pt_mod)
1726 return gerepileupto(av, a);
1727 else
1728 {
1729 GEN mod = gmael(T, lg(T)-1, 1);
1730 gerepileall(av, 2, &a, &mod);
1731 *pt_mod = mod;
1732 return a;
1733 }
1734 }
1735
1736 GEN
FpXV_chinese(GEN A,GEN P,GEN p,GEN * pt_mod)1737 FpXV_chinese(GEN A, GEN P, GEN p, GEN *pt_mod)
1738 {
1739 pari_sp av = avma;
1740 GEN s = producttree_scheme(lg(P)-1);
1741 GEN T = FpXV_producttree(P, s, p);
1742 GEN R = FpXV_chinesetree(P, T, s, p);
1743 GEN a = FpXV_chinese_tree(A, P, T, R, s, p);
1744 return gc_chinese(av, T, a, pt_mod);
1745 }
1746
1747 /***********************************************************************/
1748 /** **/
1749 /** FpXQ **/
1750 /** **/
1751 /***********************************************************************/
1752
1753 /* FpXQ are elements of Fp[X]/(T), represented by FpX*/
1754
1755 GEN
FpXQ_red(GEN x,GEN T,GEN p)1756 FpXQ_red(GEN x, GEN T, GEN p)
1757 {
1758 GEN z = FpX_red(x,p);
1759 return FpX_rem(z, T,p);
1760 }
1761
1762 GEN
FpXQ_mul(GEN x,GEN y,GEN T,GEN p)1763 FpXQ_mul(GEN x,GEN y,GEN T,GEN p)
1764 {
1765 GEN z = FpX_mul(x,y,p);
1766 return FpX_rem(z, T, p);
1767 }
1768
1769 GEN
FpXQ_sqr(GEN x,GEN T,GEN p)1770 FpXQ_sqr(GEN x, GEN T, GEN p)
1771 {
1772 GEN z = FpX_sqr(x,p);
1773 return FpX_rem(z, T, p);
1774 }
1775
1776 /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
1777 * return lift(1 / (x mod (p,pol))) */
1778 GEN
FpXQ_invsafe(GEN x,GEN y,GEN p)1779 FpXQ_invsafe(GEN x, GEN y, GEN p)
1780 {
1781 GEN V, z = FpX_extgcd(get_FpX_mod(y), x, p, NULL, &V);
1782 if (degpol(z)) return NULL;
1783 z = Fp_invsafe(gel(z,2), p);
1784 if (!z) return NULL;
1785 return FpX_Fp_mul(V, z, p);
1786 }
1787
1788 GEN
FpXQ_inv(GEN x,GEN T,GEN p)1789 FpXQ_inv(GEN x,GEN T,GEN p)
1790 {
1791 pari_sp av = avma;
1792 GEN U = FpXQ_invsafe(x, T, p);
1793 if (!U) pari_err_INV("FpXQ_inv",x);
1794 return gerepileupto(av, U);
1795 }
1796
1797 GEN
FpXQ_div(GEN x,GEN y,GEN T,GEN p)1798 FpXQ_div(GEN x,GEN y,GEN T,GEN p)
1799 {
1800 pari_sp av = avma;
1801 return gerepileupto(av, FpXQ_mul(x,FpXQ_inv(y,T,p),T,p));
1802 }
1803
1804 static GEN
_FpXQ_add(void * data,GEN x,GEN y)1805 _FpXQ_add(void *data, GEN x, GEN y)
1806 {
1807 (void) data;
1808 return ZX_add(x, y);
1809 }
1810 static GEN
_FpXQ_sub(void * data,GEN x,GEN y)1811 _FpXQ_sub(void *data, GEN x, GEN y)
1812 {
1813 (void) data;
1814 return ZX_sub(x, y);
1815 }
1816 static GEN
_FpXQ_cmul(void * data,GEN P,long a,GEN x)1817 _FpXQ_cmul(void *data, GEN P, long a, GEN x)
1818 {
1819 (void) data;
1820 return ZX_Z_mul(x, gel(P,a+2));
1821 }
1822 static GEN
_FpXQ_sqr(void * data,GEN x)1823 _FpXQ_sqr(void *data, GEN x)
1824 {
1825 struct _FpXQ *D = (struct _FpXQ*)data;
1826 return FpXQ_sqr(x, D->T, D->p);
1827 }
1828 static GEN
_FpXQ_mul(void * data,GEN x,GEN y)1829 _FpXQ_mul(void *data, GEN x, GEN y)
1830 {
1831 struct _FpXQ *D = (struct _FpXQ*)data;
1832 return FpXQ_mul(x,y, D->T, D->p);
1833 }
1834 static GEN
_FpXQ_zero(void * data)1835 _FpXQ_zero(void *data)
1836 {
1837 struct _FpXQ *D = (struct _FpXQ*)data;
1838 return pol_0(get_FpX_var(D->T));
1839 }
1840 static GEN
_FpXQ_one(void * data)1841 _FpXQ_one(void *data)
1842 {
1843 struct _FpXQ *D = (struct _FpXQ*)data;
1844 return pol_1(get_FpX_var(D->T));
1845 }
1846 static GEN
_FpXQ_red(void * data,GEN x)1847 _FpXQ_red(void *data, GEN x)
1848 {
1849 struct _FpXQ *D = (struct _FpXQ*)data;
1850 return FpX_red(x,D->p);
1851 }
1852
1853 static struct bb_algebra FpXQ_algebra = { _FpXQ_red, _FpXQ_add, _FpXQ_sub,
1854 _FpXQ_mul, _FpXQ_sqr, _FpXQ_one, _FpXQ_zero };
1855
1856 const struct bb_algebra *
get_FpXQ_algebra(void ** E,GEN T,GEN p)1857 get_FpXQ_algebra(void **E, GEN T, GEN p)
1858 {
1859 GEN z = new_chunk(sizeof(struct _FpXQ));
1860 struct _FpXQ *e = (struct _FpXQ *) z;
1861 e->T = FpX_get_red(T, p);
1862 e->p = p; *E = (void*)e;
1863 return &FpXQ_algebra;
1864 }
1865
1866 static struct bb_algebra FpX_algebra = { _FpXQ_red, _FpXQ_add, _FpXQ_sub,
1867 _FpX_mul, _FpX_sqr, _FpXQ_one, _FpXQ_zero };
1868
1869 const struct bb_algebra *
get_FpX_algebra(void ** E,GEN p,long v)1870 get_FpX_algebra(void **E, GEN p, long v)
1871 {
1872 GEN z = new_chunk(sizeof(struct _FpXQ));
1873 struct _FpXQ *e = (struct _FpXQ *) z;
1874 e->T = pol_x(v);
1875 e->p = p; *E = (void*)e;
1876 return &FpX_algebra;
1877 }
1878
1879 /* x,pol in Z[X], p in Z, n in Z, compute lift(x^n mod (p, pol)) */
1880 GEN
FpXQ_pow(GEN x,GEN n,GEN T,GEN p)1881 FpXQ_pow(GEN x, GEN n, GEN T, GEN p)
1882 {
1883 struct _FpXQ D;
1884 pari_sp av;
1885 long s = signe(n);
1886 GEN y;
1887 if (!s) return pol_1(varn(x));
1888 if (is_pm1(n)) /* +/- 1 */
1889 return (s < 0)? FpXQ_inv(x,T,p): FpXQ_red(x,T,p);
1890 av = avma;
1891 if (!is_bigint(p))
1892 {
1893 ulong pp = to_Flxq(&x, &T, p);
1894 y = Flxq_pow(x, n, T, pp);
1895 return Flx_to_ZX_inplace(gerepileuptoleaf(av, y));
1896 }
1897 if (s < 0) x = FpXQ_inv(x,T,p);
1898 D.p = p; D.T = FpX_get_red(T,p);
1899 y = gen_pow_i(x, n, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul);
1900 return gerepilecopy(av, y);
1901 }
1902
1903 GEN /*Assume n is very small*/
FpXQ_powu(GEN x,ulong n,GEN T,GEN p)1904 FpXQ_powu(GEN x, ulong n, GEN T, GEN p)
1905 {
1906 struct _FpXQ D;
1907 pari_sp av;
1908 GEN y;
1909 if (!n) return pol_1(varn(x));
1910 if (n==1) return FpXQ_red(x,T,p);
1911 av = avma;
1912 if (!is_bigint(p))
1913 {
1914 ulong pp = to_Flxq(&x, &T, p);
1915 y = Flxq_powu(x, n, T, pp);
1916 return Flx_to_ZX_inplace(gerepileuptoleaf(av, y));
1917 }
1918 D.T = FpX_get_red(T, p); D.p = p;
1919 y = gen_powu_i(x, n, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul);
1920 return gerepilecopy(av, y);
1921 }
1922
1923 /* generates the list of powers of x of degree 0,1,2,...,l*/
1924 GEN
FpXQ_powers(GEN x,long l,GEN T,GEN p)1925 FpXQ_powers(GEN x, long l, GEN T, GEN p)
1926 {
1927 struct _FpXQ D;
1928 int use_sqr;
1929 if (l>2 && lgefint(p) == 3) {
1930 pari_sp av = avma;
1931 ulong pp = to_Flxq(&x, &T, p);
1932 GEN z = FlxV_to_ZXV(Flxq_powers(x, l, T, pp));
1933 return gerepileupto(av, z);
1934 }
1935 use_sqr = 2*degpol(x)>=get_FpX_degree(T);
1936 D.T = FpX_get_red(T,p); D.p = p;
1937 return gen_powers(x, l, use_sqr, (void*)&D, &_FpXQ_sqr, &_FpXQ_mul,&_FpXQ_one);
1938 }
1939
1940 GEN
FpXQ_matrix_pow(GEN y,long n,long m,GEN P,GEN l)1941 FpXQ_matrix_pow(GEN y, long n, long m, GEN P, GEN l)
1942 {
1943 return RgXV_to_RgM(FpXQ_powers(y,m-1,P,l),n);
1944 }
1945
1946 GEN
FpX_Frobenius(GEN T,GEN p)1947 FpX_Frobenius(GEN T, GEN p)
1948 {
1949 return FpXQ_pow(pol_x(get_FpX_var(T)), p, T, p);
1950 }
1951
1952 GEN
FpX_matFrobenius(GEN T,GEN p)1953 FpX_matFrobenius(GEN T, GEN p)
1954 {
1955 long n = get_FpX_degree(T);
1956 return FpXQ_matrix_pow(FpX_Frobenius(T, p), n, n, T, p);
1957 }
1958
1959 GEN
FpX_FpXQV_eval(GEN Q,GEN x,GEN T,GEN p)1960 FpX_FpXQV_eval(GEN Q, GEN x, GEN T, GEN p)
1961 {
1962 struct _FpXQ D;
1963 D.T = FpX_get_red(T,p); D.p = p;
1964 return gen_bkeval_powers(Q,degpol(Q),x,(void*)&D,&FpXQ_algebra,_FpXQ_cmul);
1965 }
1966
1967 GEN
FpX_FpXQ_eval(GEN Q,GEN x,GEN T,GEN p)1968 FpX_FpXQ_eval(GEN Q, GEN x, GEN T, GEN p)
1969 {
1970 struct _FpXQ D;
1971 int use_sqr;
1972 if (lgefint(p) == 3)
1973 {
1974 pari_sp av = avma;
1975 ulong pp = to_Flxq(&x, &T, p);
1976 GEN z = Flx_Flxq_eval(ZX_to_Flx(Q, pp), x, T, pp);
1977 return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
1978 }
1979 use_sqr = 2*degpol(x) >= get_FpX_degree(T);
1980 D.T = FpX_get_red(T,p); D.p = p;
1981 return gen_bkeval(Q,degpol(Q),x,use_sqr,(void*)&D,&FpXQ_algebra,_FpXQ_cmul);
1982 }
1983
1984 GEN
FpXC_FpXQV_eval(GEN x,GEN v,GEN T,GEN p)1985 FpXC_FpXQV_eval(GEN x, GEN v, GEN T, GEN p)
1986 { pari_APPLY_type(t_COL, FpX_FpXQV_eval(gel(x,i), v, T, p)) }
1987
1988 GEN
FpXM_FpXQV_eval(GEN x,GEN v,GEN T,GEN p)1989 FpXM_FpXQV_eval(GEN x, GEN v, GEN T, GEN p)
1990 { pari_APPLY_same(FpXC_FpXQV_eval(gel(x,i), v, T, p)) }
1991
1992 GEN
FpXC_FpXQ_eval(GEN x,GEN F,GEN T,GEN p)1993 FpXC_FpXQ_eval(GEN x, GEN F, GEN T, GEN p)
1994 {
1995 long d = brent_kung_optpow(RgXV_maxdegree(x), lg(x)-1, 1);
1996 GEN Fp = FpXQ_powers(F, d, T, p);
1997 return FpXC_FpXQV_eval(x, Fp, T, p);
1998 }
1999
2000 GEN
FpXQ_autpowers(GEN aut,long f,GEN T,GEN p)2001 FpXQ_autpowers(GEN aut, long f, GEN T, GEN p)
2002 {
2003 pari_sp av = avma;
2004 long n = get_FpX_degree(T);
2005 long i, nautpow = brent_kung_optpow(n-1,f-2,1);
2006 long v = get_FpX_var(T);
2007 GEN autpow, V;
2008 T = FpX_get_red(T, p);
2009 autpow = FpXQ_powers(aut, nautpow,T,p);
2010 V = cgetg(f + 2, t_VEC);
2011 gel(V,1) = pol_x(v); if (f==0) return gerepileupto(av, V);
2012 gel(V,2) = gcopy(aut);
2013 for (i = 3; i <= f+1; i++)
2014 gel(V,i) = FpX_FpXQV_eval(gel(V,i-1),autpow,T,p);
2015 return gerepileupto(av, V);
2016 }
2017
2018 static GEN
FpXQ_autpow_sqr(void * E,GEN x)2019 FpXQ_autpow_sqr(void *E, GEN x)
2020 {
2021 struct _FpXQ *D = (struct _FpXQ*)E;
2022 return FpX_FpXQ_eval(x, x, D->T, D->p);
2023 }
2024
2025 static GEN
FpXQ_autpow_msqr(void * E,GEN x)2026 FpXQ_autpow_msqr(void *E, GEN x)
2027 {
2028 struct _FpXQ *D = (struct _FpXQ*)E;
2029 return FpX_FpXQV_eval(FpXQ_autpow_sqr(E, x), D->aut, D->T, D->p);
2030 }
2031
2032 GEN
FpXQ_autpow(GEN x,ulong n,GEN T,GEN p)2033 FpXQ_autpow(GEN x, ulong n, GEN T, GEN p)
2034 {
2035 pari_sp av = avma;
2036 struct _FpXQ D;
2037 long d;
2038 if (n==0) return FpX_rem(pol_x(varn(x)), T, p);
2039 if (n==1) return FpX_rem(x, T, p);
2040 D.T = FpX_get_red(T, p); D.p = p;
2041 d = brent_kung_optpow(degpol(T), hammingl(n)-1, 1);
2042 D.aut = FpXQ_powers(x, d, T, p);
2043 x = gen_powu_fold(x,n,(void*)&D,FpXQ_autpow_sqr,FpXQ_autpow_msqr);
2044 return gerepilecopy(av, x);
2045 }
2046
2047 static GEN
FpXQ_auttrace_mul(void * E,GEN x,GEN y)2048 FpXQ_auttrace_mul(void *E, GEN x, GEN y)
2049 {
2050 struct _FpXQ *D = (struct _FpXQ*)E;
2051 GEN T = D->T, p = D->p;
2052 GEN phi1 = gel(x,1), a1 = gel(x,2);
2053 GEN phi2 = gel(y,1), a2 = gel(y,2);
2054 ulong d = brent_kung_optpow(maxss(degpol(phi2),degpol(a2)),2,1);
2055 GEN V1 = FpXQ_powers(phi1, d, T, p);
2056 GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
2057 GEN aphi = FpX_FpXQV_eval(a2, V1, T, p);
2058 GEN a3 = FpX_add(a1, aphi, p);
2059 return mkvec2(phi3, a3);
2060 }
2061
2062 static GEN
FpXQ_auttrace_sqr(void * E,GEN x)2063 FpXQ_auttrace_sqr(void *E, GEN x)
2064 { return FpXQ_auttrace_mul(E, x, x); }
2065
2066 GEN
FpXQ_auttrace(GEN x,ulong n,GEN T,GEN p)2067 FpXQ_auttrace(GEN x, ulong n, GEN T, GEN p)
2068 {
2069 pari_sp av = avma;
2070 struct _FpXQ D;
2071 D.T = FpX_get_red(T, p); D.p = p;
2072 x = gen_powu_i(x,n,(void*)&D,FpXQ_auttrace_sqr,FpXQ_auttrace_mul);
2073 return gerepilecopy(av, x);
2074 }
2075
2076 static GEN
FpXQ_autsum_mul(void * E,GEN x,GEN y)2077 FpXQ_autsum_mul(void *E, GEN x, GEN y)
2078 {
2079 struct _FpXQ *D = (struct _FpXQ*)E;
2080 GEN T = D->T, p = D->p;
2081 GEN phi1 = gel(x,1), a1 = gel(x,2);
2082 GEN phi2 = gel(y,1), a2 = gel(y,2);
2083 ulong d = brent_kung_optpow(maxss(degpol(phi2),degpol(a2)),2,1);
2084 GEN V1 = FpXQ_powers(phi1, d, T, p);
2085 GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
2086 GEN aphi = FpX_FpXQV_eval(a2, V1, T, p);
2087 GEN a3 = FpXQ_mul(a1, aphi, T, p);
2088 return mkvec2(phi3, a3);
2089 }
2090 static GEN
FpXQ_autsum_sqr(void * E,GEN x)2091 FpXQ_autsum_sqr(void *E, GEN x)
2092 { return FpXQ_autsum_mul(E, x, x); }
2093
2094 GEN
FpXQ_autsum(GEN x,ulong n,GEN T,GEN p)2095 FpXQ_autsum(GEN x, ulong n, GEN T, GEN p)
2096 {
2097 pari_sp av = avma;
2098 struct _FpXQ D;
2099 D.T = FpX_get_red(T, p); D.p = p;
2100 x = gen_powu_i(x,n,(void*)&D,FpXQ_autsum_sqr,FpXQ_autsum_mul);
2101 return gerepilecopy(av, x);
2102 }
2103
2104 static GEN
FpXQM_autsum_mul(void * E,GEN x,GEN y)2105 FpXQM_autsum_mul(void *E, GEN x, GEN y)
2106 {
2107 struct _FpXQ *D = (struct _FpXQ*)E;
2108 GEN T = D->T, p = D->p;
2109 GEN phi1 = gel(x,1), a1 = gel(x,2);
2110 GEN phi2 = gel(y,1), a2 = gel(y,2);
2111 long g = lg(a2)-1, dT = get_FpX_degree(T);
2112 ulong d = brent_kung_optpow(dT-1, g*g+1, 1);
2113 GEN V1 = FpXQ_powers(phi1, d, T, p);
2114 GEN phi3 = FpX_FpXQV_eval(phi2, V1, T, p);
2115 GEN aphi = FpXM_FpXQV_eval(a2, V1, T, p);
2116 GEN a3 = FqM_mul(a1, aphi, T, p);
2117 return mkvec2(phi3, a3);
2118 }
2119 static GEN
FpXQM_autsum_sqr(void * E,GEN x)2120 FpXQM_autsum_sqr(void *E, GEN x)
2121 { return FpXQM_autsum_mul(E, x, x); }
2122
2123 GEN
FpXQM_autsum(GEN x,ulong n,GEN T,GEN p)2124 FpXQM_autsum(GEN x, ulong n, GEN T, GEN p)
2125 {
2126 pari_sp av = avma;
2127 struct _FpXQ D;
2128 D.T = FpX_get_red(T, p); D.p = p;
2129 x = gen_powu_i(x, n, (void*)&D, FpXQM_autsum_sqr, FpXQM_autsum_mul);
2130 return gerepilecopy(av, x);
2131 }
2132
2133 static long
bounded_order(GEN p,GEN b,long k)2134 bounded_order(GEN p, GEN b, long k)
2135 {
2136 long i;
2137 GEN a=modii(p,b);
2138 for(i=1;i<k;i++)
2139 {
2140 if (equali1(a))
2141 return i;
2142 a = Fp_mul(a,p,b);
2143 }
2144 return 0;
2145 }
2146
2147 /*
2148 n = (p^d-a)\b
2149 b = bb*p^vb
2150 p^k = 1 [bb]
2151 d = m*k+r+vb
2152 u = (p^k-1)/bb;
2153 v = (p^(r+vb)-a)/b;
2154 w = (p^(m*k)-1)/(p^k-1)
2155 n = p^r*w*u+v
2156 w*u = p^vb*(p^(m*k)-1)/b
2157 n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b
2158 */
2159
2160 static GEN
FpXQ_pow_Frobenius(GEN x,GEN n,GEN aut,GEN T,GEN p)2161 FpXQ_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, GEN p)
2162 {
2163 pari_sp av=avma;
2164 long d = get_FpX_degree(T);
2165 GEN an = absi_shallow(n), z, q;
2166 if (cmpii(an,p)<0 || cmpis(an,d)<=0) return FpXQ_pow(x, n, T, p);
2167 q = powiu(p, d);
2168 if (dvdii(q, n))
2169 {
2170 long vn = logint(an,p);
2171 GEN autvn = vn==1 ? aut: FpXQ_autpow(aut,vn,T,p);
2172 z = FpX_FpXQ_eval(x,autvn,T,p);
2173 } else
2174 {
2175 GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
2176 GEN bb, u, v, autk;
2177 long vb = Z_pvalrem(b,p,&bb);
2178 long m, r, k = is_pm1(bb) ? 1 : bounded_order(p,bb,d);
2179 if (!k || d-vb<k) return FpXQ_pow(x,n, T, p);
2180 m = (d-vb)/k; r = (d-vb)%k;
2181 u = diviiexact(subiu(powiu(p,k),1),bb);
2182 v = diviiexact(subii(powiu(p,r+vb),a),b);
2183 autk = k==1 ? aut: FpXQ_autpow(aut,k,T,p);
2184 if (r)
2185 {
2186 GEN autr = r==1 ? aut: FpXQ_autpow(aut,r,T,p);
2187 z = FpX_FpXQ_eval(x,autr,T,p);
2188 } else z = x;
2189 if (m > 1) z = gel(FpXQ_autsum(mkvec2(autk, z), m, T, p), 2);
2190 if (!is_pm1(u)) z = FpXQ_pow(z, u, T, p);
2191 if (signe(v)) z = FpXQ_mul(z, FpXQ_pow(x, v, T, p), T, p);
2192 }
2193 return gerepileupto(av,signe(n)>0 ? z : FpXQ_inv(z,T,p));
2194 }
2195
2196 /* assume T irreducible mod p */
2197 int
FpXQ_issquare(GEN x,GEN T,GEN p)2198 FpXQ_issquare(GEN x, GEN T, GEN p)
2199 {
2200 pari_sp av;
2201 if (lg(x) == 2 || absequalui(2, p)) return 1;
2202 if (lg(x) == 3) return Fq_issquare(gel(x,2), T, p);
2203 av = avma; /* Ng = g^((q-1)/(p-1)) */
2204 return gc_bool(av, kronecker(FpXQ_norm(x,T,p), p) != -1);
2205 }
2206 int
Fp_issquare(GEN x,GEN p)2207 Fp_issquare(GEN x, GEN p)
2208 { return absequalui(2, p) || kronecker(x, p) != -1; }
2209 /* assume T irreducible mod p */
2210 int
Fq_issquare(GEN x,GEN T,GEN p)2211 Fq_issquare(GEN x, GEN T, GEN p)
2212 {
2213 if (typ(x) != t_INT) return FpXQ_issquare(x, T, p);
2214 return (T && ! odd(get_FpX_degree(T))) || Fp_issquare(x, p);
2215 }
2216
2217 long
Fq_ispower(GEN x,GEN K,GEN T,GEN p)2218 Fq_ispower(GEN x, GEN K, GEN T, GEN p)
2219 {
2220 pari_sp av = avma;
2221 long d;
2222 GEN Q;
2223 if (equaliu(K,2)) return Fq_issquare(x, T, p);
2224 if (!T) return Fp_ispower(x, K, p);
2225 d = get_FpX_degree(T);
2226 if (typ(x) == t_INT && !umodui(d, K)) return 1;
2227 Q = subiu(powiu(p,d), 1);
2228 Q = diviiexact(Q, gcdii(Q, K));
2229 d = gequal1(Fq_pow(x, Q, T,p));
2230 return gc_long(av, d);
2231 }
2232
2233 /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
2234 GEN
Fp_FpXQ_log(GEN a,GEN g,GEN o,GEN T,GEN p)2235 Fp_FpXQ_log(GEN a, GEN g, GEN o, GEN T, GEN p)
2236 {
2237 pari_sp av = avma;
2238 GEN q,n_q,ord,ordp, op;
2239
2240 if (equali1(a)) return gen_0;
2241 /* p > 2 */
2242
2243 ordp = subiu(p, 1); /* even */
2244 ord = get_arith_Z(o);
2245 if (!ord) ord = T? subiu(powiu(p, get_FpX_degree(T)), 1): ordp;
2246 if (equalii(a, ordp)) /* -1 */
2247 return gerepileuptoint(av, shifti(ord,-1));
2248 ordp = gcdii(ordp,ord);
2249 op = typ(o)==t_MAT ? famat_Z_gcd(o,ordp) : ordp;
2250
2251 q = NULL;
2252 if (T)
2253 { /* we want < g > = Fp^* */
2254 if (!equalii(ord,ordp)) {
2255 q = diviiexact(ord,ordp);
2256 g = FpXQ_pow(g,q,T,p);
2257 }
2258 g = constant_coeff(g);
2259 }
2260 n_q = Fp_log(a,g,op,p);
2261 if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
2262 if (q) n_q = mulii(q, n_q);
2263 return gerepileuptoint(av, n_q);
2264 }
2265
2266 static GEN
_FpXQ_pow(void * data,GEN x,GEN n)2267 _FpXQ_pow(void *data, GEN x, GEN n)
2268 {
2269 struct _FpXQ *D = (struct _FpXQ*)data;
2270 return FpXQ_pow_Frobenius(x,n, D->aut, D->T, D->p);
2271 }
2272
2273 static GEN
_FpXQ_rand(void * data)2274 _FpXQ_rand(void *data)
2275 {
2276 pari_sp av=avma;
2277 struct _FpXQ *D = (struct _FpXQ*)data;
2278 GEN z;
2279 do
2280 {
2281 set_avma(av);
2282 z=random_FpX(get_FpX_degree(D->T),get_FpX_var(D->T),D->p);
2283 } while (!signe(z));
2284 return z;
2285 }
2286
2287 static GEN
_FpXQ_easylog(void * E,GEN a,GEN g,GEN ord)2288 _FpXQ_easylog(void *E, GEN a, GEN g, GEN ord)
2289 {
2290 struct _FpXQ *s=(struct _FpXQ*) E;
2291 if (degpol(a)) return NULL;
2292 return Fp_FpXQ_log(constant_coeff(a),g,ord,s->T,s->p);
2293 }
2294
2295 static const struct bb_group FpXQ_star={_FpXQ_mul,_FpXQ_pow,_FpXQ_rand,hash_GEN,ZX_equal,ZX_equal1,_FpXQ_easylog};
2296
2297 const struct bb_group *
get_FpXQ_star(void ** E,GEN T,GEN p)2298 get_FpXQ_star(void **E, GEN T, GEN p)
2299 {
2300 struct _FpXQ *e = (struct _FpXQ *) stack_malloc(sizeof(struct _FpXQ));
2301 e->T = T; e->p = p; e->aut = FpX_Frobenius(T, p);
2302 *E = (void*)e; return &FpXQ_star;
2303 }
2304
2305 GEN
FpXQ_order(GEN a,GEN ord,GEN T,GEN p)2306 FpXQ_order(GEN a, GEN ord, GEN T, GEN p)
2307 {
2308 if (lgefint(p)==3)
2309 {
2310 pari_sp av=avma;
2311 ulong pp = to_Flxq(&a, &T, p);
2312 GEN z = Flxq_order(a, ord, T, pp);
2313 return gerepileuptoint(av,z);
2314 }
2315 else
2316 {
2317 void *E;
2318 const struct bb_group *S = get_FpXQ_star(&E,T,p);
2319 return gen_order(a,ord,E,S);
2320 }
2321 }
2322
2323 GEN
FpXQ_log(GEN a,GEN g,GEN ord,GEN T,GEN p)2324 FpXQ_log(GEN a, GEN g, GEN ord, GEN T, GEN p)
2325 {
2326 pari_sp av=avma;
2327 if (lgefint(p)==3)
2328 {
2329 if (uel(p,2) == 2)
2330 {
2331 GEN z = F2xq_log(ZX_to_F2x(a), ZX_to_F2x(g), ord,
2332 ZX_to_F2x(get_FpX_mod(T)));
2333 return gerepileuptoleaf(av, z);
2334 }
2335 else
2336 {
2337 ulong pp = to_Flxq(&a, &T, p);
2338 GEN z = Flxq_log(a, ZX_to_Flx(g, pp), ord, T, pp);
2339 return gerepileuptoleaf(av, z);
2340 }
2341 }
2342 else
2343 {
2344 void *E;
2345 const struct bb_group *S = get_FpXQ_star(&E,T,p);
2346 GEN z = gen_PH_log(a,g,ord,E,S);
2347 return gerepileuptoleaf(av, z);
2348 }
2349 }
2350
2351 GEN
Fq_log(GEN a,GEN g,GEN ord,GEN T,GEN p)2352 Fq_log(GEN a, GEN g, GEN ord, GEN T, GEN p)
2353 {
2354 if (!T) return Fp_log(a,g,ord,p);
2355 if (typ(g) == t_INT)
2356 {
2357 if (typ(a) == t_POL)
2358 {
2359 if (degpol(a)) return cgetg(1,t_VEC);
2360 a = gel(a,2);
2361 }
2362 return Fp_log(a,g,ord,p);
2363 }
2364 return typ(a) == t_INT? Fp_FpXQ_log(a,g,ord,T,p): FpXQ_log(a,g,ord,T,p);
2365 }
2366
2367 GEN
FpXQ_sqrtn(GEN a,GEN n,GEN T,GEN p,GEN * zeta)2368 FpXQ_sqrtn(GEN a, GEN n, GEN T, GEN p, GEN *zeta)
2369 {
2370 pari_sp av = avma;
2371 GEN z;
2372 if (!signe(a))
2373 {
2374 long v=varn(a);
2375 if (signe(n) < 0) pari_err_INV("FpXQ_sqrtn",a);
2376 if (zeta) *zeta=pol_1(v);
2377 return pol_0(v);
2378 }
2379 if (lgefint(p)==3)
2380 {
2381 if (uel(p,2) == 2)
2382 {
2383 z = F2xq_sqrtn(ZX_to_F2x(a), n, ZX_to_F2x(get_FpX_mod(T)), zeta);
2384 if (!z) return NULL;
2385 z = F2x_to_ZX(z);
2386 if (!zeta) return gerepileuptoleaf(av, z);
2387 *zeta=F2x_to_ZX(*zeta);
2388 } else
2389 {
2390 ulong pp = to_Flxq(&a, &T, p);
2391 z = Flxq_sqrtn(a, n, T, pp, zeta);
2392 if (!z) return NULL;
2393 if (!zeta) return Flx_to_ZX_inplace(gerepileuptoleaf(av, z));
2394 z = Flx_to_ZX(z);
2395 *zeta=Flx_to_ZX(*zeta);
2396 }
2397 }
2398 else
2399 {
2400 void *E;
2401 const struct bb_group *S = get_FpXQ_star(&E,T,p);
2402 GEN o = subiu(powiu(p,get_FpX_degree(T)),1);
2403 z = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
2404 if (!z) return NULL;
2405 if (!zeta) return gerepileupto(av, z);
2406 }
2407 gerepileall(av, 2, &z,zeta);
2408 return z;
2409 }
2410
2411 GEN
FpXQ_sqrt(GEN a,GEN T,GEN p)2412 FpXQ_sqrt(GEN a, GEN T, GEN p)
2413 {
2414 return FpXQ_sqrtn(a, gen_2, T, p, NULL);
2415 }
2416
2417 GEN
FpXQ_norm(GEN x,GEN TB,GEN p)2418 FpXQ_norm(GEN x, GEN TB, GEN p)
2419 {
2420 pari_sp av = avma;
2421 GEN T = get_FpX_mod(TB);
2422 GEN y = FpX_resultant(T, x, p);
2423 GEN L = leading_coeff(T);
2424 if (gequal1(L) || signe(x)==0) return y;
2425 return gerepileupto(av, Fp_div(y, Fp_pows(L, degpol(x), p), p));
2426 }
2427
2428 GEN
FpXQ_trace(GEN x,GEN TB,GEN p)2429 FpXQ_trace(GEN x, GEN TB, GEN p)
2430 {
2431 pari_sp av = avma;
2432 GEN T = get_FpX_mod(TB);
2433 GEN dT = FpX_deriv(T,p);
2434 long n = degpol(dT);
2435 GEN z = FpXQ_mul(x, dT, TB, p);
2436 if (degpol(z)<n) return gc_const(av, gen_0);
2437 return gerepileuptoint(av, Fp_div(gel(z,2+n), gel(T,3+n),p));
2438 }
2439
2440 GEN
FpXQ_charpoly(GEN x,GEN T,GEN p)2441 FpXQ_charpoly(GEN x, GEN T, GEN p)
2442 {
2443 pari_sp ltop=avma;
2444 long vT, v = fetch_var();
2445 GEN R;
2446 T = leafcopy(get_FpX_mod(T));
2447 vT = varn(T); setvarn(T, v);
2448 x = leafcopy(x); setvarn(x, v);
2449 R = FpX_FpXY_resultant(T, deg1pol_shallow(gen_1,FpX_neg(x,p),vT),p);
2450 (void)delete_var(); return gerepileupto(ltop,R);
2451 }
2452
2453 /* Computing minimal polynomial : */
2454 /* cf Shoup 'Efficient Computation of Minimal Polynomials */
2455 /* in Algebraic Extensions of Finite Fields' */
2456
2457 /* Let v a linear form, return the linear form z->v(tau*z)
2458 that is, v*(M_tau) */
2459
2460 static GEN
FpXQ_transmul_init(GEN tau,GEN T,GEN p)2461 FpXQ_transmul_init(GEN tau, GEN T, GEN p)
2462 {
2463 GEN bht;
2464 GEN h, Tp = get_FpX_red(T, &h);
2465 long n = degpol(Tp), vT = varn(Tp);
2466 GEN ft = FpX_recipspec(Tp+2, n+1, n+1);
2467 GEN bt = FpX_recipspec(tau+2, lgpol(tau), n);
2468 setvarn(ft, vT); setvarn(bt, vT);
2469 if (h)
2470 bht = FpXn_mul(bt, h, n-1, p);
2471 else
2472 {
2473 GEN bh = FpX_div(FpX_shift(tau, n-1), T, p);
2474 bht = FpX_recipspec(bh+2, lgpol(bh), n-1);
2475 setvarn(bht, vT);
2476 }
2477 return mkvec3(bt, bht, ft);
2478 }
2479
2480 static GEN
FpXQ_transmul(GEN tau,GEN a,long n,GEN p)2481 FpXQ_transmul(GEN tau, GEN a, long n, GEN p)
2482 {
2483 pari_sp ltop = avma;
2484 GEN t1, t2, t3, vec;
2485 GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
2486 if (signe(a)==0) return pol_0(varn(a));
2487 t2 = FpX_shift(FpX_mul(bt, a, p),1-n);
2488 if (signe(bht)==0) return gerepilecopy(ltop, t2);
2489 t1 = FpX_shift(FpX_mul(ft, a, p),-n);
2490 t3 = FpXn_mul(t1, bht, n-1, p);
2491 vec = FpX_sub(t2, FpX_shift(t3, 1), p);
2492 return gerepileupto(ltop, vec);
2493 }
2494
2495 GEN
FpXQ_minpoly(GEN x,GEN T,GEN p)2496 FpXQ_minpoly(GEN x, GEN T, GEN p)
2497 {
2498 pari_sp ltop = avma;
2499 long vT, n;
2500 GEN v_x, g, tau;
2501 if (lgefint(p)==3)
2502 {
2503 ulong pp = to_Flxq(&x, &T, p);
2504 GEN g = Flxq_minpoly(x, T, pp);
2505 return gerepileupto(ltop, Flx_to_ZX(g));
2506 }
2507 vT = get_FpX_var(T);
2508 n = get_FpX_degree(T);
2509 g = pol_1(vT);
2510 tau = pol_1(vT);
2511 T = FpX_get_red(T, p);
2512 x = FpXQ_red(x, T, p);
2513 v_x = FpXQ_powers(x, usqrt(2*n), T, p);
2514 while(signe(tau) != 0)
2515 {
2516 long i, j, m, k1;
2517 GEN M, v, tr;
2518 GEN g_prime, c;
2519 if (degpol(g) == n) { tau = pol_1(vT); g = pol_1(vT); }
2520 v = random_FpX(n, vT, p);
2521 tr = FpXQ_transmul_init(tau, T, p);
2522 v = FpXQ_transmul(tr, v, n, p);
2523 m = 2*(n-degpol(g));
2524 k1 = usqrt(m);
2525 tr = FpXQ_transmul_init(gel(v_x,k1+1), T, p);
2526 c = cgetg(m+2,t_POL);
2527 c[1] = evalsigne(1)|evalvarn(vT);
2528 for (i=0; i<m; i+=k1)
2529 {
2530 long mj = minss(m-i, k1);
2531 for (j=0; j<mj; j++)
2532 gel(c,m+1-(i+j)) = FpX_dotproduct(v, gel(v_x,j+1), p);
2533 v = FpXQ_transmul(tr, v, n, p);
2534 }
2535 c = FpX_renormalize(c, m+2);
2536 /* now c contains <v,x^i> , i = 0..m-1 */
2537 M = FpX_halfgcd(pol_xn(m, vT), c, p);
2538 g_prime = gmael(M, 2, 2);
2539 if (degpol(g_prime) < 1) continue;
2540 g = FpX_mul(g, g_prime, p);
2541 tau = FpXQ_mul(tau, FpX_FpXQV_eval(g_prime, v_x, T, p), T, p);
2542 }
2543 g = FpX_normalize(g,p);
2544 return gerepilecopy(ltop,g);
2545 }
2546
2547 GEN
FpXQ_conjvec(GEN x,GEN T,GEN p)2548 FpXQ_conjvec(GEN x, GEN T, GEN p)
2549 {
2550 pari_sp av=avma;
2551 long i;
2552 long n = get_FpX_degree(T), v = varn(x);
2553 GEN M = FpX_matFrobenius(T, p);
2554 GEN z = cgetg(n+1,t_COL);
2555 gel(z,1) = RgX_to_RgC(x,n);
2556 for (i=2; i<=n; i++) gel(z,i) = FpM_FpC_mul(M,gel(z,i-1),p);
2557 gel(z,1) = x;
2558 for (i=2; i<=n; i++) gel(z,i) = RgV_to_RgX(gel(z,i),v);
2559 return gerepilecopy(av,z);
2560 }
2561
2562 /* p prime, p_1 = p-1, q = p^deg T, Lp = cofactors of some prime divisors
2563 * l_p of p-1, Lq = cofactors of some prime divisors l_q of q-1, return a
2564 * g in Fq such that
2565 * - Ng generates all l_p-Sylows of Fp^*
2566 * - g generates all l_q-Sylows of Fq^* */
2567 static GEN
gener_FpXQ_i(GEN T,GEN p,GEN p_1,GEN Lp,GEN Lq)2568 gener_FpXQ_i(GEN T, GEN p, GEN p_1, GEN Lp, GEN Lq)
2569 {
2570 pari_sp av;
2571 long vT = varn(T), f = degpol(T), l = lg(Lq);
2572 GEN F = FpX_Frobenius(T, p);
2573 int p_is_2 = is_pm1(p_1);
2574 for (av = avma;; set_avma(av))
2575 {
2576 GEN t, g = random_FpX(f, vT, p);
2577 long i;
2578 if (degpol(g) < 1) continue;
2579 if (p_is_2)
2580 t = g;
2581 else
2582 {
2583 t = FpX_resultant(T, g, p); /* Ng = g^((q-1)/(p-1)), assuming T monic */
2584 if (kronecker(t, p) == 1) continue;
2585 if (lg(Lp) > 1 && !is_gener_Fp(t, p, p_1, Lp)) continue;
2586 t = FpXQ_pow(g, shifti(p_1,-1), T, p);
2587 }
2588 for (i = 1; i < l; i++)
2589 {
2590 GEN a = FpXQ_pow_Frobenius(t, gel(Lq,i), F, T, p);
2591 if (!degpol(a) && equalii(gel(a,2), p_1)) break;
2592 }
2593 if (i == l) return g;
2594 }
2595 }
2596
2597 GEN
gener_FpXQ(GEN T,GEN p,GEN * po)2598 gener_FpXQ(GEN T, GEN p, GEN *po)
2599 {
2600 long i, j, f = get_FpX_degree(T);
2601 GEN g, Lp, Lq, p_1, q_1, N, o;
2602 pari_sp av = avma;
2603
2604 p_1 = subiu(p,1);
2605 if (f == 1) {
2606 GEN Lp, fa;
2607 o = p_1;
2608 fa = Z_factor(o);
2609 Lp = gel(fa,1);
2610 Lp = vecslice(Lp, 2, lg(Lp)-1); /* remove 2 for efficiency */
2611
2612 g = cgetg(3, t_POL);
2613 g[1] = evalsigne(1) | evalvarn(get_FpX_var(T));
2614 gel(g,2) = pgener_Fp_local(p, Lp);
2615 if (po) *po = mkvec2(o, fa);
2616 return g;
2617 }
2618 if (lgefint(p) == 3)
2619 {
2620 ulong pp = to_Flxq(NULL, &T, p);
2621 g = gener_Flxq(T, pp, po);
2622 if (!po) return Flx_to_ZX_inplace(gerepileuptoleaf(av, g));
2623 g = Flx_to_ZX(g);
2624 gerepileall(av, 2, &g, po);
2625 return g;
2626 }
2627 /* p now odd */
2628 q_1 = subiu(powiu(p,f), 1);
2629 N = diviiexact(q_1, p_1);
2630 Lp = odd_prime_divisors(p_1);
2631 for (i=lg(Lp)-1; i; i--) gel(Lp,i) = diviiexact(p_1, gel(Lp,i));
2632 o = factor_pn_1(p,f);
2633 Lq = leafcopy( gel(o, 1) );
2634 for (i = j = 1; i < lg(Lq); i++)
2635 {
2636 if (dvdii(p_1, gel(Lq,i))) continue;
2637 gel(Lq,j++) = diviiexact(N, gel(Lq,i));
2638 }
2639 setlg(Lq, j);
2640 g = gener_FpXQ_i(get_FpX_mod(T), p, p_1, Lp, Lq);
2641 if (!po) g = gerepilecopy(av, g);
2642 else {
2643 *po = mkvec2(q_1, o);
2644 gerepileall(av, 2, &g, po);
2645 }
2646 return g;
2647 }
2648
2649 GEN
gener_FpXQ_local(GEN T,GEN p,GEN L)2650 gener_FpXQ_local(GEN T, GEN p, GEN L)
2651 {
2652 GEN Lp, Lq, p_1 = subiu(p,1), q_1, N, Q;
2653 long f, i, ip, iq, l = lg(L);
2654 T = get_FpX_mod(T);
2655 f = degpol(T);
2656 q_1 = subiu(powiu(p,f), 1);
2657 N = diviiexact(q_1, p_1);
2658
2659 Q = is_pm1(p_1)? gen_1: shifti(p_1,-1);
2660 Lp = cgetg(l, t_VEC); ip = 1;
2661 Lq = cgetg(l, t_VEC); iq = 1;
2662 for (i=1; i < l; i++)
2663 {
2664 GEN a, b, ell = gel(L,i);
2665 if (absequaliu(ell,2)) continue;
2666 a = dvmdii(Q, ell, &b);
2667 if (b == gen_0)
2668 gel(Lp,ip++) = a;
2669 else
2670 gel(Lq,iq++) = diviiexact(N,ell);
2671 }
2672 setlg(Lp, ip);
2673 setlg(Lq, iq);
2674 return gener_FpXQ_i(T, p, p_1, Lp, Lq);
2675 }
2676
2677 /***********************************************************************/
2678 /** **/
2679 /** FpXn **/
2680 /** **/
2681 /***********************************************************************/
2682
2683 INLINE GEN
FpXn_red(GEN a,long n)2684 FpXn_red(GEN a, long n)
2685 { return RgXn_red_shallow(a, n); }
2686
2687 GEN
FpXn_mul(GEN a,GEN b,long n,GEN p)2688 FpXn_mul(GEN a, GEN b, long n, GEN p)
2689 {
2690 return FpX_red(ZXn_mul(a, b, n), p);
2691 }
2692
2693 GEN
FpXn_sqr(GEN a,long n,GEN p)2694 FpXn_sqr(GEN a, long n, GEN p)
2695 {
2696 return FpX_red(ZXn_sqr(a, n), p);
2697 }
2698
2699 /* (f*g) \/ x^n */
2700 static GEN
FpX_mulhigh_i(GEN f,GEN g,long n,GEN p)2701 FpX_mulhigh_i(GEN f, GEN g, long n, GEN p)
2702 {
2703 return FpX_shift(FpX_mul(f,g, p),-n);
2704 }
2705
2706 static GEN
FpXn_mulhigh(GEN f,GEN g,long n2,long n,GEN p)2707 FpXn_mulhigh(GEN f, GEN g, long n2, long n, GEN p)
2708 {
2709 GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
2710 return FpX_add(FpX_mulhigh_i(fl, g, n2, p), FpXn_mul(fh, g, n - n2, p), p);
2711 }
2712
2713 GEN
FpXn_inv(GEN f,long e,GEN p)2714 FpXn_inv(GEN f, long e, GEN p)
2715 {
2716 pari_sp av = avma, av2;
2717 ulong mask;
2718 GEN W, a;
2719 long v = varn(f), n = 1;
2720
2721 if (!signe(f)) pari_err_INV("FpXn_inv",f);
2722 a = Fp_inv(gel(f,2), p);
2723 if (e == 1) return scalarpol(a, v);
2724 else if (e == 2)
2725 {
2726 GEN b;
2727 if (degpol(f) <= 0) return scalarpol(a, v);
2728 b = Fp_neg(gel(f,3),p);
2729 if (signe(b)==0) return scalarpol(a, v);
2730 if (!is_pm1(a)) b = Fp_mul(b, Fp_sqr(a, p), p);
2731 W = deg1pol_shallow(b, a, v);
2732 return gerepilecopy(av, W);
2733 }
2734 W = scalarpol_shallow(Fp_inv(gel(f,2), p),v);
2735 mask = quadratic_prec_mask(e);
2736 av2 = avma;
2737 for (;mask>1;)
2738 {
2739 GEN u, fr;
2740 long n2 = n;
2741 n<<=1; if (mask & 1) n--;
2742 mask >>= 1;
2743 fr = FpXn_red(f, n);
2744 u = FpXn_mul(W, FpXn_mulhigh(fr, W, n2, n, p), n-n2, p);
2745 W = FpX_sub(W, FpX_shift(u, n2), p);
2746 if (gc_needed(av2,2))
2747 {
2748 if(DEBUGMEM>1) pari_warn(warnmem,"FpXn_inv, e = %ld", n);
2749 W = gerepileupto(av2, W);
2750 }
2751 }
2752 return gerepileupto(av, W);
2753 }
2754
2755 GEN
FpXn_expint(GEN h,long e,GEN p)2756 FpXn_expint(GEN h, long e, GEN p)
2757 {
2758 pari_sp av = avma, av2;
2759 long v = varn(h), n=1;
2760 GEN f = pol_1(v), g = pol_1(v);
2761 ulong mask = quadratic_prec_mask(e);
2762 av2 = avma;
2763 for (;mask>1;)
2764 {
2765 GEN u, w;
2766 long n2 = n;
2767 n<<=1; if (mask & 1) n--;
2768 mask >>= 1;
2769 u = FpXn_mul(g, FpX_mulhigh_i(f, FpXn_red(h, n2-1), n2-1, p), n-n2, p);
2770 u = FpX_add(u, FpX_shift(FpXn_red(h, n-1), 1-n2), p);
2771 w = FpXn_mul(f, FpX_integXn(u, n2-1, p), n-n2, p);
2772 f = FpX_add(f, FpX_shift(w, n2), p);
2773 if (mask<=1) break;
2774 u = FpXn_mul(g, FpXn_mulhigh(f, g, n2, n, p), n-n2, p);
2775 g = FpX_sub(g, FpX_shift(u, n2), p);
2776 if (gc_needed(av2,2))
2777 {
2778 if (DEBUGMEM>1) pari_warn(warnmem,"FpXn_exp, e = %ld", n);
2779 gerepileall(av2, 2, &f, &g);
2780 }
2781 }
2782 return gerepileupto(av, f);
2783 }
2784
2785 GEN
FpXn_exp(GEN h,long e,GEN p)2786 FpXn_exp(GEN h, long e, GEN p)
2787 {
2788 if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
2789 pari_err_DOMAIN("FpXn_exp","valuation", "<", gen_1, h);
2790 return FpXn_expint(FpX_deriv(h, p), e, p);
2791 }
2792