Line data Source code
1 : /* Copyright (C) 2004 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 with small coefficients. */
19 :
20 : static GEN
21 1021568730 : get_Flx_red(GEN T, GEN *B)
22 : {
23 1021568730 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 730292 : *B = gel(T,1); return gel(T,2);
25 : }
26 :
27 : /***********************************************************************/
28 : /** Flx **/
29 : /***********************************************************************/
30 : /* Flx objects are defined as follows:
31 : * Let l an ulong. An Flx is a t_VECSMALL:
32 : * x[0] = codeword
33 : * x[1] = evalvarn(variable number) (signe is not stored).
34 : * x[2] = a_0 x[3] = a_1, etc. with 0 <= a_i < l
35 : *
36 : * signe(x) is not valid. Use degpol(x)>0 instead. */
37 : /***********************************************************************/
38 : /** Conversion from Flx **/
39 : /***********************************************************************/
40 :
41 : GEN
42 38523117 : Flx_to_ZX(GEN z)
43 : {
44 38523117 : long i, l = lg(z);
45 38523117 : GEN x = cgetg(l,t_POL);
46 249320491 : for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
47 38523117 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
48 : }
49 :
50 : GEN
51 71583 : Flx_to_FlxX(GEN z, long sv)
52 : {
53 71583 : long i, l = lg(z);
54 71583 : GEN x = cgetg(l,t_POL);
55 278765 : for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
56 71583 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
57 : }
58 :
59 : /* same as Flx_to_ZX, in place */
60 : GEN
61 36056471 : Flx_to_ZX_inplace(GEN z)
62 : {
63 36056471 : long i, l = lg(z);
64 225787691 : for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
65 36056471 : settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
66 : }
67 :
68 : /*Flx_to_Flv=zx_to_zv*/
69 : GEN
70 67382463 : Flx_to_Flv(GEN x, long N)
71 : {
72 67382463 : GEN z = cgetg(N+1,t_VECSMALL);
73 67382463 : long i, l = lg(x)-1;
74 67382463 : x++;
75 725804634 : for (i=1; i<l ; i++) z[i]=x[i];
76 336777511 : for ( ; i<=N; i++) z[i]=0;
77 67382463 : return z;
78 : }
79 :
80 : /*Flv_to_Flx=zv_to_zx*/
81 : GEN
82 26076776 : Flv_to_Flx(GEN x, long sv)
83 : {
84 26076776 : long i, l=lg(x)+1;
85 26076776 : GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
86 26076776 : x--;
87 288934372 : for (i=2; i<l ; i++) z[i]=x[i];
88 26076776 : return Flx_renormalize(z,l);
89 : }
90 :
91 : /*Flm_to_FlxV=zm_to_zxV*/
92 : GEN
93 2772 : Flm_to_FlxV(GEN x, long sv)
94 7455 : { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
95 :
96 : /*FlxC_to_ZXC=zxC_to_ZXC*/
97 : GEN
98 104095 : FlxC_to_ZXC(GEN x)
99 527845 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
100 :
101 : /*FlxC_to_ZXC=zxV_to_ZXV*/
102 : GEN
103 610763 : FlxV_to_ZXV(GEN x)
104 2472894 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
105 :
106 : void
107 3036638 : FlxV_to_ZXV_inplace(GEN v)
108 : {
109 : long i;
110 8064529 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
111 3036638 : }
112 :
113 : /*FlxM_to_ZXM=zxM_to_ZXM*/
114 : GEN
115 2485 : FlxM_to_ZXM(GEN x)
116 8351 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
117 :
118 : GEN
119 410249 : FlxV_to_FlxX(GEN x, long v)
120 : {
121 410249 : long i, l = lg(x)+1;
122 410249 : GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
123 410249 : x--;
124 5925058 : for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
125 410249 : return FlxX_renormalize(z,l);
126 : }
127 :
128 : GEN
129 0 : FlxM_to_FlxXV(GEN x, long v)
130 0 : { pari_APPLY_type(t_COL, FlxV_to_FlxX(gel(x,i), v)) }
131 :
132 : GEN
133 0 : FlxM_Flx_add_shallow(GEN x, GEN y, ulong p)
134 : {
135 0 : long l = lg(x), i, j;
136 0 : GEN z = cgetg(l,t_MAT);
137 :
138 0 : if (l==1) return z;
139 0 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
140 0 : for (i=1; i<l; i++)
141 : {
142 0 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
143 0 : gel(z,i) = zi;
144 0 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
145 0 : gel(zi,i) = Flx_add(gel(zi,i), y, p);
146 : }
147 0 : return z;
148 : }
149 :
150 : /***********************************************************************/
151 : /** Conversion to Flx **/
152 : /***********************************************************************/
153 : /* Take an integer and return a scalar polynomial mod p, with evalvarn=vs */
154 : GEN
155 22911559 : Fl_to_Flx(ulong x, long sv) { return x? mkvecsmall2(sv, x): pol0_Flx(sv); }
156 :
157 : /* a X^d */
158 : GEN
159 947471 : monomial_Flx(ulong a, long d, long vs)
160 : {
161 : GEN P;
162 947471 : if (a==0) return pol0_Flx(vs);
163 947471 : P = const_vecsmall(d+2, 0);
164 947471 : P[1] = vs; P[d+2] = a; return P;
165 : }
166 :
167 : GEN
168 7605485 : Z_to_Flx(GEN x, ulong p, long sv)
169 : {
170 7605485 : long u = umodiu(x,p);
171 7605485 : return u? mkvecsmall2(sv, u): pol0_Flx(sv);
172 : }
173 :
174 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
175 : GEN
176 170264758 : ZX_to_Flx(GEN x, ulong p)
177 : {
178 170264758 : long i, lx = lg(x);
179 170264758 : GEN a = cgetg(lx, t_VECSMALL);
180 170264758 : a[1]=((ulong)x[1])&VARNBITS;
181 1119913808 : for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
182 170264758 : return Flx_renormalize(a,lx);
183 : }
184 :
185 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
186 : GEN
187 6493059 : zx_to_Flx(GEN x, ulong p)
188 : {
189 6493059 : long i, lx = lg(x);
190 6493059 : GEN a = cgetg(lx, t_VECSMALL);
191 6493059 : a[1] = x[1];
192 19973331 : for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
193 6493059 : return Flx_renormalize(a,lx);
194 : }
195 :
196 : ulong
197 80743024 : Rg_to_Fl(GEN x, ulong p)
198 : {
199 80743024 : switch(typ(x))
200 : {
201 56000066 : case t_INT: return umodiu(x, p);
202 480241 : case t_FRAC: {
203 480241 : ulong z = umodiu(gel(x,1), p);
204 480241 : if (!z) return 0;
205 470416 : return Fl_div(z, umodiu(gel(x,2), p), p);
206 : }
207 205973 : case t_PADIC: return padic_to_Fl(x, p);
208 24056744 : case t_INTMOD: {
209 24056744 : GEN q = gel(x,1), a = gel(x,2);
210 24056744 : if (absequaliu(q, p)) return itou(a);
211 0 : if (!dvdiu(q,p)) pari_err_MODULUS("Rg_to_Fl", q, utoipos(p));
212 0 : return umodiu(a, p);
213 : }
214 0 : default: pari_err_TYPE("Rg_to_Fl",x);
215 : return 0; /* LCOV_EXCL_LINE */
216 : }
217 : }
218 :
219 : ulong
220 1716462 : Rg_to_F2(GEN x)
221 : {
222 1716462 : switch(typ(x))
223 : {
224 281252 : case t_INT: return mpodd(x);
225 0 : case t_FRAC:
226 0 : if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
227 0 : return mpodd(gel(x,1));
228 0 : case t_PADIC:
229 0 : if (!absequaliu(padic_p(x),2)) pari_err_OP("",x, mkintmodu(1,2));
230 0 : if (valp(x) < 0) (void)Fl_inv(0,2);
231 0 : return valp(x) & 1;
232 1435210 : case t_INTMOD: {
233 1435210 : GEN q = gel(x,1), a = gel(x,2);
234 1435210 : if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
235 1435210 : return mpodd(a);
236 : }
237 0 : default: pari_err_TYPE("Rg_to_F2",x);
238 : return 0; /* LCOV_EXCL_LINE */
239 : }
240 : }
241 :
242 : GEN
243 2228650 : RgX_to_Flx(GEN x, ulong p)
244 : {
245 2228650 : long i, lx = lg(x);
246 2228650 : GEN a = cgetg(lx, t_VECSMALL);
247 2228650 : a[1]=((ulong)x[1])&VARNBITS;
248 19903308 : for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
249 2228650 : return Flx_renormalize(a,lx);
250 : }
251 :
252 : GEN
253 14 : RgXV_to_FlxV(GEN x, ulong p)
254 602 : { pari_APPLY_type(t_VEC, RgX_to_Flx(gel(x,i), p)) }
255 :
256 : /* If x is a POLMOD, assume modulus is a multiple of T. */
257 : GEN
258 3532252 : Rg_to_Flxq(GEN x, GEN T, ulong p)
259 : {
260 3532252 : long ta, tx = typ(x), v = get_Flx_var(T);
261 : ulong pi;
262 : GEN a, b;
263 3532252 : if (is_const_t(tx))
264 : {
265 3268818 : if (tx == t_FFELT) return FF_to_Flxq(x);
266 2665792 : return Fl_to_Flx(Rg_to_Fl(x, p), v);
267 : }
268 263434 : switch(tx)
269 : {
270 8576 : case t_POLMOD:
271 8576 : b = gel(x,1);
272 8576 : a = gel(x,2); ta = typ(a);
273 8576 : if (is_const_t(ta)) return Fl_to_Flx(Rg_to_Fl(a, p), v);
274 8422 : b = RgX_to_Flx(b, p); if (b[1] != v) break;
275 8422 : a = RgX_to_Flx(a, p); if (Flx_equal(b,T)) return a;
276 0 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
277 0 : if (lgpol(Flx_rem_pre(b,T,p,pi))==0) return Flx_rem_pre(a, T, p, pi);
278 0 : break;
279 254858 : case t_POL:
280 254858 : x = RgX_to_Flx(x,p);
281 254858 : if (x[1] != v) break;
282 254858 : return Flx_rem(x, T, p);
283 0 : case t_RFRAC:
284 0 : a = Rg_to_Flxq(gel(x,1), T,p);
285 0 : b = Rg_to_Flxq(gel(x,2), T,p);
286 0 : return Flxq_div(a,b, T,p);
287 : }
288 0 : pari_err_TYPE("Rg_to_Flxq",x);
289 : return NULL; /* LCOV_EXCL_LINE */
290 : }
291 :
292 : /***********************************************************************/
293 : /** Basic operation on Flx **/
294 : /***********************************************************************/
295 : /* = zx_renormalize. Similar to normalizepol, in place */
296 : GEN
297 2234725994 : Flx_renormalize(GEN /*in place*/ x, long lx)
298 : {
299 : long i;
300 2551706390 : for (i = lx-1; i>1; i--)
301 2436219990 : if (x[i]) break;
302 2234725994 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
303 2234725994 : setlg(x, i+1); return x;
304 : }
305 :
306 : GEN
307 1889643 : Flx_red(GEN z, ulong p)
308 : {
309 1889643 : long i, l = lg(z);
310 1889643 : GEN x = cgetg(l, t_VECSMALL);
311 1889643 : x[1] = z[1];
312 34189034 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
313 1889643 : return Flx_renormalize(x,l);
314 : }
315 :
316 : int
317 26518102 : Flx_equal(GEN V, GEN W)
318 : {
319 26518102 : long l = lg(V);
320 26518102 : if (lg(W) != l) return 0;
321 27172615 : while (--l > 1) /* do not compare variables, V[1] */
322 26272282 : if (V[l] != W[l]) return 0;
323 900333 : return 1;
324 : }
325 :
326 : GEN
327 2666500 : random_Flx(long d1, long vs, ulong p)
328 : {
329 2666500 : long i, d = d1+2;
330 2666500 : GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
331 18456866 : for (i=2; i<d; i++) y[i] = random_Fl(p);
332 2666500 : return Flx_renormalize(y,d);
333 : }
334 :
335 : static GEN
336 7734676 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
337 : {
338 : long i,lz;
339 : GEN z;
340 :
341 7734676 : if (ly>lx) swapspec(x,y, lx,ly);
342 7734676 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
343 115374589 : for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
344 97551794 : for ( ; i<lx; i++) z[i+2] = x[i];
345 7734676 : z[1] = 0; return Flx_renormalize(z, lz);
346 : }
347 :
348 : GEN
349 77210845 : Flx_add(GEN x, GEN y, ulong p)
350 : {
351 : long i,lz;
352 : GEN z;
353 77210845 : long lx=lg(x);
354 77210845 : long ly=lg(y);
355 77210845 : if (ly>lx) swapspec(x,y, lx,ly);
356 77210845 : lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
357 652632542 : for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
358 170790417 : for ( ; i<lx; i++) z[i] = x[i];
359 77210845 : return Flx_renormalize(z, lz);
360 : }
361 :
362 : GEN
363 9950583 : Flx_Fl_add(GEN y, ulong x, ulong p)
364 : {
365 : GEN z;
366 : long lz, i;
367 9950583 : if (!lgpol(y))
368 230308 : return Fl_to_Flx(x,y[1]);
369 9720275 : lz=lg(y);
370 9720275 : z=cgetg(lz,t_VECSMALL);
371 9720275 : z[1]=y[1];
372 9720275 : z[2] = Fl_add(y[2],x,p);
373 47174900 : for(i=3;i<lz;i++)
374 37454625 : z[i] = y[i];
375 9720275 : if (lz==3) z = Flx_renormalize(z,lz);
376 9720275 : return z;
377 : }
378 :
379 : static GEN
380 941456 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
381 : {
382 : long i,lz;
383 : GEN z;
384 :
385 941456 : if (ly <= lx)
386 : {
387 941456 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
388 58430209 : for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
389 1503744 : for ( ; i<lx; i++) z[i+2] = x[i];
390 : }
391 : else
392 : {
393 0 : lz = ly+2; z = cgetg(lz, t_VECSMALL);
394 0 : for (i=0; i<lx; i++) z[i+2] = Fl_sub(x[i],y[i],p);
395 0 : for ( ; i<ly; i++) z[i+2] = Fl_neg(y[i],p);
396 : }
397 941456 : z[1] = 0; return Flx_renormalize(z, lz);
398 : }
399 :
400 : GEN
401 141562515 : Flx_sub(GEN x, GEN y, ulong p)
402 : {
403 141562515 : long i,lz,lx = lg(x), ly = lg(y);
404 : GEN z;
405 :
406 141562515 : if (ly <= lx)
407 : {
408 93223102 : lz = lx; z = cgetg(lz, t_VECSMALL);
409 449508767 : for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
410 183495546 : for ( ; i<lx; i++) z[i] = x[i];
411 : }
412 : else
413 : {
414 48339413 : lz = ly; z = cgetg(lz, t_VECSMALL);
415 247193136 : for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
416 227941817 : for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
417 : }
418 141562515 : z[1]=x[1]; return Flx_renormalize(z, lz);
419 : }
420 :
421 : GEN
422 151980 : Flx_Fl_sub(GEN y, ulong x, ulong p)
423 : {
424 : GEN z;
425 151980 : long lz = lg(y), i;
426 151980 : if (lz==2)
427 513 : return Fl_to_Flx(Fl_neg(x, p),y[1]);
428 151467 : z = cgetg(lz, t_VECSMALL);
429 151467 : z[1] = y[1];
430 151467 : z[2] = Fl_sub(uel(y,2), x, p);
431 754576 : for(i=3; i<lz; i++)
432 603109 : z[i] = y[i];
433 151467 : if (lz==3) z = Flx_renormalize(z,lz);
434 151467 : return z;
435 : }
436 :
437 : static GEN
438 4269121 : Flx_negspec(GEN x, ulong p, long l)
439 : {
440 : long i;
441 4269121 : GEN z = cgetg(l+2, t_VECSMALL) + 2;
442 27892039 : for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
443 4269121 : return z-2;
444 : }
445 :
446 : GEN
447 4269121 : Flx_neg(GEN x, ulong p)
448 : {
449 4269121 : GEN z = Flx_negspec(x+2, p, lgpol(x));
450 4269121 : z[1] = x[1];
451 4269121 : return z;
452 : }
453 :
454 : GEN
455 1892852 : Flx_neg_inplace(GEN x, ulong p)
456 : {
457 1892852 : long i, l = lg(x);
458 56576231 : for (i=2; i<l; i++)
459 54683379 : if (x[i]) x[i] = p - x[i];
460 1892852 : return x;
461 : }
462 :
463 : GEN
464 1182839 : Flx_double(GEN y, ulong p)
465 : {
466 : long i, l;
467 1182839 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
468 9628919 : for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
469 1182839 : return Flx_renormalize(z, l);
470 : }
471 : GEN
472 551631 : Flx_triple(GEN y, ulong p)
473 : {
474 : long i, l;
475 551631 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
476 3850338 : for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
477 551631 : return Flx_renormalize(z, l);
478 : }
479 :
480 : GEN
481 20389790 : Flx_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
482 : {
483 : GEN z;
484 : long i, l;
485 20389790 : if (!x) return pol0_Flx(y[1]);
486 19358748 : z = cgetg_copy(y, &l); z[1] = y[1];
487 19358748 : if (pi==0)
488 : {
489 17070650 : if (HIGHWORD(x | p))
490 0 : for(i=2; i<l; i++) z[i] = Fl_mul(uel(y,i), x, p);
491 : else
492 101611291 : for(i=2; i<l; i++) z[i] = (uel(y,i) * x) % p;
493 : } else
494 18792113 : for(i=2; i<l; i++) z[i] = Fl_mul_pre(uel(y,i), x, p, pi);
495 19358748 : return Flx_renormalize(z, l);
496 : }
497 :
498 : GEN
499 9481826 : Flx_Fl_mul(GEN x, ulong y, ulong p)
500 9481826 : { return Flx_Fl_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
501 :
502 : GEN
503 0 : Flx_convol(GEN x, GEN y, ulong p)
504 : {
505 0 : long lx = lg(x), ly = lg(y), i;
506 : GEN z;
507 0 : if (lx < ly) swapspec(x,y, lx,ly);
508 0 : z = cgetg(ly,t_VECSMALL); z[1] = x[1];
509 0 : for (i=2; i<ly; i++) uel(z,i) = Fl_mul(uel(x,i),uel(y,i), p);
510 0 : return Flx_renormalize(z, ly);
511 : }
512 :
513 : GEN
514 11913433 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
515 : {
516 : GEN z;
517 : long i, l;
518 11913433 : z = cgetg_copy(y, &l); z[1] = y[1];
519 11913433 : if (HIGHWORD(x | p))
520 5414797 : for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
521 : else
522 27076274 : for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
523 11913433 : z[l-1] = 1; return z;
524 : }
525 :
526 : /* Return a*x^n if n>=0 and a\x^(-n) if n<0 */
527 : GEN
528 27608471 : Flx_shift(GEN a, long n)
529 : {
530 27608471 : long i, l = lg(a);
531 : GEN b;
532 27608471 : if (l==2 || !n) return Flx_copy(a);
533 27261085 : if (l+n<=2) return pol0_Flx(a[1]);
534 27043537 : b = cgetg(l+n, t_VECSMALL);
535 27043537 : b[1] = a[1];
536 27043537 : if (n < 0)
537 75586465 : for (i=2-n; i<l; i++) b[i+n] = a[i];
538 : else
539 : {
540 53253451 : for (i=0; i<n; i++) b[2+i] = 0;
541 154488862 : for (i=2; i<l; i++) b[i+n] = a[i];
542 : }
543 27043537 : return b;
544 : }
545 :
546 : GEN
547 67778867 : Flx_normalize(GEN z, ulong p)
548 : {
549 67778867 : long l = lg(z)-1;
550 67778867 : ulong p1 = z[l]; /* leading term */
551 67778867 : if (p1 == 1) return z;
552 11913433 : return Flx_Fl_mul_to_monic(z, Fl_inv(p1,p), p);
553 : }
554 :
555 : /* return (x * X^d) + y. Assume d > 0, shallow if x == 0*/
556 : static GEN
557 3970806 : Flx_addshift(GEN x, GEN y, ulong p, long d)
558 : {
559 3970806 : GEN xd,yd,zd = (GEN)avma;
560 3970806 : long a,lz,ny = lgpol(y), nx = lgpol(x);
561 3970806 : long vs = x[1];
562 3970806 : if (nx == 0) return y;
563 3968963 : x += 2; y += 2; a = ny-d;
564 3968963 : if (a <= 0)
565 : {
566 85113 : lz = (a>nx)? ny+2: nx+d+2;
567 85113 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
568 1733797 : while (xd > x) *--zd = *--xd;
569 85113 : x = zd + a;
570 165659 : while (zd > x) *--zd = 0;
571 : }
572 : else
573 : {
574 3883850 : xd = new_chunk(d); yd = y+d;
575 3883850 : x = Flx_addspec(x,yd,p, nx,a);
576 3883850 : lz = (a>nx)? ny+2: lg(x)+d;
577 143661019 : x += 2; while (xd > x) *--zd = *--xd;
578 : }
579 65216718 : while (yd > y) *--zd = *--yd;
580 3968963 : *--zd = vs;
581 3968963 : *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
582 : }
583 :
584 : /* shift polynomial + GC; do not set evalvarn*/
585 : static GEN
586 649261259 : Flx_shiftip(pari_sp av, GEN x, long v)
587 : {
588 649261259 : long i, lx = lg(x), ly;
589 : GEN y;
590 649261259 : if (!v || lx==2) return gc_leaf(av, x);
591 177214786 : ly = lx + v; /* result length */
592 177214786 : (void)new_chunk(ly); /* check that result fits */
593 177214786 : x += lx; y = (GEN)av;
594 1283482635 : for (i = 2; i<lx; i++) *--y = *--x;
595 711793554 : for (i = 0; i< v; i++) *--y = 0;
596 177214786 : y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
597 177214786 : return gc_const((pari_sp)y, y);
598 : }
599 :
600 : static long
601 2389611307 : get_Fl_threshold(ulong p, long mul, long mul2)
602 : {
603 2389611307 : return SMALL_ULONG(p) ? mul: mul2;
604 : }
605 :
606 : #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
607 : #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
608 : #define LLQUARTWORD(x) ((x) & QUARTMASK)
609 : #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
610 : #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
611 : #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
612 : INLINE long
613 8971924 : maxbitcoeffpol(ulong p, long n)
614 : {
615 8971924 : GEN z = muliu(sqru(p - 1), n);
616 8971924 : long b = expi(z) + 1;
617 : /* only do expensive bit-packing if it saves at least 1 limb */
618 8971924 : if (b <= BITS_IN_QUARTULONG)
619 : {
620 856941 : if (nbits2nlong(n*b) == (n + 3)>>2)
621 110872 : b = BITS_IN_QUARTULONG;
622 : }
623 8114983 : else if (b <= BITS_IN_HALFULONG)
624 : {
625 1696703 : if (nbits2nlong(n*b) == (n + 1)>>1)
626 5682 : b = BITS_IN_HALFULONG;
627 : }
628 : else
629 : {
630 6418280 : long l = lgefint(z) - 2;
631 6418280 : if (nbits2nlong(n*b) == n*l)
632 360226 : b = l*BITS_IN_LONG;
633 : }
634 8971924 : return b;
635 : }
636 :
637 : INLINE ulong
638 3482393121 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
639 : { /* Assume OK_ULONG*/
640 3482393121 : ulong p1 = 0;
641 : long i;
642 16559987766 : for (i=a; i<b; i++)
643 13077594645 : if (y[i])
644 : {
645 10996458190 : p1 += y[i] * x[-i];
646 10996458190 : if (p1 & HIGHBIT) p1 %= p;
647 : }
648 3482393121 : return p1 % p;
649 : }
650 :
651 : INLINE ulong
652 1230369909 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
653 : {
654 1230369909 : ulong p1 = 0;
655 : long i;
656 3938517383 : for (i=a; i<b; i++)
657 2708147474 : if (y[i])
658 2665531050 : p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
659 1230369909 : return p1;
660 : }
661 :
662 : /* assume nx >= ny > 0 */
663 : static GEN
664 356430863 : Flx_mulspec_basecase(GEN x, GEN y, ulong p, ulong pi, long nx, long ny)
665 : {
666 : long i,lz,nz;
667 : GEN z;
668 :
669 356430863 : lz = nx+ny+1; nz = lz-2;
670 356430863 : z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
671 356430863 : if (!pi)
672 : {
673 1180783471 : for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
674 740130501 : for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
675 919048167 : for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
676 : }
677 : else
678 : {
679 336567207 : for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
680 230316949 : for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
681 241871648 : for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
682 : }
683 356430863 : z -= 2; return Flx_renormalize(z, lz);
684 : }
685 :
686 : static GEN
687 14449 : int_to_Flx(GEN z, ulong p)
688 : {
689 14449 : long i, l = lgefint(z);
690 14449 : GEN x = cgetg(l, t_VECSMALL);
691 1236849 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
692 14449 : return Flx_renormalize(x, l);
693 : }
694 :
695 : INLINE GEN
696 11437 : Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
697 : {
698 11437 : GEN z=muliispec(a,b,na,nb);
699 11437 : return int_to_Flx(z,p);
700 : }
701 :
702 : static GEN
703 556176 : Flx_to_int_halfspec(GEN a, long na)
704 : {
705 : long j;
706 556176 : long n = (na+1)>>1UL;
707 556176 : GEN V = cgetipos(2+n);
708 : GEN w;
709 1513393 : for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
710 957217 : *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
711 556176 : if (j<na)
712 407228 : *w = a[j];
713 556176 : return V;
714 : }
715 :
716 : static GEN
717 733977 : int_to_Flx_half(GEN z, ulong p)
718 : {
719 : long i;
720 733977 : long lx = (lgefint(z)-2)*2+2;
721 733977 : GEN w, x = cgetg(lx, t_VECSMALL);
722 2293472 : for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
723 : {
724 1559495 : x[i] = LOWWORD((ulong)*w)%p;
725 1559495 : x[i+1] = HIGHWORD((ulong)*w)%p;
726 : }
727 733977 : return Flx_renormalize(x, lx);
728 : }
729 :
730 : static GEN
731 5543 : Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
732 : {
733 5543 : GEN A = Flx_to_int_halfspec(a,na);
734 5543 : GEN B = Flx_to_int_halfspec(b,nb);
735 5543 : GEN z = mulii(A,B);
736 5543 : return int_to_Flx_half(z,p);
737 : }
738 :
739 : static GEN
740 211232 : Flx_to_int_quartspec(GEN a, long na)
741 : {
742 : long j;
743 211232 : long n = (na+3)>>2UL;
744 211232 : GEN V = cgetipos(2+n);
745 : GEN w;
746 4625177 : for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
747 4413945 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
748 211232 : switch (na-j)
749 : {
750 118717 : case 3:
751 118717 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
752 118717 : break;
753 35407 : case 2:
754 35407 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
755 35407 : break;
756 28553 : case 1:
757 28553 : *w = a[j];
758 28553 : break;
759 28555 : case 0:
760 28555 : break;
761 : }
762 211232 : return V;
763 : }
764 :
765 : static GEN
766 110872 : int_to_Flx_quart(GEN z, ulong p)
767 : {
768 : long i;
769 110872 : long lx = (lgefint(z)-2)*4+2;
770 110872 : GEN w, x = cgetg(lx, t_VECSMALL);
771 5127515 : for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
772 : {
773 5016643 : x[i] = LLQUARTWORD((ulong)*w)%p;
774 5016643 : x[i+1] = HLQUARTWORD((ulong)*w)%p;
775 5016643 : x[i+2] = LHQUARTWORD((ulong)*w)%p;
776 5016643 : x[i+3] = HHQUARTWORD((ulong)*w)%p;
777 : }
778 110872 : return Flx_renormalize(x, lx);
779 : }
780 :
781 : static GEN
782 100360 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
783 : {
784 100360 : GEN A = Flx_to_int_quartspec(a,na);
785 100360 : GEN B = Flx_to_int_quartspec(b,nb);
786 100360 : GEN z = mulii(A,B);
787 100360 : return int_to_Flx_quart(z,p);
788 : }
789 :
790 : /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
791 : static GEN
792 681479 : Flx_eval2BILspec(GEN x, long k, long l)
793 : {
794 681479 : long i, lz = k*l, ki;
795 681479 : GEN pz = cgetipos(2+lz);
796 19498113 : for (i=0; i < lz; i++)
797 18816634 : *int_W(pz,i) = 0UL;
798 10089796 : for (i=0, ki=0; i<l; i++, ki+=k)
799 9408317 : *int_W(pz,ki) = x[i];
800 681479 : return int_normalize(pz,0);
801 : }
802 :
803 : static GEN
804 348707 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
805 : {
806 348707 : long i, offset, lm = lgefint(x)-2, l = d+3;
807 348707 : ulong pi = get_Fl_red(p);
808 348707 : GEN pol = cgetg(l, t_VECSMALL);
809 348707 : pol[1] = 0;
810 9536550 : for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
811 9187843 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
812 348707 : if (offset < lm)
813 273562 : pol[i+2] = (*int_W(x,offset)) % p;
814 348707 : return Flx_renormalize(pol,l);
815 : }
816 :
817 : static GEN
818 0 : Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
819 : {
820 0 : long i, offset, lm = lgefint(x)-2, l = d+3;
821 0 : ulong pi = get_Fl_red(p);
822 0 : GEN pol = cgetg(l, t_VECSMALL);
823 0 : pol[1] = 0;
824 0 : for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
825 0 : pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
826 0 : *int_W(x,offset), p, pi);
827 0 : if (offset+1 < lm)
828 0 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
829 0 : else if (offset < lm)
830 0 : pol[i+2] = (*int_W(x,offset)) % p;
831 0 : return Flx_renormalize(pol,l);
832 : }
833 :
834 : static GEN
835 345777 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
836 : {
837 345777 : return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
838 : }
839 :
840 : static GEN
841 332313 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
842 : {
843 332313 : pari_sp av = avma;
844 332313 : GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
845 332313 : return gc_upto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
846 : }
847 :
848 : static GEN
849 22387649 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
850 : GEN y;
851 : long i;
852 22387649 : if (l == 0)
853 3794548 : return gen_0;
854 18593101 : y = cgetg(l + 1, t_VECSMALL);
855 894011326 : for(i = 1; i <= l; i++)
856 875418225 : y[i] = x[l - i];
857 18593101 : return nv_fromdigits_2k(y, b);
858 : }
859 :
860 : /* assume b < BITS_IN_LONG */
861 : static GEN
862 6446145 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
863 6446145 : GEN v = binary_2k_nv(z, b), x;
864 6446145 : long i, l = lg(v) + 1;
865 6446145 : x = cgetg(l, t_VECSMALL);
866 687499108 : for (i = 2; i < l; i++)
867 681052963 : x[i] = v[l - i] % p;
868 6446145 : return Flx_renormalize(x, l);
869 : }
870 :
871 : static GEN
872 5979729 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
873 5979729 : GEN v = binary_2k(z, b), x, y;
874 5979729 : long i, l = lg(v) + 1, ly;
875 5979729 : x = cgetg(l, t_VECSMALL);
876 253053653 : for (i = 2; i < l; i++) {
877 247073924 : y = gel(v, l - i);
878 247073924 : ly = lgefint(y);
879 247073924 : switch (ly) {
880 6290663 : case 2: x[i] = 0; break;
881 31496598 : case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
882 193253268 : case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
883 32066790 : case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
884 16033395 : *int_W_lg(y, 0, ly), p, pi); break;
885 0 : default: x[i] = umodiu(gel(v, l - i), p);
886 : }
887 : }
888 5979729 : return Flx_renormalize(x, l);
889 : }
890 :
891 : static GEN
892 7768105 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
893 : {
894 : GEN C, D;
895 7768105 : pari_sp av = avma;
896 7768105 : A = kron_pack_Flx_spec_bits(A, b, lA);
897 7768105 : B = kron_pack_Flx_spec_bits(B, b, lB);
898 7768105 : C = gc_INT(av, mulii(A, B));
899 7768105 : if (b < BITS_IN_LONG)
900 2204096 : D = kron_unpack_Flx_bits_narrow(C, b, p);
901 : else
902 : {
903 5564009 : ulong pi = get_Fl_red(p);
904 5564009 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
905 : }
906 7768105 : return D;
907 : }
908 :
909 : static GEN
910 727039 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
911 : {
912 : GEN C, D;
913 727039 : A = kron_pack_Flx_spec_bits(A, b, lA);
914 727039 : C = sqri(A);
915 727039 : if (b < BITS_IN_LONG)
916 478671 : D = kron_unpack_Flx_bits_narrow(C, b, p);
917 : else
918 : {
919 248368 : ulong pi = get_Fl_red(p);
920 248368 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
921 : }
922 727039 : return D;
923 : }
924 :
925 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
926 : * b+2 were sent instead. na, nb = number of terms of a, b.
927 : * Only c, c0, c1, c2 are genuine GEN.
928 : */
929 : static GEN
930 399610311 : Flx_mulspec(GEN a, GEN b, ulong p, ulong pi, long na, long nb)
931 : {
932 : GEN a0,c,c0;
933 399610311 : long n0, n0a, i, v = 0;
934 : pari_sp av;
935 :
936 506397598 : while (na && !a[0]) { a++; na--; v++; }
937 587153769 : while (nb && !b[0]) { b++; nb--; v++; }
938 399610311 : if (na < nb) swapspec(a,b, na,nb);
939 399610311 : if (!nb) return pol0_Flx(0);
940 :
941 366596312 : av = avma;
942 366596312 : if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
943 : {
944 8217758 : long m = maxbitcoeffpol(p,nb);
945 8217758 : switch (m)
946 : {
947 100360 : case BITS_IN_QUARTULONG:
948 100360 : return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
949 5543 : case BITS_IN_HALFULONG:
950 5543 : return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
951 11437 : case BITS_IN_LONG:
952 11437 : return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
953 332313 : case 2*BITS_IN_LONG:
954 332313 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
955 0 : case 3*BITS_IN_LONG:
956 0 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
957 7768105 : default:
958 7768105 : return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
959 : }
960 : }
961 358378554 : if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
962 356430863 : return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,pi,na,nb), v);
963 1947691 : i=(na>>1); n0=na-i; na=i;
964 1947691 : a0=a+n0; n0a=n0;
965 2764332 : while (n0a && !a[n0a-1]) n0a--;
966 :
967 1947691 : if (nb > n0)
968 : {
969 : GEN b0,c1,c2;
970 : long n0b;
971 :
972 1892852 : nb -= n0; b0 = b+n0; n0b = n0;
973 3023871 : while (n0b && !b[n0b-1]) n0b--;
974 1892852 : c = Flx_mulspec(a,b,p,pi,n0a,n0b);
975 1892852 : c0 = Flx_mulspec(a0,b0,p,pi,na,nb);
976 :
977 1892852 : c2 = Flx_addspec(a0,a,p,na,n0a);
978 1892852 : c1 = Flx_addspec(b0,b,p,nb,n0b);
979 :
980 1892852 : c1 = Flx_mul_pre(c1,c2,p,pi);
981 1892852 : c2 = Flx_add(c0,c,p);
982 :
983 1892852 : c2 = Flx_neg_inplace(c2,p);
984 1892852 : c2 = Flx_add(c1,c2,p);
985 1892852 : c0 = Flx_addshift(c0,c2, p, n0);
986 : }
987 : else
988 : {
989 54839 : c = Flx_mulspec(a,b,p,pi,n0a,nb);
990 54839 : c0 = Flx_mulspec(a0,b,p,pi,na,nb);
991 : }
992 1947691 : c0 = Flx_addshift(c0,c,p,n0);
993 1947691 : return Flx_shiftip(av,c0, v);
994 : }
995 :
996 : GEN
997 393575947 : Flx_mul_pre(GEN x, GEN y, ulong p, ulong pi)
998 : {
999 393575947 : GEN z = Flx_mulspec(x+2,y+2,p, pi, lgpol(x),lgpol(y));
1000 393575947 : z[1] = x[1]; return z;
1001 : }
1002 : GEN
1003 26819185 : Flx_mul(GEN x, GEN y, ulong p)
1004 26819185 : { return Flx_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1005 :
1006 : static GEN
1007 281845640 : Flx_sqrspec_basecase(GEN x, ulong p, ulong pi, long nx)
1008 : {
1009 : long i, lz, nz;
1010 : ulong p1;
1011 : GEN z;
1012 :
1013 281845640 : if (!nx) return pol0_Flx(0);
1014 281845640 : lz = (nx << 1) + 1, nz = lz-2;
1015 281845640 : z = cgetg(lz, t_VECSMALL) + 2;
1016 281845640 : if (!pi)
1017 : {
1018 217398404 : z[0] = x[0]*x[0]%p;
1019 931216851 : for (i=1; i<nx; i++)
1020 : {
1021 713818447 : p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1022 713818447 : p1 <<= 1;
1023 713818447 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1024 713818447 : z[i] = p1 % p;
1025 : }
1026 931216851 : for ( ; i<nz; i++)
1027 : {
1028 713818447 : p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1029 713818447 : p1 <<= 1;
1030 713818447 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1031 713818447 : z[i] = p1 % p;
1032 : }
1033 : }
1034 : else
1035 : {
1036 64447236 : z[0] = Fl_sqr_pre(x[0], p, pi);
1037 417297627 : for (i=1; i<nx; i++)
1038 : {
1039 352850391 : p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1040 352850391 : p1 = Fl_add(p1, p1, p);
1041 352850391 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1042 352850391 : z[i] = p1;
1043 : }
1044 417297627 : for ( ; i<nz; i++)
1045 : {
1046 352850391 : p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1047 352850391 : p1 = Fl_add(p1, p1, p);
1048 352850391 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1049 352850391 : z[i] = p1;
1050 : }
1051 : }
1052 281845640 : z -= 2; return Flx_renormalize(z, lz);
1053 : }
1054 :
1055 : static GEN
1056 3012 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
1057 : {
1058 3012 : GEN z=sqrispec(a,na);
1059 3012 : return int_to_Flx(z,p);
1060 : }
1061 :
1062 : static GEN
1063 139 : Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
1064 : {
1065 139 : GEN z = sqri(Flx_to_int_halfspec(a,na));
1066 139 : return int_to_Flx_half(z,p);
1067 : }
1068 :
1069 : static GEN
1070 10512 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1071 : {
1072 10512 : GEN z = sqri(Flx_to_int_quartspec(a,na));
1073 10512 : return int_to_Flx_quart(z,p);
1074 : }
1075 :
1076 : static GEN
1077 13464 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1078 : {
1079 13464 : pari_sp av = avma;
1080 13464 : GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1081 13464 : return gc_upto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1082 : }
1083 :
1084 : static GEN
1085 282911709 : Flx_sqrspec(GEN a, ulong p, ulong pi, long na)
1086 : {
1087 : GEN a0, c, c0;
1088 282911709 : long n0, n0a, i, v = 0, m;
1089 : pari_sp av;
1090 :
1091 404947111 : while (na && !a[0]) { a++; na--; v += 2; }
1092 282911709 : if (!na) return pol0_Flx(0);
1093 :
1094 282664947 : av = avma;
1095 282664947 : if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1096 : {
1097 754166 : m = maxbitcoeffpol(p,na);
1098 754166 : switch(m)
1099 : {
1100 10512 : case BITS_IN_QUARTULONG:
1101 10512 : return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
1102 139 : case BITS_IN_HALFULONG:
1103 139 : return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
1104 3012 : case BITS_IN_LONG:
1105 3012 : return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1106 13464 : case 2*BITS_IN_LONG:
1107 13464 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
1108 0 : case 3*BITS_IN_LONG:
1109 0 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
1110 727039 : default:
1111 727039 : return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1112 : }
1113 : }
1114 281910781 : if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1115 281845640 : return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,pi,na), v);
1116 65141 : i=(na>>1); n0=na-i; na=i;
1117 65141 : a0=a+n0; n0a=n0;
1118 80228 : while (n0a && !a[n0a-1]) n0a--;
1119 :
1120 65141 : c = Flx_sqrspec(a,p,pi,n0a);
1121 65141 : c0= Flx_sqrspec(a0,p,pi,na);
1122 65141 : if (p == 2) n0 *= 2;
1123 : else
1124 : {
1125 65122 : GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1126 65122 : t = Flx_sqr_pre(t,p,pi);
1127 65122 : c1= Flx_add(c0,c, p);
1128 65122 : c1= Flx_sub(t, c1, p);
1129 65122 : c0 = Flx_addshift(c0,c1,p,n0);
1130 : }
1131 65141 : c0 = Flx_addshift(c0,c,p,n0);
1132 65141 : return Flx_shiftip(av,c0,v);
1133 : }
1134 :
1135 : GEN
1136 282781427 : Flx_sqr_pre(GEN x, ulong p, ulong pi)
1137 : {
1138 282781427 : GEN z = Flx_sqrspec(x+2,p, pi, lgpol(x));
1139 282781427 : z[1] = x[1]; return z;
1140 : }
1141 : GEN
1142 394177 : Flx_sqr(GEN x, ulong p)
1143 394177 : { return Flx_sqr_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1144 :
1145 : GEN
1146 4020 : Flx_powu_pre(GEN x, ulong n, ulong p, ulong pi)
1147 : {
1148 4020 : GEN y = pol1_Flx(x[1]), z;
1149 : ulong m;
1150 4020 : if (n == 0) return y;
1151 4020 : m = n; z = x;
1152 : for (;;)
1153 : {
1154 13037 : if (m&1UL) y = Flx_mul_pre(y,z, p, pi);
1155 13037 : m >>= 1; if (!m) return y;
1156 9017 : z = Flx_sqr_pre(z, p, pi);
1157 : }
1158 : }
1159 : GEN
1160 0 : Flx_powu(GEN x, ulong n, ulong p)
1161 : {
1162 0 : if (n == 0) return pol1_Flx(x[1]);
1163 0 : return Flx_powu_pre(x, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
1164 : }
1165 :
1166 : GEN
1167 13965 : Flx_halve(GEN y, ulong p)
1168 : {
1169 : GEN z;
1170 : long i, l;
1171 13965 : z = cgetg_copy(y, &l); z[1] = y[1];
1172 57625 : for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
1173 13965 : return z;
1174 : }
1175 :
1176 : static GEN
1177 7414319 : Flx_recipspec(GEN x, long l, long n)
1178 : {
1179 : long i;
1180 7414319 : GEN z=cgetg(n+2,t_VECSMALL)+2;
1181 124654067 : for(i=0; i<l; i++)
1182 117239748 : z[n-i-1] = x[i];
1183 16161185 : for( ; i<n; i++)
1184 8746866 : z[n-i-1] = 0;
1185 7414319 : return Flx_renormalize(z-2,n+2);
1186 : }
1187 :
1188 : GEN
1189 0 : Flx_recip(GEN x)
1190 : {
1191 0 : GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
1192 0 : z[1]=x[1];
1193 0 : return z;
1194 : }
1195 :
1196 : /* Return P(x * h) */
1197 : GEN
1198 0 : Flx_unscale(GEN P, ulong h, ulong p)
1199 : {
1200 : long i, l;
1201 0 : ulong hi = 1UL;
1202 0 : GEN Q = cgetg_copy(P, &l);
1203 0 : Q[1] = P[1];
1204 0 : if (l == 2) return Q;
1205 0 : uel(Q,2) = uel(P,2);
1206 0 : for (i=3; i<l; i++)
1207 : {
1208 0 : hi = Fl_mul(hi, h, p);
1209 0 : uel(Q,i) = Fl_mul(uel(P,i), hi, p);
1210 : }
1211 0 : return Q;
1212 : }
1213 : /* Return h^degpol(P) P(x / h) */
1214 : GEN
1215 1117 : Flx_rescale(GEN P, ulong h, ulong p)
1216 : {
1217 1117 : long i, l = lg(P);
1218 1117 : GEN Q = cgetg(l,t_VECSMALL);
1219 1117 : ulong hi = h;
1220 1117 : Q[l-1] = P[l-1];
1221 12538 : for (i=l-2; i>=2; i--)
1222 : {
1223 12538 : Q[i] = Fl_mul(P[i], hi, p);
1224 12538 : if (i == 2) break;
1225 11421 : hi = Fl_mul(hi,h, p);
1226 : }
1227 1117 : Q[1] = P[1]; return Q;
1228 : }
1229 :
1230 : /* x/polrecip(P)+O(x^n); allow pi = 0 */
1231 : static GEN
1232 134967 : Flx_invBarrett_basecase(GEN T, ulong p, ulong pi)
1233 : {
1234 134967 : long i, l=lg(T)-1, lr=l-1, k;
1235 134967 : GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1236 134967 : r[2] = 1;
1237 134967 : if (!pi)
1238 784007 : for (i=3;i<lr;i++)
1239 : {
1240 776645 : ulong u = uel(T, l-i+2);
1241 46900361 : for (k=3; k<i; k++)
1242 46123716 : { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1243 776645 : r[i] = Fl_neg(u % p, p);
1244 : }
1245 : else
1246 2122222 : for (i=3;i<lr;i++)
1247 : {
1248 1994617 : ulong u = Fl_neg(uel(T,l-i+2), p);
1249 59830442 : for (k=3; k<i; k++)
1250 : {
1251 57835825 : ulong t = Fl_neg(uel(T,l-i+k), p);
1252 57835825 : u = Fl_addmul_pre(u, t, uel(r,k), p, pi);
1253 : }
1254 1994617 : r[i] = u;
1255 : }
1256 134967 : return Flx_renormalize(r,lr);
1257 : }
1258 :
1259 : /* Return new lgpol */
1260 : static long
1261 2236844 : Flx_lgrenormalizespec(GEN x, long lx)
1262 : {
1263 : long i;
1264 7711508 : for (i = lx-1; i>=0; i--)
1265 7710602 : if (x[i]) break;
1266 2236844 : return i+1;
1267 : }
1268 : /* allow pi = 0 */
1269 : static GEN
1270 24044 : Flx_invBarrett_Newton(GEN T, ulong p, ulong pi)
1271 : {
1272 24044 : long nold, lx, lz, lq, l = degpol(T), lQ;
1273 24044 : GEN q, y, z, x = zero_zv(l+1) + 2;
1274 24044 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1275 : pari_sp av;
1276 :
1277 24044 : y = T+2;
1278 24044 : q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1279 24044 : av = avma;
1280 : /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1281 :
1282 : /* initialize */
1283 24044 : x[0] = Fl_inv(q[0], p);
1284 24044 : if (lQ>1 && q[1])
1285 5861 : {
1286 5861 : ulong u = q[1];
1287 5861 : if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1288 5861 : x[1] = p - u; lx = 2;
1289 : }
1290 : else
1291 18183 : lx = 1;
1292 24044 : nold = 1;
1293 166034 : for (; mask > 1; set_avma(av))
1294 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1295 141990 : long i, lnew, nnew = nold << 1;
1296 :
1297 141990 : if (mask & 1) nnew--;
1298 141990 : mask >>= 1;
1299 :
1300 141990 : lnew = nnew + 1;
1301 141990 : lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1302 141990 : z = Flx_mulspec(x, q, p, pi, lx, lq); /* FIXME: high product */
1303 141990 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1304 141990 : z += 2;
1305 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1306 327752 : for (i = nold; i < lz; i++) if (z[i]) break;
1307 141990 : nold = nnew;
1308 141990 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1309 :
1310 : /* z + i represents (x*q - 1) / t^i */
1311 106094 : lz = Flx_lgrenormalizespec (z+i, lz-i);
1312 106094 : z = Flx_mulspec(x, z+i, p, pi, lx, lz); /* FIXME: low product */
1313 106094 : lz = lgpol(z); z += 2;
1314 106094 : if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1315 :
1316 106094 : lx = lz+ i;
1317 106094 : y = x + i; /* x -= z * t^i, in place */
1318 1106928 : for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1319 : }
1320 24044 : x -= 2; setlg(x, lx + 2); x[1] = T[1];
1321 24044 : return x;
1322 : }
1323 :
1324 : /* allow pi = 0 */
1325 : static GEN
1326 160314 : Flx_invBarrett_pre(GEN T, ulong p, ulong pi)
1327 : {
1328 160314 : pari_sp ltop = avma;
1329 160314 : long l = lgpol(T);
1330 : GEN r;
1331 160314 : if (l < 3) return pol0_Flx(T[1]);
1332 159011 : if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1333 : {
1334 134967 : ulong c = T[l+1];
1335 134967 : if (c != 1)
1336 : {
1337 98118 : ulong ci = Fl_inv(c,p);
1338 98118 : T = Flx_Fl_mul_pre(T, ci, p, pi);
1339 98118 : r = Flx_invBarrett_basecase(T, p, pi);
1340 98118 : r = Flx_Fl_mul_pre(r, ci, p, pi);
1341 : }
1342 : else
1343 36849 : r = Flx_invBarrett_basecase(T, p, pi);
1344 : }
1345 : else
1346 24044 : r = Flx_invBarrett_Newton(T, p, pi);
1347 159011 : return gc_leaf(ltop, r);
1348 : }
1349 : GEN
1350 0 : Flx_invBarrett(GEN T, ulong p)
1351 0 : { return Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1352 :
1353 : /* allow pi = 0 */
1354 : GEN
1355 98338446 : Flx_get_red_pre(GEN T, ulong p, ulong pi)
1356 : {
1357 98338446 : if (typ(T)!=t_VECSMALL
1358 98301254 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1359 : Flx_BARRETT2_LIMIT))
1360 98330417 : return T;
1361 8029 : retmkvec2(Flx_invBarrett_pre(T, p, pi),T);
1362 : }
1363 : GEN
1364 14550265 : Flx_get_red(GEN T, ulong p)
1365 : {
1366 14550265 : if (typ(T)!=t_VECSMALL
1367 14550158 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1368 : Flx_BARRETT2_LIMIT))
1369 14544871 : return T;
1370 5394 : retmkvec2(Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)),T);
1371 : }
1372 :
1373 : /* separate from Flx_divrem for maximal speed. */
1374 : static GEN
1375 815529340 : Flx_rem_basecase(GEN x, GEN y, ulong p, ulong pi)
1376 : {
1377 : pari_sp av;
1378 : GEN z, c;
1379 : long dx,dy,dy1,dz,i,j;
1380 : ulong p1,inv;
1381 815529340 : long vs=x[1];
1382 :
1383 815529340 : dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1384 779390587 : dx = degpol(x);
1385 779390587 : dz = dx-dy; if (dz < 0) return Flx_copy(x);
1386 779390587 : x += 2; y += 2;
1387 779390587 : inv = y[dy];
1388 779390587 : if (inv != 1UL) inv = Fl_inv(inv,p);
1389 932516058 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1390 :
1391 779390581 : c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1392 779390581 : z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1393 :
1394 779390581 : if (!pi)
1395 : {
1396 494986694 : z[dz] = (inv*x[dx]) % p;
1397 1842765221 : for (i=dx-1; i>=dy; --i)
1398 : {
1399 1347778527 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1400 10684676435 : for (j=i-dy1; j<=i && j<=dz; j++)
1401 : {
1402 9336897908 : p1 += z[j]*y[i-j];
1403 9336897908 : if (p1 & HIGHBIT) p1 %= p;
1404 : }
1405 1347778527 : p1 %= p;
1406 1347778527 : z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1407 : }
1408 3389567223 : for (i=0; i<dy; i++)
1409 : {
1410 2894580529 : p1 = z[0]*y[i];
1411 14866059210 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1412 : {
1413 11971478681 : p1 += z[j]*y[i-j];
1414 11971478681 : if (p1 & HIGHBIT) p1 %= p;
1415 : }
1416 2894580529 : c[i] = Fl_sub(x[i], p1%p, p);
1417 : }
1418 : }
1419 : else
1420 : {
1421 284403887 : z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1422 873565007 : for (i=dx-1; i>=dy; --i)
1423 : {
1424 589161120 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1425 2496183719 : for (j=i-dy1; j<=i && j<=dz; j++)
1426 1907022599 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1427 589161120 : z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1428 : }
1429 2094520114 : for (i=0; i<dy; i++)
1430 : {
1431 1810116227 : p1 = Fl_mul_pre(z[0],y[i],p,pi);
1432 4919012910 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1433 3108896683 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1434 1810116227 : c[i] = Fl_sub(x[i], p1, p);
1435 : }
1436 : }
1437 951075641 : i = dy-1; while (i>=0 && !c[i]) i--;
1438 779390581 : set_avma(av); return Flx_renormalize(c-2, i+3);
1439 : }
1440 :
1441 : /* as FpX_divrem but working only on ulong types.
1442 : * if relevant, *pr is the last object on stack */
1443 : static GEN
1444 61532362 : Flx_divrem_basecase(GEN x, GEN y, ulong p, ulong pi, GEN *pr)
1445 : {
1446 : GEN z,q,c;
1447 : long dx,dy,dy1,dz,i,j;
1448 : ulong p1,inv;
1449 61532362 : long sv=x[1];
1450 :
1451 61532362 : dy = degpol(y);
1452 61532362 : if (dy<0) pari_err_INV("Flx_divrem",y);
1453 61532362 : if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p, pi);
1454 61531961 : if (!dy)
1455 : {
1456 7141049 : if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1457 7141049 : if (y[2] == 1UL) return Flx_copy(x);
1458 5127147 : return Flx_Fl_mul_pre(x, Fl_inv(y[2], p), p, pi);
1459 : }
1460 54390912 : dx = degpol(x);
1461 54390912 : dz = dx-dy;
1462 54390912 : if (dz < 0)
1463 : {
1464 1062385 : q = pol0_Flx(sv);
1465 1062385 : if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1466 1062385 : return q;
1467 : }
1468 53328527 : x += 2;
1469 53328527 : y += 2;
1470 53328527 : z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1471 53328527 : inv = uel(y, dy);
1472 53328527 : if (inv != 1UL) inv = Fl_inv(inv,p);
1473 78687355 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1474 :
1475 53328527 : if (SMALL_ULONG(p))
1476 : {
1477 51433320 : z[dz] = (inv*x[dx]) % p;
1478 131170505 : for (i=dx-1; i>=dy; --i)
1479 : {
1480 79737185 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1481 262072122 : for (j=i-dy1; j<=i && j<=dz; j++)
1482 : {
1483 182334937 : p1 += z[j]*y[i-j];
1484 182334937 : if (p1 & HIGHBIT) p1 %= p;
1485 : }
1486 79737185 : p1 %= p;
1487 79737185 : z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1488 : }
1489 : }
1490 : else
1491 : {
1492 1895207 : z[dz] = Fl_mul(inv, x[dx], p);
1493 9302516 : for (i=dx-1; i>=dy; --i)
1494 : { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1495 7407309 : p1 = p - uel(x,i);
1496 26553917 : for (j=i-dy1; j<=i && j<=dz; j++)
1497 19146608 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1498 7407309 : z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1499 : }
1500 : }
1501 53328527 : q = Flx_renormalize(z-2, dz+3);
1502 53328527 : if (!pr) return q;
1503 :
1504 25701686 : c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1505 25701686 : if (SMALL_ULONG(p))
1506 : {
1507 213643893 : for (i=0; i<dy; i++)
1508 : {
1509 189592259 : p1 = (ulong)z[0]*y[i];
1510 453086646 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1511 : {
1512 263494387 : p1 += (ulong)z[j]*y[i-j];
1513 263494387 : if (p1 & HIGHBIT) p1 %= p;
1514 : }
1515 189592259 : c[i] = Fl_sub(x[i], p1%p, p);
1516 : }
1517 : }
1518 : else
1519 : {
1520 16113377 : for (i=0; i<dy; i++)
1521 : {
1522 14463325 : p1 = Fl_mul(z[0],y[i],p);
1523 50363696 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1524 35900371 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1525 14463325 : c[i] = Fl_sub(x[i], p1, p);
1526 : }
1527 : }
1528 34883555 : i=dy-1; while (i>=0 && !c[i]) i--;
1529 25701686 : c = Flx_renormalize(c-2, i+3);
1530 25701686 : if (pr == ONLY_DIVIDES)
1531 455 : { if (lg(c) != 2) return NULL; }
1532 : else
1533 25701231 : *pr = c;
1534 25701539 : return q;
1535 : }
1536 :
1537 : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1538 : * and mg is the Barrett inverse of T. */
1539 : static GEN
1540 949442 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1541 : {
1542 : GEN q, r;
1543 949442 : long lt = degpol(T); /*We discard the leading term*/
1544 : long ld, lm, lT, lmg;
1545 949442 : ld = l-lt;
1546 949442 : lm = minss(ld, lgpol(mg));
1547 949442 : lT = Flx_lgrenormalizespec(T+2,lt);
1548 949442 : lmg = Flx_lgrenormalizespec(mg+2,lm);
1549 949442 : q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1550 949442 : q = Flx_mulspec(q+2,mg+2,p,pi,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1551 949442 : q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1552 949442 : if (!pr) return q;
1553 941456 : r = Flx_mulspec(q+2,T+2,p,pi,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1554 941456 : r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1555 941456 : if (pr == ONLY_REM) return r;
1556 442136 : *pr = r; return q;
1557 : }
1558 :
1559 : static GEN
1560 635049 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1561 : {
1562 635049 : GEN q = NULL, r = Flx_copy(x);
1563 635049 : long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1564 : long i;
1565 635049 : if (l <= lt)
1566 : {
1567 0 : if (pr == ONLY_REM) return Flx_copy(x);
1568 0 : if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1569 0 : if (pr) *pr = Flx_copy(x);
1570 0 : return pol0_Flx(v);
1571 : }
1572 635049 : if (lt <= 1)
1573 1303 : return Flx_divrem_basecase(x,T,p,pi,pr);
1574 633746 : if (pr != ONLY_REM && l>lm)
1575 29291 : { q = zero_zv(l-lt+1); q[1] = T[1]; }
1576 951087 : while (l>lm)
1577 : {
1578 317341 : GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,pi,&zr);
1579 317341 : long lz = lgpol(zr);
1580 317341 : if (pr != ONLY_REM)
1581 : {
1582 65542 : long lq = lgpol(zq);
1583 947971 : for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1584 : }
1585 4617241 : for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1586 317341 : l = l-lm+lz;
1587 : }
1588 633746 : if (pr == ONLY_REM)
1589 : {
1590 499363 : if (l > lt)
1591 499320 : r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi,ONLY_REM);
1592 : else
1593 43 : r = Flx_renormalize(r, l+2);
1594 499363 : r[1] = v; return r;
1595 : }
1596 134383 : if (l > lt)
1597 : {
1598 132781 : GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi, pr ? &r: NULL);
1599 132781 : if (!q) q = zq;
1600 : else
1601 : {
1602 27689 : long lq = lgpol(zq);
1603 166223 : for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1604 : }
1605 : }
1606 1602 : else if (pr)
1607 1550 : r = Flx_renormalize(r, l+2);
1608 134383 : q[1] = v; q = Flx_renormalize(q, lg(q));
1609 134383 : if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1610 134383 : if (pr) { r[1] = v; *pr = r; }
1611 134383 : return q;
1612 : }
1613 :
1614 : /* allow pi = 0 (SMALL_ULONG) */
1615 : GEN
1616 78297209 : Flx_divrem_pre(GEN x, GEN T, ulong p, ulong pi, GEN *pr)
1617 : {
1618 : GEN B, y;
1619 : long dy, dx, d;
1620 78297209 : if (pr==ONLY_REM) return Flx_rem_pre(x, T, p, pi);
1621 61666344 : y = get_Flx_red(T, &B);
1622 61666344 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1623 61666344 : if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1624 61531059 : return Flx_divrem_basecase(x,y,p,pi,pr);
1625 : else
1626 : {
1627 135285 : pari_sp av = avma;
1628 135285 : GEN mg = B? B: Flx_invBarrett_pre(y, p, pi);
1629 135285 : GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pi,pr);
1630 135285 : if (!q1) return gc_NULL(av);
1631 135285 : if (!pr || pr==ONLY_DIVIDES) return gc_leaf(av, q1);
1632 126692 : return gc_all(av, 2, &q1, pr);
1633 : }
1634 : }
1635 : GEN
1636 29837359 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1637 29837359 : { return Flx_divrem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p), pr); }
1638 :
1639 : GEN
1640 958143912 : Flx_rem_pre(GEN x, GEN T, ulong p, ulong pi)
1641 : {
1642 958143912 : GEN B, y = get_Flx_red(T, &B);
1643 958143912 : long d = degpol(x) - degpol(y);
1644 958143912 : if (d < 0) return Flx_copy(x);
1645 816028703 : if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1646 815528939 : return Flx_rem_basecase(x,y,p, pi);
1647 : else
1648 : {
1649 499764 : pari_sp av=avma;
1650 499764 : GEN mg = B ? B: Flx_invBarrett_pre(y, p, pi);
1651 499764 : GEN r = Flx_divrem_Barrett(x, mg, y, p, pi, ONLY_REM);
1652 499764 : return gc_leaf(av, r);
1653 : }
1654 : }
1655 : GEN
1656 42360736 : Flx_rem(GEN x, GEN T, ulong p)
1657 42360736 : { return Flx_rem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1658 :
1659 : /* reduce T mod (X^n - 1, p). Shallow function */
1660 : GEN
1661 4920861 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1662 : {
1663 4920861 : long i, j, L = lg(T), l = n+2;
1664 : GEN S;
1665 4920861 : if (L <= l || n & ~LGBITS) return T;
1666 3583 : S = cgetg(l, t_VECSMALL);
1667 3583 : S[1] = T[1];
1668 15896 : for (i = 2; i < l; i++) S[i] = T[i];
1669 10029 : for (j = 2; i < L; i++) {
1670 6446 : S[j] = Fl_add(S[j], T[i], p);
1671 6446 : if (++j == l) j = 2;
1672 : }
1673 3583 : return Flx_renormalize(S, l);
1674 : }
1675 : /* reduce T mod (X^n + 1, p). Shallow function */
1676 : GEN
1677 31667 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
1678 : {
1679 31667 : long i, j, L = lg(T), l = n+2, s = -1;
1680 : GEN S;
1681 31667 : if (L <= l || n & ~LGBITS) return T;
1682 2794 : S = cgetg(l, t_VECSMALL);
1683 2794 : S[1] = T[1];
1684 12831 : for (i = 2; i < l; i++) S[i] = T[i];
1685 7541 : for (j = 2; i < L; i++) {
1686 4747 : S[j] = s==-1 ? Fl_sub(S[j], T[i], p): Fl_add(S[j], T[i], p);
1687 4747 : if (++j == l) { j = 2; s = -s; }
1688 : }
1689 2794 : return Flx_renormalize(S, l);
1690 : }
1691 :
1692 : struct _Flxq {
1693 : GEN aut, T;
1694 : ulong p, pi;
1695 : };
1696 : /* allow pi = 0 */
1697 : static void
1698 70189208 : set_Flxq_pre(struct _Flxq *D, GEN T, ulong p, ulong pi)
1699 : {
1700 70189208 : D->p = p;
1701 70189208 : D->pi = pi;
1702 70189208 : D->T = Flx_get_red_pre(T, p, pi);
1703 70189208 : }
1704 : static void
1705 68811 : set_Flxq(struct _Flxq *D, GEN T, ulong p)
1706 68811 : { set_Flxq_pre(D, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1707 :
1708 : struct _Flx {
1709 : long v;
1710 : ulong p, pi;
1711 : };
1712 :
1713 : static GEN
1714 0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1715 : {
1716 0 : struct _Flx *D = (struct _Flx*) E;
1717 0 : return Flx_divrem_pre(x, y, D->p, D->pi, r);
1718 : }
1719 : static GEN
1720 0 : _Flx_add(void * E, GEN x, GEN y) {
1721 0 : struct _Flx *D = (struct _Flx*) E;
1722 0 : return Flx_add(x, y, D->p);
1723 : }
1724 : static GEN
1725 0 : _Flx_sub(void * E, GEN x, GEN y) {
1726 0 : struct _Flx *D = (struct _Flx*) E;
1727 0 : return Flx_sub(x, y, D->p);
1728 : }
1729 : static GEN
1730 13263924 : _Flx_mul(void *E, GEN x, GEN y) {
1731 13263924 : struct _Flx *D = (struct _Flx*) E;
1732 13263924 : return Flx_mul_pre(x, y, D->p, D->pi);
1733 : }
1734 : static GEN
1735 0 : _Flx_sqr(void *E, GEN x) {
1736 0 : struct _Flx *D = (struct _Flx*) E;
1737 0 : return Flx_sqr_pre(x, D->p, D->pi);
1738 : }
1739 :
1740 : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1741 :
1742 : GEN
1743 0 : Flx_digits(GEN x, GEN T, ulong p)
1744 : {
1745 : struct _Flx D;
1746 0 : long d = get_Flx_degree(T), n = (lgpol(x)+d-1)/d;
1747 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1748 0 : return gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1749 : }
1750 :
1751 : GEN
1752 0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1753 : {
1754 : struct _Flx D;
1755 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1756 0 : return gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1757 : }
1758 :
1759 : long
1760 4584270 : Flx_val(GEN x)
1761 : {
1762 4584270 : long i, l=lg(x);
1763 4584270 : if (l==2) return LONG_MAX;
1764 4593060 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1765 4584270 : return i-2;
1766 : }
1767 : long
1768 31189832 : Flx_valrem(GEN x, GEN *Z)
1769 : {
1770 31189832 : long v, i, l=lg(x);
1771 : GEN y;
1772 31189832 : if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1773 33370936 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1774 31189832 : v = i-2;
1775 31189832 : if (v == 0) { *Z = x; return 0; }
1776 1025466 : l -= v;
1777 1025466 : y = cgetg(l, t_VECSMALL); y[1] = x[1];
1778 2638577 : for (i=2; i<l; i++) y[i] = x[i+v];
1779 1025466 : *Z = y; return v;
1780 : }
1781 :
1782 : GEN
1783 23061986 : Flx_deriv(GEN z, ulong p)
1784 : {
1785 23061986 : long i,l = lg(z)-1;
1786 : GEN x;
1787 23061986 : if (l < 2) l = 2;
1788 23061986 : x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1789 23061986 : if (HIGHWORD(l | p))
1790 63204159 : for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1791 : else
1792 88502149 : for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1793 23061986 : return Flx_renormalize(x,l);
1794 : }
1795 :
1796 : static GEN
1797 423966 : Flx_integXn(GEN x, long n, ulong p)
1798 : {
1799 423966 : long i, lx = lg(x);
1800 : GEN y;
1801 423966 : if (lx == 2) return Flx_copy(x);
1802 414152 : y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1803 2106983 : for (i=2; i<lx; i++)
1804 : {
1805 1692831 : ulong xi = uel(x,i);
1806 1692831 : if (xi == 0)
1807 13352 : uel(y,i) = 0;
1808 : else
1809 : {
1810 1679479 : ulong j = n+i-1;
1811 1679479 : ulong d = ugcd(j, xi);
1812 1679479 : if (d==1)
1813 1023392 : uel(y,i) = Fl_div(xi, j, p);
1814 : else
1815 656087 : uel(y,i) = Fl_div(xi/d, j/d, p);
1816 : }
1817 : }
1818 414152 : return Flx_renormalize(y, lx);;
1819 : }
1820 :
1821 : GEN
1822 0 : Flx_integ(GEN x, ulong p)
1823 : {
1824 0 : long i, lx = lg(x);
1825 : GEN y;
1826 0 : if (lx == 2) return Flx_copy(x);
1827 0 : y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1828 0 : uel(y,2) = 0;
1829 0 : for (i=3; i<=lx; i++)
1830 0 : uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1831 0 : return Flx_renormalize(y, lx+1);;
1832 : }
1833 :
1834 : /* assume p prime */
1835 : GEN
1836 660758 : Flx_diff1(GEN P, ulong p)
1837 : {
1838 660758 : return Flx_sub(Flx_translate1(P, p), P, p);
1839 : }
1840 :
1841 : GEN
1842 421892 : Flx_deflate(GEN x0, long d)
1843 : {
1844 : GEN z, y, x;
1845 421892 : long i,id, dy, dx = degpol(x0);
1846 421892 : if (d == 1 || dx <= 0) return Flx_copy(x0);
1847 358220 : dy = dx/d;
1848 358220 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1849 358220 : z = y + 2;
1850 358220 : x = x0+ 2;
1851 1164920 : for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1852 358220 : return y;
1853 : }
1854 :
1855 : GEN
1856 161027 : Flx_inflate(GEN x0, long d)
1857 : {
1858 161027 : long i, id, dy, dx = degpol(x0);
1859 161027 : GEN x = x0 + 2, z, y;
1860 161027 : if (dx <= 0) return Flx_copy(x0);
1861 159965 : dy = dx*d;
1862 159965 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1863 159965 : z = y + 2;
1864 8999283 : for (i=0; i<=dy; i++) z[i] = 0;
1865 4380153 : for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1866 159965 : return y;
1867 : }
1868 :
1869 : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1870 : GEN
1871 149183 : Flx_splitting(GEN p, long k)
1872 : {
1873 149183 : long n = degpol(p), v = p[1], m, i, j, l;
1874 : GEN r;
1875 :
1876 149183 : m = n/k;
1877 149183 : r = cgetg(k+1,t_VEC);
1878 707319 : for(i=1; i<=k; i++)
1879 : {
1880 558136 : gel(r,i) = cgetg(m+3, t_VECSMALL);
1881 558136 : mael(r,i,1) = v;
1882 : }
1883 4596150 : for (j=1, i=0, l=2; i<=n; i++)
1884 : {
1885 4446967 : mael(r,j,l) = p[2+i];
1886 4446967 : if (j==k) { j=1; l++; } else j++;
1887 : }
1888 707319 : for(i=1; i<=k; i++)
1889 558136 : gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1890 149183 : return r;
1891 : }
1892 :
1893 : /* ux + vy */
1894 : static GEN
1895 360140 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p, ulong pi)
1896 360140 : { return Flx_add(Flx_mul_pre(u,x, p,pi), Flx_mul_pre(v,y, p,pi), p); }
1897 :
1898 : static GEN
1899 26703 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p, ulong pi)
1900 : {
1901 26703 : GEN res = cgetg(3, t_COL);
1902 26703 : gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p, pi);
1903 26703 : gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p, pi);
1904 26703 : return res;
1905 : }
1906 :
1907 : #if 0
1908 : static GEN
1909 : FlxM_mul2_old(GEN M, GEN N, ulong p)
1910 : {
1911 : GEN res = cgetg(3, t_MAT);
1912 : gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1913 : gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1914 : return res;
1915 : }
1916 : #endif
1917 : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1918 : static GEN
1919 7297 : FlxM_mul2(GEN A, GEN B, ulong p, ulong pi)
1920 : {
1921 7297 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1922 7297 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1923 7297 : GEN M1 = Flx_mul_pre(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p, pi);
1924 7297 : GEN M2 = Flx_mul_pre(Flx_add(A21,A22, p), B11, p, pi);
1925 7297 : GEN M3 = Flx_mul_pre(A11, Flx_sub(B12,B22, p), p, pi);
1926 7297 : GEN M4 = Flx_mul_pre(A22, Flx_sub(B21,B11, p), p, pi);
1927 7297 : GEN M5 = Flx_mul_pre(Flx_add(A11,A12, p), B22, p, pi);
1928 7297 : GEN M6 = Flx_mul_pre(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p, pi);
1929 7297 : GEN M7 = Flx_mul_pre(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p, pi);
1930 7297 : GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1931 7297 : GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1932 7297 : retmkmat22(Flx_add(T1,T2, p), Flx_add(M3,M5, p),
1933 : Flx_add(M2,M4, p), Flx_add(T3,T4, p));
1934 : }
1935 :
1936 : /* Return [0,1;1,-q]*M */
1937 : static GEN
1938 7125 : Flx_FlxM_qmul(GEN q, GEN M, ulong p, ulong pi)
1939 : {
1940 7125 : GEN u = Flx_mul_pre(gcoeff(M,2,1), q, p, pi);
1941 7125 : GEN v = Flx_mul_pre(gcoeff(M,2,2), q, p, pi);
1942 7125 : retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
1943 : Flx_sub(gcoeff(M,1,1), u, p), Flx_sub(gcoeff(M,1,2), v, p));
1944 : }
1945 :
1946 : static GEN
1947 955 : matid2_FlxM(long v)
1948 955 : { retmkmat22(pol1_Flx(v),pol0_Flx(v),pol0_Flx(v),pol1_Flx(v)); }
1949 :
1950 : static GEN
1951 13 : matJ2_FlxM(long v)
1952 13 : { retmkmat22(pol0_Flx(v),pol1_Flx(v),pol1_Flx(v),pol0_Flx(v)); }
1953 :
1954 : struct Flx_res
1955 : {
1956 : ulong res, lc;
1957 : long deg0, deg1, off;
1958 : };
1959 :
1960 : INLINE void
1961 9405 : Flx_halfres_update_pre(long da, long db, long dr, ulong p, ulong pi, struct Flx_res *res)
1962 : {
1963 9405 : if (dr >= 0)
1964 : {
1965 9405 : if (res->lc != 1)
1966 : {
1967 7596 : if (pi)
1968 : {
1969 3127 : res->lc = Fl_powu_pre(res->lc, da - dr, p, pi);
1970 3127 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1971 : } else
1972 : {
1973 4469 : res->lc = Fl_powu(res->lc, da - dr, p);
1974 4469 : res->res = Fl_mul(res->res, res->lc, p);
1975 : }
1976 : }
1977 9405 : if (both_odd(da + res->off, db + res->off))
1978 63 : res->res = Fl_neg(res->res, p);
1979 : } else
1980 : {
1981 0 : if (db == 0)
1982 : {
1983 0 : if (res->lc != 1)
1984 : {
1985 0 : if (pi)
1986 : {
1987 0 : res->lc = Fl_powu_pre(res->lc, da, p, pi);
1988 0 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1989 : } else
1990 : {
1991 0 : res->lc = Fl_powu(res->lc, da, p);
1992 0 : res->res = Fl_mul(res->res, res->lc, p);
1993 : }
1994 : }
1995 : } else
1996 0 : res->res = 0;
1997 : }
1998 9405 : }
1999 :
2000 : static GEN
2001 1099347 : Flx_halfres_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *pa, GEN *pb, struct Flx_res *res)
2002 : {
2003 1099347 : pari_sp av = avma;
2004 : GEN u, u1, v, v1, M;
2005 1099347 : long vx = a[1], n = lgpol(a)>>1;
2006 1099347 : u1 = v = pol0_Flx(vx);
2007 1099347 : u = v1 = pol1_Flx(vx);
2008 6393020 : while (lgpol(b)>n)
2009 : {
2010 : GEN r, q;
2011 5293673 : q = Flx_divrem_pre(a,b,p,pi, &r);
2012 5293673 : if (res)
2013 : {
2014 8362 : long da = degpol(a), db=degpol(b), dr = degpol(r);
2015 8362 : res->lc = b[db+2];
2016 8362 : if (dr >= n)
2017 7133 : Flx_halfres_update_pre(da, db, dr, p, pi, res);
2018 : else
2019 : {
2020 1229 : res->deg0 = da;
2021 1229 : res->deg1 = db;
2022 : }
2023 : }
2024 5293673 : a = b; b = r; swap(u,u1); swap(v,v1);
2025 5293673 : u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
2026 5293673 : v1 = Flx_sub(v1, Flx_mul(v, q, p), p);
2027 5293673 : if (gc_needed(av,2))
2028 : {
2029 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
2030 0 : (void)gc_all(av,6, &a,&b,&u1,&v1,&u,&v);
2031 : }
2032 : }
2033 1099347 : M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
2034 1099347 : return gc_all(av,3, &M, pa, pb);
2035 : }
2036 :
2037 : static GEN Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res);
2038 :
2039 : static GEN
2040 20513 : Flx_halfres_split(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2041 : {
2042 20513 : pari_sp av = avma;
2043 : GEN R, S, T, V1, V2;
2044 : GEN x1, y1, r, q;
2045 20513 : long l = lgpol(x), n = l>>1, k;
2046 20513 : if (lgpol(y) <= n)
2047 915 : { *a = Flx_copy(x); *b = Flx_copy(y); return matid2_FlxM(x[1]); }
2048 19598 : if (res)
2049 : {
2050 3263 : res->lc = Flx_lead(y);
2051 3263 : res->deg0 -= n;
2052 3263 : res->deg1 -= n;
2053 3263 : res->off += n;
2054 : }
2055 19598 : R = Flx_halfres_i(Flx_shift(x,-n),Flx_shift(y,-n),p,pi,a,b,res);
2056 19598 : if (res)
2057 : {
2058 3263 : res->off -= n;
2059 3263 : res->deg0 += n;
2060 3263 : res->deg1 += n;
2061 : }
2062 19598 : V1 = FlxM_Flx_mul2(R, Flxn_red(x,n), Flxn_red(y,n), p, pi);
2063 19598 : x1 = Flx_add(Flx_shift(*a,n), gel(V1,1), p);
2064 19598 : y1 = Flx_add(Flx_shift(*b,n), gel(V1,2), p);
2065 19598 : if (lgpol(y1) <= n)
2066 12493 : { *a = x1; *b = y1; return gc_all(av, 3, &R, a, b); }
2067 7105 : k = 2*n-degpol(y1);
2068 7105 : q = Flx_divrem_pre(x1, y1, p, pi, &r);
2069 7105 : if (res)
2070 : {
2071 1043 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
2072 1043 : if (dy1 < degpol(y))
2073 185 : Flx_halfres_update_pre(res->deg0, res->deg1, dy1, p, pi, res);
2074 1043 : res->lc = uel(y1, dy1+2);
2075 1043 : res->deg0 = dx1;
2076 1043 : res->deg1 = dy1;
2077 1043 : if (dr >= n)
2078 : {
2079 1043 : Flx_halfres_update_pre(dx1, dy1, dr, p, pi, res);
2080 1043 : res->deg0 = dy1;
2081 1043 : res->deg1 = dr;
2082 : }
2083 1043 : res->deg0 -= k;
2084 1043 : res->deg1 -= k;
2085 1043 : res->off += k;
2086 : }
2087 7105 : S = Flx_halfres_i(Flx_shift(y1,-k), Flx_shift(r,-k), p, pi, a, b, res);
2088 7105 : if (res)
2089 : {
2090 1043 : res->deg0 += k;
2091 1043 : res->deg1 += k;
2092 1043 : res->off -= k;
2093 : }
2094 7105 : T = FlxM_mul2(S, Flx_FlxM_qmul(q, R, p,pi), p, pi);
2095 7105 : V2 = FlxM_Flx_mul2(S, Flxn_red(y1,k), Flxn_red(r,k), p, pi);
2096 7105 : *a = Flx_add(Flx_shift(*a,k), gel(V2,1), p);
2097 7105 : *b = Flx_add(Flx_shift(*b,k), gel(V2,2), p);
2098 7105 : return gc_all(av, 3, &T, a, b);
2099 : }
2100 :
2101 : static GEN
2102 1119860 : Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2103 : {
2104 1119860 : if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
2105 1099347 : return Flx_halfres_basecase(x, y, p, pi, a, b, res);
2106 20513 : return Flx_halfres_split(x, y, p, pi, a, b, res);
2107 : }
2108 :
2109 : static GEN
2110 1092113 : Flx_halfgcd_all_i(GEN x, GEN y, ulong p, ulong pi, GEN *pa, GEN *pb)
2111 : {
2112 : GEN a, b, R;
2113 1092113 : R = Flx_halfres_i(x, y, p, pi, &a, &b, NULL);
2114 1092113 : if (pa) *pa = a;
2115 1092113 : if (pb) *pb = b;
2116 1092113 : return R;
2117 : }
2118 :
2119 : /* Return M in GL_2(Fl[X]) such that:
2120 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
2121 : */
2122 :
2123 : GEN
2124 1092113 : Flx_halfgcd_all_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b)
2125 : {
2126 : pari_sp av;
2127 : GEN R, q, r;
2128 1092113 : long lx = lgpol(x), ly = lgpol(y);
2129 1092113 : if (!lx)
2130 : {
2131 0 : if (a) *a = Flx_copy(y);
2132 0 : if (b) *b = Flx_copy(x);
2133 0 : return matJ2_FlxM(x[1]);
2134 : }
2135 1092113 : if (ly < lx) return Flx_halfgcd_all_i(x, y, p, pi, a, b);
2136 8512 : av = avma;
2137 8512 : q = Flx_divrem(y,x,p,&r);
2138 8512 : R = Flx_halfgcd_all_i(x, r, p, pi, a, b);
2139 8512 : gcoeff(R,1,1) = Flx_sub(gcoeff(R,1,1), Flx_mul_pre(q,gcoeff(R,1,2), p,pi), p);
2140 8512 : gcoeff(R,2,1) = Flx_sub(gcoeff(R,2,1), Flx_mul_pre(q,gcoeff(R,2,2), p,pi), p);
2141 8512 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
2142 : }
2143 :
2144 : GEN
2145 154 : Flx_halfgcd_all(GEN x, GEN y, ulong p, GEN *a, GEN *b)
2146 154 : { return Flx_halfgcd_all_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), a, b); }
2147 :
2148 : GEN
2149 879237 : Flx_halfgcd_pre(GEN x, GEN y, ulong p, ulong pi)
2150 879237 : { return Flx_halfgcd_all_pre(x, y, p, pi, NULL, NULL); }
2151 :
2152 : GEN
2153 0 : Flx_halfgcd(GEN x, GEN y, ulong p)
2154 0 : { return Flx_halfgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2155 :
2156 : /*Do not garbage collect*/
2157 : static GEN
2158 87735564 : Flx_gcd_basecase(GEN a, GEN b, ulong p, ulong pi)
2159 : {
2160 87735564 : pari_sp av = avma;
2161 87735564 : ulong iter = 0;
2162 87735564 : if (lg(b) > lg(a)) swap(a, b);
2163 301756128 : while (lgpol(b))
2164 : {
2165 214020570 : GEN c = Flx_rem_pre(a,b,p,pi);
2166 214020564 : iter++; a = b; b = c;
2167 214020564 : if (gc_needed(av,2))
2168 : {
2169 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
2170 0 : (void)gc_all(av,2, &a,&b);
2171 : }
2172 : }
2173 87735558 : return iter < 2 ? Flx_copy(a) : a;
2174 : }
2175 :
2176 : GEN
2177 89494437 : Flx_gcd_pre(GEN x, GEN y, ulong p, ulong pi)
2178 : {
2179 89494437 : pari_sp av = avma;
2180 : long lim;
2181 89494437 : if (!lgpol(x)) return Flx_copy(y);
2182 87735564 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2183 87735885 : while (lgpol(y) >= lim)
2184 : {
2185 321 : if (lgpol(y)<=(lgpol(x)>>1))
2186 : {
2187 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2188 0 : x = y; y = r;
2189 : }
2190 321 : (void) Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2191 321 : if (gc_needed(av,2))
2192 : {
2193 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
2194 0 : (void)gc_all(av,2,&x,&y);
2195 : }
2196 : }
2197 87735564 : return gc_leaf(av, Flx_gcd_basecase(x,y,p,pi));
2198 : }
2199 : GEN
2200 36623665 : Flx_gcd(GEN x, GEN y, ulong p)
2201 36623665 : { return Flx_gcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2202 :
2203 : int
2204 9229580 : Flx_is_squarefree(GEN z, ulong p)
2205 : {
2206 9229580 : pari_sp av = avma;
2207 9229580 : GEN d = Flx_gcd(z, Flx_deriv(z,p), p);
2208 9229580 : return gc_bool(av, degpol(d) == 0);
2209 : }
2210 :
2211 : static long
2212 127014 : Flx_is_smooth_squarefree(GEN f, long r, ulong p, ulong pi)
2213 : {
2214 127014 : pari_sp av = avma;
2215 : long i;
2216 127014 : GEN sx = polx_Flx(f[1]), a = sx;
2217 127014 : for(i=1;;i++)
2218 : {
2219 536410 : if (degpol(f)<=r) return gc_long(av,1);
2220 515617 : a = Flxq_powu_pre(Flx_rem_pre(a,f,p,pi), p, f, p, pi);
2221 515617 : if (Flx_equal(a, sx)) return gc_long(av,1);
2222 511155 : if (i==r) return gc_long(av,0);
2223 409396 : f = Flx_div_pre(f, Flx_gcd_pre(Flx_sub(a,sx,p),f,p,pi),p,pi);
2224 : }
2225 : }
2226 :
2227 : static long
2228 8208 : Flx_is_l_pow(GEN x, ulong p)
2229 : {
2230 8208 : ulong i, lx = lgpol(x);
2231 16396 : for (i=1; i<lx; i++)
2232 14711 : if (x[i+2] && i%p) return 0;
2233 1685 : return 1;
2234 : }
2235 :
2236 : int
2237 127014 : Flx_is_smooth_pre(GEN g, long r, ulong p, ulong pi)
2238 : {
2239 : while (1)
2240 8208 : {
2241 127014 : GEN f = Flx_gcd_pre(g, Flx_deriv(g, p), p, pi);
2242 127014 : if (!Flx_is_smooth_squarefree(Flx_div_pre(g, f, p, pi), r, p, pi))
2243 101759 : return 0;
2244 25255 : if (degpol(f)==0) return 1;
2245 8208 : g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2246 : }
2247 : }
2248 : int
2249 74256 : Flx_is_smooth(GEN g, long r, ulong p)
2250 74256 : { return Flx_is_smooth_pre(g, r, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2251 :
2252 : static GEN
2253 6255875 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2254 : {
2255 6255875 : pari_sp av=avma;
2256 : GEN u,v,u1,v1;
2257 6255875 : long vx = a[1];
2258 6255875 : v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2259 6255875 : if (ptu) { u = pol1_Flx(vx); u1 = pol0_Flx(vx); }
2260 27596328 : while (lgpol(b))
2261 : {
2262 21340453 : GEN r, q = Flx_divrem_pre(a,b,p,pi, &r);
2263 21340453 : a = b; b = r;
2264 21340453 : if (ptu)
2265 : {
2266 2201089 : swap(u,u1);
2267 2201089 : u1 = Flx_sub(u1, Flx_mul_pre(u, q, p, pi), p);
2268 : }
2269 21340453 : swap(v,v1);
2270 21340453 : v1 = Flx_sub(v1, Flx_mul_pre(v, q, p, pi), p);
2271 21340453 : if (gc_needed(av,2))
2272 : {
2273 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(a));
2274 0 : (void)gc_all(av,ptu ? 6: 4, &a,&b,&v,&v1,&u,&u1);
2275 : }
2276 : }
2277 6255875 : if (ptu) *ptu = u;
2278 6255875 : *ptv = v;
2279 6255875 : return a;
2280 : }
2281 :
2282 : static GEN
2283 122234 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2284 : {
2285 : GEN u, v;
2286 122234 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2287 122234 : GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
2288 122234 : long i, n = 0, vs = x[1];
2289 333319 : while (lgpol(y) >= lim)
2290 : {
2291 211085 : if (lgpol(y)<=(lgpol(x)>>1))
2292 : {
2293 26 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2294 26 : x = y; y = r;
2295 26 : gel(V,++n) = mkmat22(pol0_Flx(vs),pol1_Flx(vs),pol1_Flx(vs),Flx_neg(q,p));
2296 : } else
2297 211059 : gel(V,++n) = Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2298 : }
2299 122234 : y = Flx_extgcd_basecase(x,y,p,pi,&u,&v);
2300 211085 : for (i = n; i>1; i--)
2301 : {
2302 88851 : GEN R = gel(V,i);
2303 88851 : GEN u1 = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2304 88851 : GEN v1 = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2305 88851 : u = u1; v = v1;
2306 : }
2307 : {
2308 122234 : GEN R = gel(V,1);
2309 122234 : if (ptu)
2310 6692 : *ptu = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2311 122234 : *ptv = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2312 : }
2313 122234 : return y;
2314 : }
2315 :
2316 : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2317 : * ux + vy = gcd (mod p) */
2318 : GEN
2319 6255875 : Flx_extgcd_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2320 : {
2321 6255875 : pari_sp av = avma;
2322 : GEN d;
2323 6255875 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2324 6255875 : if (lgpol(y) >= lim)
2325 122234 : d = Flx_extgcd_halfgcd(x, y, p, pi, ptu, ptv);
2326 : else
2327 6133641 : d = Flx_extgcd_basecase(x, y, p, pi, ptu, ptv);
2328 6255875 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
2329 : }
2330 : GEN
2331 861529 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2332 861529 : { return Flx_extgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptu, ptv); }
2333 :
2334 : static GEN
2335 1044 : Flx_halfres_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, ulong *r)
2336 : {
2337 : struct Flx_res res;
2338 : GEN R;
2339 : long dB;
2340 :
2341 1044 : res.res = *r;
2342 1044 : res.lc = Flx_lead(y);
2343 1044 : res.deg0 = degpol(x);
2344 1044 : res.deg1 = degpol(y);
2345 1044 : res.off = 0;
2346 1044 : R = Flx_halfres_i(x, y, p, pi, a, b, &res);
2347 1044 : dB = degpol(*b);
2348 1044 : if (dB < degpol(y))
2349 1044 : Flx_halfres_update_pre(res.deg0, res.deg1, dB, p, pi, &res);
2350 1044 : *r = res.res;
2351 1044 : return R;
2352 : }
2353 :
2354 : static ulong
2355 14609815 : Flx_resultant_basecase_pre(GEN a, GEN b, ulong p, ulong pi)
2356 : {
2357 : pari_sp av;
2358 : long da,db,dc;
2359 14609815 : ulong lb, res = 1UL;
2360 : GEN c;
2361 :
2362 14609815 : da = degpol(a);
2363 14609815 : db = degpol(b);
2364 14609815 : if (db > da)
2365 : {
2366 0 : swapspec(a,b, da,db);
2367 0 : if (both_odd(da,db)) res = p-res;
2368 : }
2369 14609815 : else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2370 14609815 : av = avma;
2371 120050461 : while (db)
2372 : {
2373 105446786 : lb = b[db+2];
2374 105446786 : c = Flx_rem_pre(a,b, p,pi);
2375 105446786 : a = b; b = c; dc = degpol(c);
2376 105446786 : if (dc < 0) return gc_long(av,0);
2377 :
2378 105440646 : if (both_odd(da,db)) res = p - res;
2379 105440646 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, da - dc, p, pi), p);
2380 105440646 : if (gc_needed(av,2))
2381 : {
2382 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant (da = %ld)",da);
2383 0 : (void)gc_all(av,2, &a,&b);
2384 : }
2385 105440646 : da = db; /* = degpol(a) */
2386 105440646 : db = dc; /* = degpol(b) */
2387 : }
2388 14603675 : return gc_ulong(av, Fl_mul(res, Fl_powu_pre(b[2], da, p, pi), p));
2389 : }
2390 :
2391 : ulong
2392 14611257 : Flx_resultant_pre(GEN x, GEN y, ulong p, ulong pi)
2393 : {
2394 14611257 : pari_sp av = avma;
2395 : long lim;
2396 14611257 : ulong res = 1;
2397 14611257 : long dx = degpol(x), dy = degpol(y);
2398 14611257 : if (dx < 0 || dy < 0) return 0;
2399 14609815 : if (dx < dy)
2400 : {
2401 1045004 : swap(x,y);
2402 1045004 : if (both_odd(dx, dy))
2403 1906 : res = Fl_neg(res, p);
2404 : }
2405 14609815 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2406 14610667 : while (lgpol(y) >= lim)
2407 : {
2408 852 : if (lgpol(y)<=(lgpol(x)>>1))
2409 : {
2410 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2411 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2412 0 : ulong ly = y[dy+2];
2413 0 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2414 0 : if (both_odd(dx, dy))
2415 0 : res = Fl_neg(res, p);
2416 0 : x = y; y = r;
2417 : }
2418 852 : (void) Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2419 852 : if (gc_needed(av,2))
2420 : {
2421 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_res (y = %ld)",degpol(y));
2422 0 : (void)gc_all(av,2,&x,&y);
2423 : }
2424 : }
2425 14609815 : return gc_ulong(av, Fl_mul(res, Flx_resultant_basecase_pre(x, y, p, pi), p));
2426 : }
2427 :
2428 : ulong
2429 4738003 : Flx_resultant(GEN a, GEN b, ulong p)
2430 4738003 : { return Flx_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2431 :
2432 : /* If resultant is 0, *ptU and *ptV are not set */
2433 : static ulong
2434 53 : Flx_extresultant_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2435 : {
2436 53 : GEN z,q,u,v, x = a, y = b;
2437 53 : ulong lb, res = 1UL;
2438 53 : pari_sp av = avma;
2439 : long dx, dy, dz;
2440 53 : long vs = a[1];
2441 :
2442 53 : u = pol0_Flx(vs);
2443 53 : v = pol1_Flx(vs); /* v = 1 */
2444 53 : dx = degpol(x);
2445 53 : dy = degpol(y);
2446 764 : while (dy)
2447 : { /* b u = x (a), b v = y (a) */
2448 711 : lb = y[dy+2];
2449 711 : q = Flx_divrem_pre(x,y, p, pi, &z);
2450 711 : x = y; y = z; /* (x,y) = (y, x - q y) */
2451 711 : dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2452 711 : z = Flx_sub(u, Flx_mul_pre(q,v, p, pi), p);
2453 711 : u = v; v = z; /* (u,v) = (v, u - q v) */
2454 :
2455 711 : if (both_odd(dx,dy)) res = p - res;
2456 711 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, dx-dz, p, pi), p);
2457 711 : dx = dy; /* = degpol(x) */
2458 711 : dy = dz; /* = degpol(y) */
2459 : }
2460 53 : res = Fl_mul(res, Fl_powu_pre(y[2], dx, p, pi), p);
2461 53 : lb = Fl_mul(res, Fl_inv(y[2],p), p);
2462 53 : v = gc_leaf(av, Flx_Fl_mul_pre(v, lb, p, pi));
2463 53 : av = avma;
2464 53 : u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul_pre(b,v,p,pi), p);
2465 53 : u = gc_leaf(av, Flx_div_pre(u,a,p,pi)); /* = (res - b v) / a */
2466 53 : *ptU = u;
2467 53 : *ptV = v; return res;
2468 : }
2469 :
2470 : ulong
2471 53 : Flx_extresultant_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2472 : {
2473 53 : pari_sp av=avma;
2474 : GEN u, v, R;
2475 53 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2476 53 : ulong res = 1, res1;
2477 53 : long dx = degpol(x), dy = degpol(y);
2478 53 : if (dy > dx)
2479 : {
2480 13 : swap(x,y); lswap(dx,dy);
2481 13 : if (both_odd(dx,dy)) res = p-res;
2482 13 : R = matJ2_FlxM(x[1]);
2483 40 : } else R = matid2_FlxM(x[1]);
2484 53 : if (dy < 0) return 0;
2485 245 : while (lgpol(y) >= lim)
2486 : {
2487 : GEN M;
2488 192 : if (lgpol(y)<=(lgpol(x)>>1))
2489 : {
2490 20 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2491 20 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2492 20 : ulong ly = y[dy+2];
2493 20 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2494 20 : if (both_odd(dx, dy))
2495 0 : res = Fl_neg(res, p);
2496 20 : x = y; y = r;
2497 20 : R = Flx_FlxM_qmul(q, R, p,pi);
2498 : }
2499 192 : M = Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2500 192 : if (!res) return gc_ulong(av, 0);
2501 192 : R = FlxM_mul2(M, R, p, pi);
2502 192 : (void)gc_all(av,3,&x,&y,&R);
2503 : }
2504 53 : res1 = Flx_extresultant_basecase(x,y,p,pi,&u,&v);
2505 53 : if (!res1) return gc_ulong(av, 0);
2506 53 : *ptU = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi), res, p, pi);
2507 53 : *ptV = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi), res, p, pi);
2508 53 : (void)gc_all(av, 2, ptU, ptV);
2509 53 : return Fl_mul(res1,res,p);
2510 : }
2511 :
2512 : ulong
2513 53 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2514 53 : { return Flx_extresultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptU, ptV); }
2515 :
2516 : /* allow pi = 0 (SMALL_ULONG) */
2517 : ulong
2518 49647471 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2519 : {
2520 49647471 : ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2521 :
2522 49647471 : if (lx == 1) return 0;
2523 46599599 : x++;
2524 46599599 : if (pi)
2525 : {
2526 : LOCAL_OVERFLOW;
2527 : LOCAL_HIREMAINDER;
2528 46528689 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2529 117316938 : while (++i < lx)
2530 : {
2531 70788249 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2532 70788249 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2533 : }
2534 81325 : return v1? remlll_pre(v1, h1, l1, p, pi)
2535 46610014 : : remll_pre(h1, l1, p, pi);
2536 : }
2537 : else
2538 : {
2539 70910 : l1 = x[i] * y[i];
2540 30946658 : while (++i < lx) { l1 += x[i] * y[i]; if (l1 & HIGHBIT) l1 %= p; }
2541 70910 : return l1 % p;
2542 : }
2543 : }
2544 :
2545 : /* allow pi = 0 (SMALL_ULONG) */
2546 : ulong
2547 136438260 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2548 : {
2549 136438260 : long i, n = degpol(x);
2550 : ulong t;
2551 136438260 : if (n <= 0) return n? 0: x[2];
2552 39669967 : if (n > 15)
2553 : {
2554 180213 : pari_sp av = avma;
2555 180213 : GEN v = Fl_powers_pre(y, n, p, pi);
2556 180213 : return gc_ulong(av, Flx_eval_powers_pre(x, v, p, pi));
2557 : }
2558 39489754 : i = n+2; t = x[i];
2559 39489754 : if (pi)
2560 : {
2561 137110086 : for (i--; i>=2; i--) t = Fl_addmul_pre(uel(x, i), t, y, p, pi);
2562 38270447 : return t;
2563 : }
2564 2895297 : for (i--; i>=2; i--) t = (t * y + x[i]) % p;
2565 1219307 : return t %= p;
2566 : }
2567 : ulong
2568 20507460 : Flx_eval(GEN x, ulong y, ulong p)
2569 20507460 : { return Flx_eval_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2570 :
2571 : ulong
2572 3719 : Flv_prod_pre(GEN x, ulong p, ulong pi)
2573 : {
2574 3719 : pari_sp ltop = avma;
2575 : GEN v;
2576 3719 : long i,k,lx = lg(x);
2577 3719 : if (lx == 1) return 1UL;
2578 3719 : if (lx == 2) return uel(x,1);
2579 3152 : v = cgetg(1+(lx << 1), t_VECSMALL);
2580 3152 : k = 1;
2581 26952 : for (i=1; i<lx-1; i+=2)
2582 23800 : uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2583 3152 : if (i < lx) uel(v,k++) = uel(x,i);
2584 13043 : while (k > 2)
2585 : {
2586 9891 : lx = k; k = 1;
2587 33691 : for (i=1; i<lx-1; i+=2)
2588 23800 : uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2589 9891 : if (i < lx) uel(v,k++) = uel(v,i);
2590 : }
2591 3152 : return gc_ulong(ltop, uel(v,1));
2592 : }
2593 :
2594 : ulong
2595 0 : Flv_prod(GEN v, ulong p)
2596 : {
2597 0 : return Flv_prod_pre(v, p, get_Fl_red(p));
2598 : }
2599 :
2600 : GEN
2601 0 : FlxV_prod(GEN V, ulong p)
2602 : {
2603 : struct _Flx D;
2604 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2605 0 : return gen_product(V, (void *)&D, &_Flx_mul);
2606 : }
2607 :
2608 : static GEN
2609 0 : _Flx_pow(void* E, GEN x, GEN y)
2610 0 : { struct _Flx *D = (struct _Flx *)E; return Flx_powu(x, itou(y), D->p); }
2611 : static GEN
2612 0 : _Flx_one(void *E)
2613 0 : { struct _Flx *D = (struct _Flx *)E; return pol1_Flx(D->v); }
2614 :
2615 : GEN
2616 0 : FlxV_factorback(GEN f, GEN e, ulong p, long v)
2617 : {
2618 : struct _Flx D;
2619 0 : D.p = p; D.v = v;
2620 0 : D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2621 0 : return gen_factorback(f, e, (void *)&D, &_Flx_mul, &_Flx_pow, &_Flx_one);
2622 : }
2623 :
2624 : /* compute prod (x - a[i]) */
2625 : GEN
2626 904279 : Flv_roots_to_pol(GEN a, ulong p, long vs)
2627 : {
2628 : struct _Flx D;
2629 904279 : long i,k,lx = lg(a);
2630 : GEN p1;
2631 904279 : if (lx == 1) return pol1_Flx(vs);
2632 904279 : p1 = cgetg(lx, t_VEC);
2633 15013125 : for (k=1,i=1; i<lx-1; i+=2)
2634 14108846 : gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2635 14108846 : Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2636 904279 : if (i < lx)
2637 59357 : gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2638 904279 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2639 904279 : setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2640 : }
2641 :
2642 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2643 : INLINE void
2644 21753054 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2645 : {
2646 21753054 : pari_sp av = avma;
2647 21753054 : long n = lg(w), i;
2648 : ulong u;
2649 : GEN c;
2650 :
2651 21753054 : if (n == 1) return;
2652 21753054 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2653 93673980 : for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2654 21753054 : i = n-1; u = Fl_inv(c[i], p);
2655 93673980 : for ( ; i > 1; --i)
2656 : {
2657 71920926 : ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2658 71920926 : u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2659 : }
2660 21753054 : v[1] = u; set_avma(av);
2661 : }
2662 :
2663 : void
2664 20037772 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2665 :
2666 : GEN
2667 9646 : Flv_inv_pre(GEN w, ulong p, ulong pi)
2668 9646 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2669 :
2670 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2671 : INLINE void
2672 52847 : Flv_inv_indir(GEN w, GEN v, ulong p)
2673 : {
2674 52847 : pari_sp av = avma;
2675 52847 : long n = lg(w), i;
2676 : ulong u;
2677 : GEN c;
2678 :
2679 52847 : if (n == 1) return;
2680 52847 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2681 1827626 : for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2682 52847 : i = n-1; u = Fl_inv(c[i], p);
2683 1827626 : for ( ; i > 1; --i)
2684 : {
2685 1774779 : ulong t = Fl_mul(u, c[i-1], p);
2686 1774779 : u = Fl_mul(u, w[i], p); v[i] = t;
2687 : }
2688 52847 : v[1] = u; set_avma(av);
2689 : }
2690 : static void
2691 1758483 : Flv_inv_i(GEN v, GEN w, ulong p)
2692 : {
2693 1758483 : if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2694 1705636 : else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2695 1758483 : }
2696 : void
2697 12017 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2698 : GEN
2699 1746466 : Flv_inv(GEN w, ulong p)
2700 1746466 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2701 :
2702 : GEN
2703 40744193 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2704 : {
2705 40744193 : long l = lg(a), i;
2706 : GEN a0, z0, z;
2707 40744193 : if (l <= 3)
2708 : {
2709 0 : if (rem) *rem = l == 2? 0: a[2];
2710 0 : return zero_Flx(a[1]);
2711 : }
2712 40744193 : z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2713 40744193 : a0 = a + l-1;
2714 40744193 : z0 = z + l-2; *z0 = *a0--;
2715 40744193 : if (SMALL_ULONG(p))
2716 : {
2717 94815787 : for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2718 : {
2719 68944366 : ulong t = (*a0-- + x * *z0--) % p;
2720 68944366 : *z0 = (long)t;
2721 : }
2722 25871421 : if (rem) *rem = (*a0 + x * *z0) % p;
2723 : }
2724 : else
2725 : {
2726 56099130 : for (i=l-3; i>1; i--)
2727 : {
2728 41226358 : ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2729 41226358 : *z0 = (long)t;
2730 : }
2731 14872772 : if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2732 : }
2733 40744193 : return z;
2734 : }
2735 :
2736 : /* xa, ya = t_VECSMALL */
2737 : static GEN
2738 1747672 : Flv_producttree(GEN xa, GEN s, ulong p, ulong pi, long vs)
2739 : {
2740 1747672 : long n = lg(xa)-1;
2741 1747672 : long m = n==1 ? 1: expu(n-1)+1;
2742 1747672 : long i, j, k, ls = lg(s);
2743 1747672 : GEN T = cgetg(m+1, t_VEC);
2744 1747672 : GEN t = cgetg(ls, t_VEC);
2745 11664058 : for (j=1, k=1; j<ls; k+=s[j++])
2746 9916386 : gel(t, j) = s[j] == 1 ?
2747 9916386 : mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2748 3308194 : mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2749 3308194 : Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2750 1747672 : gel(T,1) = t;
2751 4780837 : for (i=2; i<=m; i++)
2752 : {
2753 3033165 : GEN u = gel(T, i-1);
2754 3033165 : long n = lg(u)-1;
2755 3033165 : GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2756 11201879 : for (j=1, k=1; k<n; j++, k+=2)
2757 8168714 : gel(t, j) = Flx_mul_pre(gel(u, k), gel(u, k+1), p, pi);
2758 3033165 : gel(T, i) = t;
2759 : }
2760 1747672 : return T;
2761 : }
2762 :
2763 : static GEN
2764 1788391 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p, ulong pi)
2765 : {
2766 : long i,j,k;
2767 1788391 : long m = lg(T)-1;
2768 1788391 : GEN R = cgetg(lg(xa), t_VECSMALL);
2769 1788391 : GEN Tp = cgetg(m+1, t_VEC), t;
2770 1788391 : gel(Tp, m) = mkvec(P);
2771 5009016 : for (i=m-1; i>=1; i--)
2772 : {
2773 3220625 : GEN u = gel(T, i), v = gel(Tp, i+1);
2774 3220625 : long n = lg(u)-1;
2775 3220625 : t = cgetg(n+1, t_VEC);
2776 12433928 : for (j=1, k=1; k<n; j++, k+=2)
2777 : {
2778 9213303 : gel(t, k) = Flx_rem_pre(gel(v, j), gel(u, k), p, pi);
2779 9213303 : gel(t, k+1) = Flx_rem_pre(gel(v, j), gel(u, k+1), p, pi);
2780 : }
2781 3220625 : gel(Tp, i) = t;
2782 : }
2783 : {
2784 1788391 : GEN u = gel(T, i+1), v = gel(Tp, i+1);
2785 1788391 : long n = lg(u)-1;
2786 12790085 : for (j=1, k=1; j<=n; j++)
2787 : {
2788 11001694 : long c, d = degpol(gel(u,j));
2789 25572570 : for (c=1; c<=d; c++, k++) R[k] = Flx_eval_pre(gel(v, j), xa[k], p, pi);
2790 : }
2791 1788391 : return gc_const((pari_sp)R, R);
2792 : }
2793 : }
2794 :
2795 : static GEN
2796 2648470 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, ulong pi, long vs)
2797 : {
2798 2648470 : pari_sp av = avma;
2799 2648470 : long m = lg(T)-1;
2800 2648470 : long i, j, k, ls = lg(s);
2801 2648470 : GEN Tp = cgetg(m+1, t_VEC);
2802 2648470 : GEN t = cgetg(ls, t_VEC);
2803 32778212 : for (j=1, k=1; j<ls; k+=s[j++])
2804 30129742 : if (s[j]==2)
2805 : {
2806 10539845 : ulong a = Fl_mul(ya[k], R[k], p);
2807 10539845 : ulong b = Fl_mul(ya[k+1], R[k+1], p);
2808 10539845 : gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2809 10539845 : Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2810 10539845 : gel(t, j) = Flx_renormalize(gel(t, j), 4);
2811 : }
2812 : else
2813 19589897 : gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2814 2648470 : gel(Tp, 1) = t;
2815 9633648 : for (i=2; i<=m; i++)
2816 : {
2817 6985178 : GEN u = gel(T, i-1);
2818 6985178 : GEN t = cgetg(lg(gel(T,i)), t_VEC);
2819 6985178 : GEN v = gel(Tp, i-1);
2820 6985178 : long n = lg(v)-1;
2821 34466450 : for (j=1, k=1; k<n; j++, k+=2)
2822 27481272 : gel(t, j) = Flx_add(Flx_mul_pre(gel(u, k), gel(v, k+1), p, pi),
2823 27481272 : Flx_mul_pre(gel(u, k+1), gel(v, k), p, pi), p);
2824 6985178 : gel(Tp, i) = t;
2825 : }
2826 2648470 : return gc_leaf(av, gmael(Tp,m,1));
2827 : }
2828 :
2829 : GEN
2830 0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2831 : {
2832 0 : pari_sp av = avma;
2833 0 : GEN s = producttree_scheme(lg(xa)-1);
2834 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2835 0 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2836 0 : return gc_leaf(av, Flx_Flv_multieval_tree(P, xa, T, p, pi));
2837 : }
2838 :
2839 : static GEN
2840 2478 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p, ulong pi)
2841 45675 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p, pi)) }
2842 :
2843 : GEN
2844 2478 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2845 : {
2846 2478 : pari_sp av = avma;
2847 2478 : GEN s = producttree_scheme(lg(xa)-1);
2848 2478 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2849 2478 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2850 2478 : return gc_upto(av, FlxV_Flv_multieval_tree(P, xa, T, p, pi));
2851 : }
2852 :
2853 : GEN
2854 1485852 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2855 : {
2856 1485852 : pari_sp av = avma;
2857 1485852 : GEN s = producttree_scheme(lg(xa)-1);
2858 1485852 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2859 1485852 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2860 1485852 : long m = lg(T)-1;
2861 1485852 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2862 1485852 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2863 1485852 : return gc_leaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, pi, vs));
2864 : }
2865 :
2866 : GEN
2867 105897 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2868 : {
2869 105897 : pari_sp av = avma;
2870 105897 : GEN s = producttree_scheme(lg(xa)-1);
2871 105897 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2872 105897 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2873 105897 : long i, m = lg(T)-1, l = lg(ya)-1;
2874 105897 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2875 105897 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2876 105897 : GEN M = cgetg(l+1, t_VEC);
2877 1268515 : for (i=1; i<=l; i++)
2878 1162618 : gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, pi, vs);
2879 105897 : return gc_upto(av, M);
2880 : }
2881 :
2882 : GEN
2883 153445 : Flv_invVandermonde(GEN L, ulong den, ulong p)
2884 : {
2885 153445 : pari_sp av = avma;
2886 153445 : long i, n = lg(L);
2887 : GEN M, R;
2888 153445 : GEN s = producttree_scheme(n-1);
2889 153445 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2890 153445 : GEN tree = Flv_producttree(L, s, p, pi, 0);
2891 153445 : long m = lg(tree)-1;
2892 153445 : GEN T = gmael(tree, m, 1);
2893 153445 : R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p, pi), p);
2894 153445 : if (den!=1) R = Flv_Fl_mul(R, den, p);
2895 153445 : M = cgetg(n, t_MAT);
2896 603585 : for (i = 1; i < n; i++)
2897 : {
2898 450140 : GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2899 450140 : gel(M,i) = Flx_to_Flv(P, n-1);
2900 : }
2901 153445 : return gc_GEN(av, M);
2902 : }
2903 :
2904 : /***********************************************************************/
2905 : /** Flxq **/
2906 : /***********************************************************************/
2907 : /* Flxq objects are Flx modulo another Flx called q. */
2908 :
2909 : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2910 : GEN
2911 194806362 : Flxq_mul_pre(GEN x,GEN y,GEN T,ulong p,ulong pi)
2912 194806362 : { return Flx_rem_pre(Flx_mul_pre(x,y,p,pi),T,p,pi); }
2913 : GEN
2914 18090061 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2915 18090061 : { return Flxq_mul_pre(x,y,T,p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2916 :
2917 : GEN
2918 281870995 : Flxq_sqr_pre(GEN x,GEN T,ulong p,ulong pi)
2919 281870995 : { return Flx_rem_pre(Flx_sqr_pre(x, p,pi), T, p,pi); }
2920 : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2921 : GEN
2922 2760620 : Flxq_sqr(GEN x,GEN T,ulong p)
2923 2760620 : { return Flxq_sqr_pre(x,T,p,SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2924 :
2925 : static GEN
2926 1549491 : _Flxq_red(void *E, GEN x)
2927 1549491 : { struct _Flxq *s = (struct _Flxq *)E;
2928 1549491 : return Flx_rem_pre(x, s->T, s->p, s->pi); }
2929 : static GEN
2930 389834 : _Flxq_add(void *E, GEN x, GEN y)
2931 389834 : { struct _Flxq *s = (struct _Flxq *)E;
2932 389834 : return Flx_add(x,y,s->p); }
2933 : static GEN
2934 0 : _Flxq_sub(void *E, GEN x, GEN y)
2935 0 : { struct _Flxq *s = (struct _Flxq *)E;
2936 0 : return Flx_sub(x,y,s->p); }
2937 : static GEN
2938 275780571 : _Flxq_sqr(void *data, GEN x)
2939 : {
2940 275780571 : struct _Flxq *D = (struct _Flxq*)data;
2941 275780571 : return Flxq_sqr_pre(x, D->T, D->p, D->pi);
2942 : }
2943 : static GEN
2944 149810385 : _Flxq_mul(void *data, GEN x, GEN y)
2945 : {
2946 149810385 : struct _Flxq *D = (struct _Flxq*)data;
2947 149810385 : return Flxq_mul_pre(x,y, D->T, D->p, D->pi);
2948 : }
2949 : static GEN
2950 22642779 : _Flxq_one(void *data)
2951 : {
2952 22642779 : struct _Flxq *D = (struct _Flxq*)data;
2953 22642779 : return pol1_Flx(get_Flx_var(D->T));
2954 : }
2955 : static GEN
2956 0 : _Flxq_zero(void *data)
2957 : {
2958 0 : struct _Flxq *D = (struct _Flxq*)data;
2959 0 : return pol0_Flx(get_Flx_var(D->T));
2960 : }
2961 : static GEN
2962 23260288 : _Flxq_powu_i(struct _Flxq *D, GEN x, ulong n)
2963 23260288 : { return gen_powu_i(x, n, (void*)D, &_Flxq_sqr, &_Flxq_mul); }
2964 : static GEN
2965 68 : _Flxq_powu(struct _Flxq *D, GEN x, ulong n)
2966 68 : { pari_sp av = avma; return gc_leaf(av, _Flxq_powu_i(D, x, n)); }
2967 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2968 : GEN
2969 24517634 : Flxq_powu_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
2970 : {
2971 : pari_sp av;
2972 : struct _Flxq D;
2973 24517634 : switch(n)
2974 : {
2975 0 : case 0: return pol1_Flx(get_Flx_var(T));
2976 279256 : case 1: return Flx_copy(x);
2977 978158 : case 2: return Flxq_sqr_pre(x, T, p, pi);
2978 : }
2979 23260220 : av = avma; set_Flxq_pre(&D, T, p, pi);
2980 23260220 : return gc_leaf(av, _Flxq_powu_i(&D, x, n));
2981 : }
2982 : GEN
2983 500641 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2984 500641 : { return Flxq_powu_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2985 :
2986 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2987 : GEN
2988 23678551 : Flxq_pow_pre(GEN x, GEN n, GEN T, ulong p, ulong pi)
2989 : {
2990 23678551 : pari_sp av = avma;
2991 : struct _Flxq D;
2992 : GEN y;
2993 23678551 : long s = signe(n);
2994 23678551 : if (!s) return pol1_Flx(get_Flx_var(T));
2995 23600839 : if (s < 0) x = Flxq_inv_pre(x,T,p,pi);
2996 23600839 : if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2997 23079305 : set_Flxq_pre(&D, T, p, pi);
2998 23079305 : y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2999 23079305 : return gc_leaf(av, y);
3000 : }
3001 : GEN
3002 934136 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
3003 934136 : { return Flxq_pow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3004 :
3005 : GEN
3006 28 : Flxq_pow_init_pre(GEN x, GEN n, long k, GEN T, ulong p, ulong pi)
3007 : {
3008 28 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3009 28 : return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
3010 : }
3011 : GEN
3012 0 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
3013 0 : { return Flxq_pow_init_pre(x, n, k, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3014 :
3015 : GEN
3016 4393 : Flxq_pow_table_pre(GEN R, GEN n, GEN T, ulong p, ulong pi)
3017 : {
3018 4393 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3019 4393 : return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
3020 : }
3021 : GEN
3022 0 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
3023 0 : { return Flxq_pow_table_pre(R, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3024 :
3025 : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
3026 : * not stack clean. */
3027 : GEN
3028 5394346 : Flxq_invsafe_pre(GEN x, GEN T, ulong p, ulong pi)
3029 : {
3030 5394346 : GEN V, z = Flx_extgcd_pre(get_Flx_mod(T), x, p, pi, NULL, &V);
3031 : ulong iz;
3032 5394346 : if (degpol(z)) return NULL;
3033 5393682 : iz = Fl_inv(uel(z,2), p);
3034 5393682 : return Flx_Fl_mul_pre(V, iz, p, pi);
3035 : }
3036 : GEN
3037 671082 : Flxq_invsafe(GEN x, GEN T, ulong p)
3038 671082 : { return Flxq_invsafe_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3039 :
3040 : GEN
3041 4240179 : Flxq_inv_pre(GEN x, GEN T, ulong p, ulong pi)
3042 : {
3043 4240179 : pari_sp av=avma;
3044 4240179 : GEN U = Flxq_invsafe_pre(x, T, p, pi);
3045 4240179 : if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
3046 4240151 : return gc_leaf(av, U);
3047 : }
3048 : GEN
3049 335914 : Flxq_inv(GEN x, GEN T, ulong p)
3050 335914 : { return Flxq_inv_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3051 :
3052 : GEN
3053 2323066 : Flxq_div_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
3054 : {
3055 2323066 : pari_sp av = avma;
3056 2323066 : return gc_leaf(av, Flxq_mul_pre(x,Flxq_inv_pre(y,T,p,pi),T,p,pi));
3057 : }
3058 : GEN
3059 1170208 : Flxq_div(GEN x, GEN y, GEN T, ulong p)
3060 1170208 : { return Flxq_div_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3061 :
3062 : GEN
3063 22638386 : Flxq_powers_pre(GEN x, long l, GEN T, ulong p, ulong pi)
3064 : {
3065 : struct _Flxq D;
3066 22638386 : long d = degpol(x), dT = get_Flx_degree(T);
3067 22638386 : if (d >= dT) { x = Flx_rem_pre(x, T, p, pi); d = degpol(x); }
3068 22638386 : set_Flxq_pre(&D, T, p, pi);
3069 22638386 : return gen_powers(x, l, 2*d>=dT, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
3070 : }
3071 : GEN
3072 232435 : Flxq_powers(GEN x, long l, GEN T, ulong p)
3073 232435 : { return Flxq_powers_pre(x, l, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3074 :
3075 : GEN
3076 171249 : Flxq_matrix_pow_pre(GEN y, long n, long m, GEN P, ulong l, ulong li)
3077 171249 : { return FlxV_to_Flm(Flxq_powers_pre(y,m-1,P,l,li),n); }
3078 : GEN
3079 399 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
3080 399 : { return Flxq_matrix_pow_pre(y, n, m, P, l, SMALL_ULONG(l)? 0: get_Fl_red(l)); }
3081 :
3082 : GEN
3083 14056189 : Flx_Frobenius_pre(GEN T, ulong p, ulong pi)
3084 14056189 : { return Flxq_powu_pre(polx_Flx(get_Flx_var(T)), p, T, p, pi); }
3085 : GEN
3086 87438 : Flx_Frobenius(GEN T, ulong p)
3087 87438 : { return Flx_Frobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3088 :
3089 : GEN
3090 86900 : Flx_matFrobenius_pre(GEN T, ulong p, ulong pi)
3091 : {
3092 86900 : long n = get_Flx_degree(T);
3093 86900 : return Flxq_matrix_pow_pre(Flx_Frobenius_pre(T, p, pi), n, n, T, p, pi);
3094 : }
3095 : GEN
3096 0 : Flx_matFrobenius(GEN T, ulong p)
3097 0 : { return Flx_matFrobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3098 :
3099 : static GEN
3100 13198433 : Flx_blocks_Flm(GEN P, long n, long m)
3101 : {
3102 13198433 : GEN z = cgetg(m+1,t_MAT);
3103 13198433 : long i,j, k=2, l = lg(P);
3104 37914241 : for(i=1; i<=m; i++)
3105 : {
3106 24715808 : GEN zi = cgetg(n+1,t_VECSMALL);
3107 24715808 : gel(z,i) = zi;
3108 114182971 : for(j=1; j<=n; j++)
3109 89467163 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3110 : }
3111 13198433 : return z;
3112 : }
3113 :
3114 : GEN
3115 529859 : Flx_blocks(GEN P, long n, long m)
3116 : {
3117 529859 : GEN z = cgetg(m+1,t_VEC);
3118 529859 : long i,j, k=2, l = lg(P);
3119 1589577 : for(i=1; i<=m; i++)
3120 : {
3121 1059718 : GEN zi = cgetg(n+2,t_VECSMALL);
3122 1059718 : zi[1] = P[1];
3123 1059718 : gel(z,i) = zi;
3124 7196204 : for(j=2; j<n+2; j++)
3125 6136486 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3126 1059718 : zi = Flx_renormalize(zi, n+2);
3127 : }
3128 529859 : return z;
3129 : }
3130 :
3131 : static GEN
3132 13198433 : FlxV_to_Flm_lg(GEN x, long m, long n)
3133 : {
3134 : long i;
3135 13198433 : GEN y = cgetg(n+1, t_MAT);
3136 62344001 : for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
3137 13198433 : return y;
3138 : }
3139 :
3140 : /* allow pi = 0 (SMALL_ULONG) */
3141 : GEN
3142 13454003 : Flx_FlxqV_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3143 : {
3144 13454003 : pari_sp btop, av = avma;
3145 13454003 : long sv = get_Flx_var(T), m = get_Flx_degree(T);
3146 13454003 : long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
3147 : GEN A, B, C, S, g;
3148 13454003 : if (lQ == 0) return pol0_Flx(sv);
3149 13198433 : if (lQ <= l)
3150 : {
3151 6538402 : n = l;
3152 6538402 : d = 1;
3153 : }
3154 : else
3155 : {
3156 6660031 : n = l-1;
3157 6660031 : d = (lQ+n-1)/n;
3158 : }
3159 13198433 : A = FlxV_to_Flm_lg(x, m, n);
3160 13198433 : B = Flx_blocks_Flm(Q, n, d);
3161 13198433 : C = gc_upto(av, Flm_mul(A, B, p));
3162 13198433 : g = gel(x, l);
3163 13198433 : if (pi && SMALL_ULONG(p)) pi = 0;
3164 13198433 : T = Flx_get_red_pre(T, p, pi);
3165 13198433 : btop = avma;
3166 13198433 : S = Flv_to_Flx(gel(C, d), sv);
3167 24715808 : for (i = d-1; i>0; i--)
3168 : {
3169 11517375 : S = Flx_add(Flxq_mul_pre(S, g, T, p, pi), Flv_to_Flx(gel(C,i), sv), p);
3170 11517375 : if (gc_needed(btop,1))
3171 0 : S = gc_leaf(btop, S);
3172 : }
3173 13198433 : return gc_leaf(av, S);
3174 : }
3175 : GEN
3176 5124 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
3177 5124 : { return Flx_FlxqV_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3178 :
3179 : /* allow pi = 0 (SMALL_ULONG) */
3180 : GEN
3181 2533419 : Flx_Flxq_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3182 : {
3183 2533419 : pari_sp av = avma;
3184 : GEN z, V;
3185 2533419 : long d = degpol(Q), rtd;
3186 2533419 : if (d < 0) return pol0_Flx(get_Flx_var(T));
3187 2533328 : rtd = (long) sqrt((double)d);
3188 2533328 : T = Flx_get_red_pre(T, p, pi);
3189 2533328 : V = Flxq_powers_pre(x, rtd, T, p, pi);
3190 2533328 : z = Flx_FlxqV_eval_pre(Q, V, T, p, pi);
3191 2533328 : return gc_upto(av, z);
3192 : }
3193 : GEN
3194 804900 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
3195 804900 : { return Flx_Flxq_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3196 :
3197 : /* allow pi = 0 (SMALL_ULONG) */
3198 : GEN
3199 0 : FlxC_FlxqV_eval_pre(GEN x, GEN v, GEN T, ulong p, ulong pi)
3200 0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval_pre(gel(x,i), v, T, p, pi)) }
3201 : GEN
3202 0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
3203 0 : { return FlxC_FlxqV_eval_pre(x, v, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3204 :
3205 : /* allow pi = 0 (SMALL_ULONG) */
3206 : GEN
3207 0 : FlxC_Flxq_eval_pre(GEN x, GEN F, GEN T, ulong p, ulong pi)
3208 : {
3209 0 : long d = brent_kung_optpow(get_Flx_degree(T)-1,lg(x)-1,1);
3210 0 : GEN Fp = Flxq_powers_pre(F, d, T, p, pi);
3211 0 : return FlxC_FlxqV_eval_pre(x, Fp, T, p, pi);
3212 : }
3213 : GEN
3214 0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
3215 0 : { return FlxC_Flxq_eval_pre(x, F, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3216 :
3217 : struct bb_algebra Flxq_algebra = { _Flxq_red, _Flxq_add, _Flxq_sub,
3218 : _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
3219 :
3220 : const struct bb_algebra *
3221 0 : get_Flxq_algebra(void **E, GEN T, ulong p)
3222 : {
3223 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3224 0 : GEN z = new_chunk(sizeof(struct _Flxq));
3225 0 : struct _Flxq *e = (struct _Flxq *) z;
3226 0 : e->T = Flx_get_red(T, p);
3227 0 : e->p = p;
3228 0 : e->pi = pi; *E = (void*)e;
3229 0 : return &Flxq_algebra;
3230 : }
3231 :
3232 : static GEN
3233 0 : _Flx_red(void *E, GEN x)
3234 0 : { (void) E; return x; }
3235 : static GEN
3236 0 : _Flx_zero(void *E)
3237 0 : { struct _Flx *D = (struct _Flx *)E; return pol0_Flx(D->v); }
3238 :
3239 : static struct bb_algebra Flx_algebra = { _Flx_red, _Flx_add, _Flx_sub,
3240 : _Flx_mul, _Flx_sqr, _Flx_one, _Flx_zero };
3241 :
3242 : const struct bb_algebra *
3243 0 : get_Flx_algebra(void **E, ulong p, long v)
3244 : {
3245 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3246 0 : GEN z = new_chunk(sizeof(struct _Flx));
3247 0 : struct _Flx *e = (struct _Flx *) z;
3248 0 : e->p = p; e->pi = pi; e->v = v;
3249 0 : *E = (void*)e; return &Flx_algebra;
3250 : }
3251 :
3252 : static GEN
3253 77454 : Flxq_autpow_sqr(void *E, GEN x)
3254 : {
3255 77454 : struct _Flxq *D = (struct _Flxq*)E;
3256 77454 : return Flx_Flxq_eval_pre(x, x, D->T, D->p, D->pi);
3257 : }
3258 : static GEN
3259 46259 : Flxq_autpow_msqr(void *E, GEN x)
3260 : {
3261 46259 : struct _Flxq *D = (struct _Flxq*)E;
3262 46259 : return Flx_FlxqV_eval_pre(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p, D->pi);
3263 : }
3264 :
3265 : GEN
3266 98538 : Flxq_autpow_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3267 : {
3268 98538 : pari_sp av = avma;
3269 : struct _Flxq D;
3270 : long d;
3271 98538 : if (n==0) return Flx_rem_pre(polx_Flx(x[1]), T, p, pi);
3272 98531 : if (n==1) return Flx_rem_pre(x, T, p, pi);
3273 61508 : set_Flxq_pre(&D, T, p, pi);
3274 61508 : d = brent_kung_optpow(get_Flx_degree(T), hammingu(n)-1, 1);
3275 61508 : D.aut = Flxq_powers_pre(x, d, T, p, D.pi);
3276 61508 : x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
3277 61508 : return gc_GEN(av, x);
3278 : }
3279 : GEN
3280 7 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
3281 7 : { return Flxq_autpow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3282 :
3283 : GEN
3284 1667 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
3285 : {
3286 1667 : long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
3287 : ulong i, pi;
3288 1667 : pari_sp av = avma;
3289 1667 : GEN xp, V = cgetg(l+2,t_VEC);
3290 1667 : gel(V,1) = polx_Flx(vT); if (l==0) return V;
3291 1667 : gel(V,2) = gcopy(x); if (l==1) return V;
3292 1667 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3293 1667 : T = Flx_get_red_pre(T, p, pi);
3294 1667 : d = brent_kung_optpow(dT-1, l-1, 1);
3295 1667 : xp = Flxq_powers_pre(x, d, T, p, pi);
3296 6998 : for(i = 3; i < l+2; i++)
3297 5331 : gel(V,i) = Flx_FlxqV_eval_pre(gel(V,i-1), xp, T, p, pi);
3298 1667 : return gc_GEN(av, V);
3299 : }
3300 :
3301 : static GEN
3302 112033 : Flxq_autsum_mul(void *E, GEN x, GEN y)
3303 : {
3304 112033 : struct _Flxq *D = (struct _Flxq*)E;
3305 112033 : GEN T = D->T;
3306 112033 : ulong p = D->p, pi = D->pi;
3307 112033 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3308 112033 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3309 112033 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3310 112033 : GEN V2 = Flxq_powers_pre(phi2, d, T, p, pi);
3311 112033 : GEN phi3 = Flx_FlxqV_eval_pre(phi1, V2, T, p, pi);
3312 112033 : GEN aphi = Flx_FlxqV_eval_pre(a1, V2, T, p, pi);
3313 112033 : GEN a3 = Flxq_mul_pre(aphi, a2, T, p, pi);
3314 112033 : return mkvec2(phi3, a3);
3315 : }
3316 : static GEN
3317 104860 : Flxq_autsum_sqr(void *E, GEN x)
3318 104860 : { return Flxq_autsum_mul(E, x, x); }
3319 :
3320 : static GEN
3321 98539 : Flxq_autsum_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3322 : {
3323 98539 : pari_sp av = avma;
3324 98539 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3325 98539 : x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
3326 98539 : return gc_GEN(av, x);
3327 : }
3328 : GEN
3329 0 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
3330 0 : { return Flxq_autsum_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3331 :
3332 : static GEN
3333 799629 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
3334 : {
3335 799629 : struct _Flxq *D = (struct _Flxq*)E;
3336 799629 : GEN T = D->T;
3337 799629 : ulong p = D->p, pi = D->pi;
3338 799629 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3339 799629 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3340 799629 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3341 799629 : GEN V1 = Flxq_powers_pre(phi1, d, T, p, pi);
3342 799629 : GEN phi3 = Flx_FlxqV_eval_pre(phi2, V1, T, p, pi);
3343 799629 : GEN aphi = Flx_FlxqV_eval_pre(a2, V1, T, p, pi);
3344 799629 : GEN a3 = Flx_add(a1, aphi, p);
3345 799629 : return mkvec2(phi3, a3);
3346 : }
3347 :
3348 : static GEN
3349 668191 : Flxq_auttrace_sqr(void *E, GEN x)
3350 668191 : { return Flxq_auttrace_mul(E, x, x); }
3351 :
3352 : GEN
3353 978018 : Flxq_auttrace_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3354 : {
3355 978018 : pari_sp av = avma;
3356 : struct _Flxq D;
3357 978018 : set_Flxq_pre(&D, T, p, pi);
3358 978018 : x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
3359 978018 : return gc_GEN(av, x);
3360 : }
3361 : GEN
3362 0 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
3363 0 : { return Flxq_auttrace_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3364 :
3365 : static long
3366 393685 : bounded_order(ulong p, GEN b, long k)
3367 : {
3368 393685 : GEN a = modii(utoipos(p), b);
3369 : long i;
3370 805183 : for(i = 1; i < k; i++)
3371 : {
3372 510658 : if (equali1(a)) return i;
3373 411498 : a = modii(muliu(a,p),b);
3374 : }
3375 294525 : return 0;
3376 : }
3377 :
3378 : /* n = (p^d-a)\b
3379 : * b = bb*p^vb
3380 : * p^k = 1 [bb]
3381 : * d = m*k+r+vb
3382 : * u = (p^k-1)/bb;
3383 : * v = (p^(r+vb)-a)/b;
3384 : * w = (p^(m*k)-1)/(p^k-1)
3385 : * n = p^r*w*u+v
3386 : * w*u = p^vb*(p^(m*k)-1)/b
3387 : * n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b */
3388 : static GEN
3389 22629129 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p, ulong pi)
3390 : {
3391 22629129 : pari_sp av=avma;
3392 22629129 : long d = get_Flx_degree(T);
3393 22629129 : GEN an = absi_shallow(n), z, q;
3394 22629129 : if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow_pre(x, n, T, p, pi);
3395 394020 : q = powuu(p, d);
3396 394020 : if (dvdii(q, n))
3397 : {
3398 315 : long vn = logint(an, utoipos(p));
3399 315 : GEN autvn = vn==1 ? aut: Flxq_autpow_pre(aut,vn,T,p,pi);
3400 315 : z = Flx_Flxq_eval_pre(x,autvn,T,p,pi);
3401 : } else
3402 : {
3403 393705 : GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
3404 : GEN bb, u, v, autk;
3405 393705 : long vb = Z_lvalrem(b,p,&bb);
3406 393705 : long m, r, k = is_pm1(bb)? 1: bounded_order(p,bb,d);
3407 393705 : if (!k || d-vb < k) return Flxq_pow_pre(x,n, T,p,pi);
3408 99173 : m = (d-vb)/k; r = (d-vb)%k;
3409 99173 : u = diviiexact(subiu(powuu(p,k),1),bb);
3410 99173 : v = diviiexact(subii(powuu(p,r+vb),a),b);
3411 99173 : autk = k==1 ? aut: Flxq_autpow_pre(aut,k,T,p,pi);
3412 99173 : if (r)
3413 : {
3414 448 : GEN autr = r==1 ? aut: Flxq_autpow_pre(aut,r,T,p,pi);
3415 448 : z = Flx_Flxq_eval_pre(x,autr,T,p,pi);
3416 98725 : } else z = x;
3417 99173 : if (m > 1) z = gel(Flxq_autsum_pre(mkvec2(autk, z), m, T, p, pi), 2);
3418 99173 : if (!is_pm1(u)) z = Flxq_pow_pre(z, u, T, p, pi);
3419 99173 : if (signe(v)) z = Flxq_mul_pre(z, Flxq_pow_pre(x, v, T, p, pi), T, p, pi);
3420 : }
3421 99488 : return gc_upto(av,signe(n)>0 ? z : Flxq_inv_pre(z,T,p,pi));
3422 : }
3423 :
3424 : static GEN
3425 22621711 : _Flxq_pow(void *data, GEN x, GEN n)
3426 : {
3427 22621711 : struct _Flxq *D = (struct _Flxq*)data;
3428 22621711 : return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p, D->pi);
3429 : }
3430 :
3431 : static GEN
3432 6291 : _Flxq_rand(void *data)
3433 : {
3434 6291 : pari_sp av=avma;
3435 6291 : struct _Flxq *D = (struct _Flxq*)data;
3436 : GEN z;
3437 : do
3438 : {
3439 6313 : set_avma(av);
3440 6313 : z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3441 6313 : } while (lgpol(z)==0);
3442 6291 : return z;
3443 : }
3444 :
3445 : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3446 : static GEN
3447 35414 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3448 : {
3449 35414 : pari_sp av = avma;
3450 : GEN q,n_q,ord,ordp, op;
3451 :
3452 35414 : if (a == 1UL) return gen_0;
3453 : /* p > 2 */
3454 :
3455 35414 : ordp = utoi(p - 1);
3456 35414 : ord = get_arith_Z(o);
3457 35414 : if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3458 35414 : if (a == p - 1) /* -1 */
3459 7746 : return gc_INT(av, shifti(ord,-1));
3460 27668 : ordp = gcdii(ordp, ord);
3461 27668 : op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3462 :
3463 27668 : q = NULL;
3464 27668 : if (T)
3465 : { /* we want < g > = Fp^* */
3466 27668 : if (!equalii(ord,ordp)) {
3467 11906 : q = diviiexact(ord,ordp);
3468 11906 : g = Flxq_pow(g,q,T,p);
3469 : }
3470 : }
3471 27668 : n_q = Fp_log(utoi(a), utoipos(uel(g,2)), op, utoipos(p));
3472 27668 : if (lg(n_q)==1) return gc_leaf(av, n_q);
3473 27668 : if (q) n_q = mulii(q, n_q);
3474 27668 : return gc_INT(av, n_q);
3475 : }
3476 :
3477 : static GEN
3478 519247 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3479 : {
3480 519247 : struct _Flxq *f = (struct _Flxq *)E;
3481 519247 : GEN T = f->T;
3482 519247 : ulong p = f->p;
3483 519247 : long d = get_Flx_degree(T);
3484 519247 : if (Flx_equal1(a)) return gen_0;
3485 359591 : if (Flx_equal(a,g)) return gen_1;
3486 174476 : if (!degpol(a))
3487 35414 : return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3488 139062 : if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3489 139034 : return NULL;
3490 28 : return Flxq_log_index(a, g, ord, T, p);
3491 : }
3492 :
3493 : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3494 :
3495 : const struct bb_group *
3496 283184 : get_Flxq_star(void **E, GEN T, ulong p)
3497 : {
3498 283184 : struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3499 283184 : e->T = T; e->p = p; e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3500 283184 : e->aut = Flx_Frobenius_pre(T, p, e->pi);
3501 283184 : *E = (void*)e; return &Flxq_star;
3502 : }
3503 :
3504 : GEN
3505 97052 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3506 : {
3507 : void *E;
3508 97052 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3509 97052 : return gen_order(a,ord,E,S);
3510 : }
3511 :
3512 : GEN
3513 164217 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3514 : {
3515 : void *E;
3516 164217 : pari_sp av = avma;
3517 164217 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3518 164217 : GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3519 164217 : if (lg(F) > 1 && Flxq_log_use_index(veclast(F), T, p))
3520 24283 : v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3521 164217 : return gc_leaf(av, gen_PH_log(a, g, v, E, S));
3522 : }
3523 :
3524 : static GEN
3525 295304 : Flxq_sumautsum_sqr(void *E, GEN xzd)
3526 : {
3527 295304 : struct _Flxq *D = (struct _Flxq*)E;
3528 295304 : pari_sp av = avma;
3529 : GEN xi, zeta, delta, xi2, zeta2, delta2, temp, xipow;
3530 295304 : GEN T = D->T;
3531 295304 : ulong d, p = D-> p, pi = D->pi;
3532 295304 : xi = gel(xzd, 1); zeta = gel(xzd, 2); delta = gel(xzd, 3);
3533 :
3534 295304 : d = brent_kung_optpow(get_Flx_degree(T)-1,3,1);
3535 295304 : xipow = Flxq_powers_pre(xi, d, T, p, pi);
3536 :
3537 295304 : xi2 = Flx_FlxqV_eval_pre(xi, xipow, T, p, pi);
3538 295304 : zeta2 = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(zeta, xipow, T, p, pi), T, p, pi);
3539 295304 : temp = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(delta, xipow, T, p, pi), T, p, pi);
3540 295304 : delta2 = Flx_add(delta, temp, p);
3541 295304 : return gc_GEN(av, mkvec3(xi2, zeta2, delta2));
3542 : }
3543 :
3544 : static GEN
3545 39750 : Flxq_sumautsum_msqr(void *E, GEN xzd)
3546 : {
3547 39750 : struct _Flxq *D = (struct _Flxq*)E;
3548 39750 : pari_sp av = avma;
3549 : GEN xii, zetai, deltai, xzd2;
3550 39750 : GEN T = D->T, xi0pow = gel(D->aut, 1), zeta0 = gel(D->aut, 2);
3551 39750 : ulong p = D-> p, pi = D->pi;
3552 39750 : xzd2 = Flxq_sumautsum_sqr(E, xzd);
3553 39750 : xii = Flx_FlxqV_eval_pre(gel(xzd2, 1), xi0pow, T, p, pi);
3554 39750 : zetai = Flxq_mul_pre(zeta0, Flx_FlxqV_eval_pre(gel(xzd2, 2), xi0pow, T, p, pi), T, p, pi);
3555 39750 : deltai = Flx_add(gel(xzd2, 3), zetai, p);
3556 :
3557 39750 : return gc_GEN(av, mkvec3(xii, zetai, deltai));
3558 : }
3559 :
3560 : /*returns a + a^(1+s) + a^(1+s+2s) + ... + a^(1+s+...+is)
3561 : where ax = [a,s] with s an automorphism */
3562 : static GEN
3563 209819 : Flxq_sumautsum_pre(GEN ax, long i, GEN T, ulong p, ulong pi) {
3564 209819 : pari_sp av = avma;
3565 : GEN a, xi, zeta, vec, res;
3566 : struct _Flxq D;
3567 : ulong d;
3568 209819 : D.T = Flx_get_red(T, p); D.p = p; D.pi = pi;
3569 209819 : a = gel(ax, 1); xi = gel(ax,2);
3570 209819 : d = brent_kung_optpow(get_Flx_degree(T)-1,2*(hammingu(i)-1),1);
3571 209819 : zeta = Flx_Flxq_eval_pre(a, xi, T, p, pi);
3572 209819 : D.aut = mkvec2(Flxq_powers_pre(xi, d, T, p, pi), zeta);
3573 :
3574 209819 : vec = gen_powu_fold(mkvec3(xi, zeta, zeta), i, (void *)&D, Flxq_sumautsum_sqr, Flxq_sumautsum_msqr);
3575 209819 : res = Flxq_mul_pre(a, Flx_add(pol1_Flx(get_Flx_var(T)), gel(vec, 3), p), T, p, pi);
3576 :
3577 209819 : return gc_GEN(av, res);
3578 : }
3579 :
3580 : /*algorithm from
3581 : Doliskani, J., & Schost, E. (2014).
3582 : Taking roots over high extensions of finite fields
3583 : https://arxiv.org/abs/1110.4350
3584 : */
3585 : static GEN
3586 37674 : Flxq_sqrtl_spec_pre(GEN z, GEN n, GEN T, ulong p, ulong pi, GEN *zetan)
3587 : {
3588 37674 : pari_sp av = avma;
3589 : GEN psn, c, b, new_z, beta, x, y, w, ax, g, zeta;
3590 37674 : long s, l, v = get_Flx_var(T), d = get_Flx_degree(T);
3591 : ulong zeta2, beta2;
3592 37674 : s = itos(Fp_order(utoi(p), stoi(d), n));
3593 37674 : if(s >= d || d % s != 0)
3594 0 : pari_err(e_MISC, "expected p's order mod n to divide the degree of T");
3595 37674 : l = d/s;
3596 37674 : if (!lgpol(z)) return pol0_Flx(get_Flx_var(T));
3597 37674 : T = Flx_get_red(T, p);
3598 37674 : ax = mkvec2(NULL, Flxq_autpow_pre(Flx_Frobenius_pre(T,p,pi), s, T, p,pi));
3599 37674 : psn = diviiexact(subiu(powuu(p, s), 1), n);
3600 : do {
3601 41682 : do c = random_Flx(d, v, p); while (!lgpol(c));
3602 41154 : new_z = Flxq_mul_pre(z, Flxq_pow_pre(c, n, T, p,pi), T, p,pi);
3603 41154 : gel(ax,1) = Flxq_pow_pre(new_z, psn, T, p,pi);
3604 :
3605 : /*If l == 2, b has to be 1 + a^((p^s-1)/n)*/
3606 41154 : if(l == 2) y = gel(ax, 1);
3607 3258 : else y = Flxq_sumautsum_pre(ax, l-2, T, p, pi);
3608 41154 : b = Flx_Fl_add(y, 1, p);
3609 41154 : } while (!lgpol(b));
3610 :
3611 37674 : x = Flxq_mul_pre(new_z, Flxq_pow_pre(b, n, T, p,pi), T, p,pi);
3612 37674 : if(s == 1) {
3613 36519 : if (degpol(x) > 0) return gc_NULL(av);
3614 36482 : beta2 = Fl_sqrtn(Flx_constant(x), umodiu(n, p), p, &zeta2);
3615 36482 : if (beta2==~0UL) return gc_NULL(av);
3616 36482 : if(zetan) *zetan = monomial_Flx(zeta2, 0, get_Flx_var(T));
3617 36482 : w = Flx_Fl_mul(Flxq_inv_pre(Flxq_mul_pre(b, c, T, p,pi), T, p,pi), beta2, p);
3618 36482 : (void)gc_all(av, zetan? 2: 1, &w, zetan);
3619 36482 : return w;
3620 : }
3621 1155 : g = Flxq_minpoly(x, T, p);
3622 1155 : if (degpol(g) > s) return gc_NULL(av);
3623 1155 : beta = Flxq_sqrtn(polx_Flx(get_Flx_var(T)), n, g, p, &zeta);
3624 1155 : if (!beta) return gc_NULL(av);
3625 :
3626 1155 : if(zetan) *zetan = Flx_Flxq_eval(zeta, x, T, p);
3627 1155 : beta = Flx_Flxq_eval(beta, x, T, p);
3628 1155 : w = Flxq_mul_pre(Flxq_inv_pre(Flxq_mul_pre(b, c, T, p,pi), T, p,pi), beta, T, p,pi);
3629 1155 : (void)gc_all(av, zetan? 2: 1, &w, zetan);
3630 1155 : return w;
3631 : }
3632 :
3633 : static GEN
3634 21915 : Flxq_sqrtn_spec_pre(GEN a, GEN n, GEN T, ulong p, ulong pi, GEN q, GEN *zetan)
3635 : {
3636 21915 : pari_sp ltop = avma;
3637 : GEN z, m, u1, u2;
3638 : int is_1;
3639 21915 : if (is_pm1(n))
3640 : {
3641 1925 : if (zetan) *zetan = pol1_Flx(get_Flx_var(T));
3642 1925 : return signe(n) < 0? Flxq_inv_pre(a, T, p,pi): gcopy(a);
3643 : }
3644 19990 : is_1 = gequal1(a);
3645 19990 : if (is_1 && !zetan) return gcopy(a);
3646 19990 : z = pol1_Flx(get_Flx_var(T));
3647 19990 : m = bezout(n,q,&u1,&u2);
3648 19990 : if (!is_pm1(m))
3649 : {
3650 19990 : GEN F = Z_factor(m);
3651 19990 : long i, j, j2 = 0; /* -Wall */
3652 : GEN y, l;
3653 19990 : pari_sp av1 = avma;
3654 40069 : for (i = nbrows(F); i; i--)
3655 : {
3656 20116 : l = gcoeff(F,i,1);
3657 20116 : j = itos(gcoeff(F,i,2));
3658 20116 : if(zetan) {
3659 111 : a = Flxq_sqrtl_spec_pre(a,l,T,p,pi,&y);
3660 148 : if (!a) return gc_NULL(ltop);
3661 111 : j--;
3662 111 : j2 = j;
3663 : }
3664 20116 : if (!is_1 && j > 0) {
3665 : do
3666 : {
3667 37451 : a = Flxq_sqrtl_spec_pre(a,l,T,p,pi,NULL);
3668 37451 : if (!a) return gc_NULL(ltop);
3669 37414 : } while (--j);
3670 : }
3671 : /*This is below finding a's root,
3672 : so we don't spend time doing this, if a is not n-th root*/
3673 20079 : if(zetan) {
3674 223 : for(; j2>0; j2--) y = Flxq_sqrtl_spec_pre(y, l, T, p,pi,NULL);
3675 111 : z = Flxq_mul_pre(z, y, T, p,pi);
3676 : }
3677 20079 : if (gc_needed(ltop,1))
3678 : { /* n can have lots of prime factors*/
3679 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxq_sqrtn_spec");
3680 0 : (void)gc_all(av1, zetan? 2: 1, &a, &z);
3681 : }
3682 : }
3683 : }
3684 :
3685 19953 : if (!equalii(m, n))
3686 434 : a = Flxq_pow_pre(a,modii(u1,q), T, p,pi);
3687 19953 : if (zetan)
3688 : {
3689 111 : *zetan = z;
3690 111 : (void)gc_all(ltop,2,&a,zetan);
3691 : }
3692 : else /* is_1 is 0: a was modified above -> gc_upto valid */
3693 19842 : a = gc_upto(ltop, a);
3694 19953 : return a;
3695 : }
3696 :
3697 : GEN
3698 23128 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3699 : {
3700 23128 : if (!lgpol(a))
3701 : {
3702 7 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3703 0 : if (zeta)
3704 0 : *zeta=pol1_Flx(get_Flx_var(T));
3705 0 : return pol0_Flx(get_Flx_var(T));
3706 : }
3707 23121 : else if(p == 2) {
3708 1206 : pari_sp av = avma;
3709 : GEN z;
3710 1206 : z = F2xq_sqrtn(Flx_to_F2x(a), n, Flx_to_F2x(get_FpX_mod(T)), zeta);
3711 1206 : if (!z) return NULL;
3712 1206 : z = F2x_to_Flx(z);
3713 1206 : if (!zeta) return gc_leaf(av, z);
3714 0 : *zeta=F2x_to_Flx(*zeta);
3715 0 : return gc_all(av, 2, &z,zeta);
3716 : }
3717 : else
3718 : {
3719 : void *E;
3720 21915 : pari_sp av = avma;
3721 21915 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3722 21915 : GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3723 : GEN m, u1, u2, l, zeta2, F, n2, z;
3724 21915 : long i, s, pi, d = get_Flx_degree(T);
3725 21915 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3726 21915 : m = bezout(n,o,&u1,&u2);
3727 21915 : F = Z_factor(m);
3728 46462 : for (i = nbrows(F); i; i--)
3729 : {
3730 24547 : l = gcoeff(F,i,1);
3731 24547 : s = itos(Fp_order(utoi(p), subiu(l, 1), l));
3732 : /*Flxq_sqrtn_spec only works if d > s and s | d
3733 : for those factors of m we use Flxq_sqrtn_spec
3734 : for the other factor we stay with gen_Shanks_sqrtn*/
3735 24547 : if(d <= s || d % s != 0) {
3736 4431 : gcoeff(F,i,2) = gen_0;
3737 : }
3738 20116 : else gcoeff(F,i,2) = stoi(Z_pval(n,l));
3739 : }
3740 21915 : F = factorback(F);
3741 21915 : z = Flxq_sqrtn_spec_pre(a,F,T, p,pi,o,zeta);
3742 21915 : if(!z) return gc_NULL(av);
3743 21878 : n2 = diviiexact(n, F);
3744 21878 : if(!gequal1(n2)) {
3745 5012 : if(zeta) zeta2 = gcopy(*zeta);
3746 5012 : z = gen_Shanks_sqrtn(z, n2, o, zeta, E, S);
3747 5012 : if (!z) return gc_NULL(av);
3748 5012 : if(zeta) *zeta = Flxq_mul_pre(*zeta, zeta2, T, p,pi);
3749 : }
3750 21878 : return gc_all(av, zeta?2:1, &z, zeta);
3751 : }
3752 : }
3753 :
3754 : GEN
3755 230608 : Flxq_sqrt_pre(GEN z, GEN T, ulong p, ulong pi)
3756 : {
3757 230608 : pari_sp av = avma;
3758 : long d;
3759 230608 : if (p==2)
3760 : {
3761 0 : GEN r = F2xq_sqrt(Flx_to_F2x(z), Flx_to_F2x(get_Flx_mod(T)));
3762 0 : return gc_upto(av, F2x_to_Flx(r));
3763 : }
3764 230608 : d = get_Flx_degree(T);
3765 230608 : if (d==2)
3766 : {
3767 65964 : GEN P = get_Flx_mod(T), s;
3768 65964 : ulong c = uel(P,2), b = uel(P,3), a = uel(P,4);
3769 65964 : ulong y = degpol(z)<1 ? 0: uel(z,3);
3770 65964 : if (a==1 && b==0)
3771 15289 : {
3772 16090 : ulong x = degpol(z)<1 ? Flx_constant(z): uel(z,2);
3773 16090 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), Fl_neg(c, p), p, pi);
3774 16090 : if (!r) return gc_NULL(av);
3775 15289 : s = mkvecsmall3(P[1], uel(r,1), uel(r,2));
3776 : }
3777 : else
3778 : {
3779 49874 : ulong b2 = Fl_halve(b, p), t = Fl_div(b2, a, p);
3780 49874 : ulong D = Fl_sub(Fl_sqr(b2, p), Fl_mul(a, c, p), p);
3781 49874 : ulong x = degpol(z)<1 ? Flx_constant(z): Fl_sub(uel(z,2), Fl_mul(uel(z,3), t, p), p);
3782 49874 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), D, p, pi);
3783 49874 : if (!r) return gc_NULL(av);
3784 47480 : s = mkvecsmall3(P[1], Fl_add(uel(r,1), Fl_mul(uel(r,2),t,p), p), uel(r,2));
3785 : }
3786 62769 : return gc_leaf(av, Flx_renormalize(s, 4));
3787 : }
3788 164644 : if (lgpol(z)<=1 && odd(d))
3789 : {
3790 11832 : pari_sp av = avma;
3791 11832 : ulong s = Fl_sqrt(Flx_constant(z), p);
3792 11832 : if (s==~0UL) return gc_NULL(av);
3793 11818 : return gc_GEN(av, Fl_to_Flx(s, get_Flx_var(T)));
3794 : } else
3795 : {
3796 : GEN c, b, new_z, x, y, w, ax;
3797 : ulong p2, beta;
3798 152812 : long v = get_Flx_var(T);
3799 152812 : if (!lgpol(z)) return pol0_Flx(v);
3800 152279 : T = Flx_get_red_pre(T, p, pi);
3801 152279 : ax = mkvec2(NULL, Flx_Frobenius_pre(T, p, pi));
3802 152279 : p2 = p >> 1; /* (p-1) / 2 */
3803 : do {
3804 207178 : do c = random_Flx(d, v, p); while (!lgpol(c));
3805 :
3806 206561 : new_z = Flxq_mul_pre(z, Flxq_sqr_pre(c, T, p, pi), T, p, pi);
3807 206561 : gel(ax, 1) = Flxq_powu_pre(new_z, p2, T, p, pi);
3808 206561 : y = Flxq_sumautsum_pre(ax, d-2, T, p, pi); /* d > 2 */
3809 206561 : b = Flx_Fl_add(y, 1UL, p);
3810 206561 : } while (!lgpol(b));
3811 :
3812 152279 : x = Flxq_mul_pre(new_z, Flxq_sqr_pre(b, T, p, pi), T, p, pi);
3813 152279 : if (degpol(x) > 0) return gc_NULL(av);
3814 145237 : beta = Fl_sqrt_pre(Flx_constant(x), p, pi);
3815 145237 : if (beta==~0UL) return gc_NULL(av);
3816 145237 : w = Flx_Fl_mul(Flxq_inv_pre(Flxq_mul_pre(b, c, T,p,pi), T,p,pi), beta, p);
3817 145237 : return gc_GEN(av, w);
3818 : }
3819 : }
3820 :
3821 : GEN
3822 230608 : Flxq_sqrt(GEN a, GEN T, ulong p)
3823 230608 : { return Flxq_sqrt_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3824 :
3825 : /* assume T irreducible mod p */
3826 : int
3827 402057 : Flxq_issquare(GEN x, GEN T, ulong p)
3828 : {
3829 402057 : if (lgpol(x) == 0 || p == 2) return 1;
3830 395457 : return krouu(Flxq_norm(x,T,p), p) == 1;
3831 : }
3832 :
3833 : /* assume T irreducible mod p */
3834 : int
3835 0 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3836 : {
3837 : pari_sp av;
3838 : GEN m;
3839 0 : if (n==1) return Flxq_issquare(x, T, p);
3840 0 : if (lgpol(x) == 0 || p == 2) return 1;
3841 0 : av = avma;
3842 0 : m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3843 0 : return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3844 : }
3845 :
3846 : GEN
3847 114702 : Flxq_lroot_fast_pre(GEN a, GEN sqx, GEN T, ulong p, ulong pi)
3848 : {
3849 114702 : pari_sp av=avma;
3850 114702 : GEN A = Flx_splitting(a,p);
3851 114702 : return gc_leaf(av, FlxqV_dotproduct_pre(A,sqx,T,p,pi));
3852 : }
3853 : GEN
3854 0 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, ulong p)
3855 0 : { return Flxq_lroot_fast_pre(a, sqx, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3856 :
3857 : GEN
3858 54306 : Flxq_lroot_pre(GEN a, GEN T, ulong p, ulong pi)
3859 : {
3860 54306 : pari_sp av=avma;
3861 54306 : long n = get_Flx_degree(T), d = degpol(a);
3862 : GEN sqx, V;
3863 54306 : if (n==1 || d==-1) return leafcopy(a);
3864 54180 : if (n==2) return Flxq_powu_pre(a, p, T, p, pi);
3865 54180 : sqx = Flxq_autpow_pre(Flx_Frobenius_pre(T, p, pi), n-1, T, p, pi);
3866 54180 : if (d==1 && a[2]==0 && a[3]==1) return gc_leaf(av, sqx);
3867 29071 : if ((ulong) d>=p)
3868 : {
3869 0 : V = Flxq_powers_pre(sqx,p-1,T,p,pi);
3870 0 : return gc_leaf(av, Flxq_lroot_fast_pre(a,V,T,p,pi));
3871 : } else
3872 29071 : return gc_leaf(av, Flx_Flxq_eval_pre(a,sqx,T,p,pi));
3873 : }
3874 : GEN
3875 0 : Flxq_lroot(GEN a, GEN T, ulong p)
3876 0 : { return Flxq_lroot_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3877 :
3878 : ulong
3879 440660 : Flxq_norm(GEN x, GEN TB, ulong p)
3880 : {
3881 440660 : GEN T = get_Flx_mod(TB);
3882 440660 : ulong y = Flx_resultant(T, x, p), L = Flx_lead(T);
3883 440660 : if (L==1 || lgpol(x)==0) return y;
3884 0 : return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3885 : }
3886 :
3887 : ulong
3888 4434 : Flxq_trace(GEN x, GEN TB, ulong p)
3889 : {
3890 4434 : pari_sp av = avma;
3891 : ulong t;
3892 4434 : GEN T = get_Flx_mod(TB);
3893 4434 : long n = degpol(T)-1;
3894 4434 : GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3895 4434 : t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3896 4434 : return gc_ulong(av, t);
3897 : }
3898 :
3899 : /*x must be reduced*/
3900 : GEN
3901 3625 : Flxq_charpoly(GEN x, GEN TB, ulong p)
3902 : {
3903 3625 : pari_sp ltop=avma;
3904 3625 : GEN T = get_Flx_mod(TB);
3905 3625 : long vs = evalvarn(fetch_var());
3906 3625 : GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3907 3625 : GEN r = Flx_FlxY_resultant(T, xm1, p);
3908 3625 : r[1] = x[1];
3909 3625 : (void)delete_var(); return gc_upto(ltop, r);
3910 : }
3911 :
3912 : /* Computing minimal polynomial : */
3913 : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3914 : /* in Algebraic Extensions of Finite Fields' */
3915 :
3916 : /* Let v a linear form, return the linear form z->v(tau*z)
3917 : that is, v*(M_tau) */
3918 :
3919 : static GEN
3920 1758474 : Flxq_transmul_init(GEN tau, GEN T, ulong p, ulong pi)
3921 : {
3922 : GEN bht;
3923 1758474 : GEN h, Tp = get_Flx_red(T, &h);
3924 1758474 : long n = degpol(Tp), vT = Tp[1];
3925 1758474 : GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3926 1758474 : GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3927 1758474 : ft[1] = vT; bt[1] = vT;
3928 1758474 : if (h)
3929 2702 : bht = Flxn_mul_pre(bt, h, n-1, p, pi);
3930 : else
3931 : {
3932 1755772 : GEN bh = Flx_div_pre(Flx_shift(tau, n-1), T, p, pi);
3933 1755772 : bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3934 1755772 : bht[1] = vT;
3935 : }
3936 1758474 : return mkvec3(bt, bht, ft);
3937 : }
3938 :
3939 : static GEN
3940 4244902 : Flxq_transmul(GEN tau, GEN a, long n, ulong p, ulong pi)
3941 : {
3942 4244902 : pari_sp ltop = avma;
3943 : GEN t1, t2, t3, vec;
3944 4244902 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3945 4244902 : if (lgpol(a)==0) return pol0_Flx(a[1]);
3946 4213754 : t2 = Flx_shift(Flx_mul_pre(bt, a, p, pi),1-n);
3947 4213754 : if (lgpol(bht)==0) return gc_leaf(ltop, t2);
3948 3180586 : t1 = Flx_shift(Flx_mul_pre(ft, a, p, pi),-n);
3949 3180586 : t3 = Flxn_mul_pre(t1, bht, n-1, p, pi);
3950 3180586 : vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3951 3180586 : return gc_leaf(ltop, vec);
3952 : }
3953 :
3954 : GEN
3955 815535 : Flxq_minpoly_pre(GEN x, GEN T, ulong p, ulong pi)
3956 : {
3957 815535 : pari_sp ltop = avma;
3958 815535 : long vT = get_Flx_var(T), n = get_Flx_degree(T);
3959 : GEN v_x;
3960 815535 : GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3961 815535 : T = Flx_get_red_pre(T, p, pi);
3962 815535 : v_x = Flxq_powers_pre(x, usqrt(2*n), T, p, pi);
3963 1694772 : while (lgpol(tau) != 0)
3964 : {
3965 : long i, j, m, k1;
3966 : GEN M, v, tr, g_prime, c;
3967 879237 : if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3968 879237 : v = random_Flx(n, vT, p);
3969 879237 : tr = Flxq_transmul_init(tau, T, p, pi);
3970 879237 : v = Flxq_transmul(tr, v, n, p, pi);
3971 879237 : m = 2*(n-degpol(g));
3972 879237 : k1 = usqrt(m);
3973 879237 : tr = Flxq_transmul_init(gel(v_x,k1+1), T, p, pi);
3974 879237 : c = cgetg(m+2,t_VECSMALL);
3975 879237 : c[1] = vT;
3976 4244902 : for (i=0; i<m; i+=k1)
3977 : {
3978 3365665 : long mj = minss(m-i, k1);
3979 13109541 : for (j=0; j<mj; j++)
3980 9743876 : uel(c,m+1-(i+j)) = Flx_dotproduct_pre(v, gel(v_x,j+1), p, pi);
3981 3365665 : v = Flxq_transmul(tr, v, n, p, pi);
3982 : }
3983 879237 : c = Flx_renormalize(c, m+2);
3984 : /* now c contains <v,x^i>, i = 0..m-1 */
3985 879237 : M = Flx_halfgcd_pre(monomial_Flx(1, m, vT), c, p, pi);
3986 879237 : g_prime = gmael(M, 2, 2);
3987 879237 : if (degpol(g_prime) < 1) continue;
3988 866729 : g = Flx_mul_pre(g, g_prime, p, pi);
3989 866729 : tau = Flxq_mul_pre(tau, Flx_FlxqV_eval_pre(g_prime, v_x, T,p,pi), T,p,pi);
3990 : }
3991 815535 : g = Flx_normalize(g,p);
3992 815535 : return gc_leaf(ltop,g);
3993 : }
3994 : GEN
3995 45978 : Flxq_minpoly(GEN x, GEN T, ulong p)
3996 45978 : { return Flxq_minpoly_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3997 :
3998 : GEN
3999 20 : Flxq_conjvec(GEN x, GEN T, ulong p)
4000 : {
4001 20 : long i, l = 1+get_Flx_degree(T);
4002 20 : GEN z = cgetg(l,t_COL);
4003 20 : struct _Flxq D; set_Flxq(&D, T, p);
4004 20 : gel(z,1) = Flx_copy(x);
4005 88 : for (i=2; i<l; i++) gel(z,i) = _Flxq_powu(&D, gel(z,i-1), p);
4006 20 : return z;
4007 : }
4008 :
4009 : GEN
4010 7201 : gener_Flxq(GEN T, ulong p, GEN *po)
4011 : {
4012 7201 : long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
4013 : ulong p_1, pi;
4014 : GEN g, L, L2, o, q, F;
4015 : pari_sp av0, av;
4016 :
4017 7201 : if (f == 1) {
4018 : GEN fa;
4019 28 : o = utoipos(p-1);
4020 28 : fa = Z_factor(o);
4021 28 : L = gel(fa,1);
4022 28 : L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
4023 28 : g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
4024 28 : if (po) *po = mkvec2(o, fa);
4025 28 : return g;
4026 : }
4027 :
4028 7173 : av0 = avma; p_1 = p - 1;
4029 7173 : q = diviuexact(subiu(powuu(p,f), 1), p_1);
4030 :
4031 7173 : L = cgetg(1, t_VECSMALL);
4032 7173 : if (p > 3)
4033 : {
4034 2371 : ulong t = p_1 >> vals(p_1);
4035 2371 : GEN P = gel(factoru(t), 1);
4036 2371 : L = cgetg_copy(P, &i);
4037 3787 : while (--i) L[i] = p_1 / P[i];
4038 : }
4039 7173 : o = factor_pn_1(utoipos(p),f);
4040 7173 : L2 = leafcopy( gel(o, 1) );
4041 19212 : for (i = j = 1; i < lg(L2); i++)
4042 : {
4043 12039 : if (umodui(p_1, gel(L2,i)) == 0) continue;
4044 6488 : gel(L2,j++) = diviiexact(q, gel(L2,i));
4045 : }
4046 7173 : setlg(L2, j); pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4047 7173 : F = Flx_Frobenius_pre(T, p, pi);
4048 17716 : for (av = avma;; set_avma(av))
4049 10543 : {
4050 : GEN tt;
4051 17716 : g = random_Flx(f, vT, p);
4052 17716 : if (degpol(g) < 1) continue;
4053 12102 : if (p == 2) tt = g;
4054 : else
4055 : {
4056 8903 : ulong t = Flxq_norm(g, T, p);
4057 8903 : if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
4058 4781 : tt = Flxq_powu_pre(g, p_1>>1, T, p, pi);
4059 : }
4060 14591 : for (i = 1; i < j; i++)
4061 : {
4062 7418 : GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p, pi);
4063 7418 : if (!degpol(a) && uel(a,2) == p_1) break;
4064 : }
4065 7980 : if (i == j) break;
4066 : }
4067 7173 : if (!po)
4068 : {
4069 187 : set_avma((pari_sp)g);
4070 187 : g = gc_leaf(av0, g);
4071 : }
4072 : else {
4073 6986 : *po = mkvec2(subiu(powuu(p,f), 1), o);
4074 6986 : (void)gc_all(av0, 2, &g, po);
4075 : }
4076 7173 : return g;
4077 : }
4078 :
4079 : static GEN
4080 366572 : _Flxq_neg(void *E, GEN x)
4081 366572 : { struct _Flxq *s = (struct _Flxq *)E;
4082 366572 : return Flx_neg(x,s->p); }
4083 :
4084 : static GEN
4085 1460401 : _Flxq_rmul(void *E, GEN x, GEN y)
4086 1460401 : { struct _Flxq *s = (struct _Flxq *)E;
4087 1460401 : return Flx_mul_pre(x,y,s->p,s->pi); }
4088 :
4089 : static GEN
4090 9460 : _Flxq_inv(void *E, GEN x)
4091 9460 : { struct _Flxq *s = (struct _Flxq *)E;
4092 9460 : return Flxq_inv(x,s->T,s->p); }
4093 :
4094 : static int
4095 69139 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
4096 :
4097 : static GEN
4098 6567 : _Flxq_s(void *E, long x)
4099 6567 : { struct _Flxq *s = (struct _Flxq *)E;
4100 6567 : ulong u = x<0 ? s->p+x: (ulong)x;
4101 6567 : return Fl_to_Flx(u, get_Flx_var(s->T));
4102 : }
4103 :
4104 : static const struct bb_field Flxq_field={_Flxq_red,_Flxq_add,_Flxq_rmul,_Flxq_neg,
4105 : _Flxq_inv,_Flxq_equal0,_Flxq_s};
4106 :
4107 68791 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
4108 : {
4109 68791 : GEN z = new_chunk(sizeof(struct _Flxq));
4110 68791 : set_Flxq((struct _Flxq *)z, T, p); *E = (void*)z; return &Flxq_field;
4111 : }
4112 :
4113 : /***********************************************************************/
4114 : /** Flxn **/
4115 : /***********************************************************************/
4116 :
4117 : GEN
4118 54830 : Flx_invLaplace(GEN x, ulong p)
4119 : {
4120 54830 : long i, d = degpol(x);
4121 : ulong t;
4122 : GEN y;
4123 54830 : if (d <= 1) return Flx_copy(x);
4124 54830 : t = Fl_inv(factorial_Fl(d, p), p);
4125 54830 : y = cgetg(d+3, t_VECSMALL);
4126 54830 : y[1] = x[1];
4127 1343364 : for (i=d; i>=2; i--)
4128 : {
4129 1288534 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
4130 1288534 : t = Fl_mul(t, i, p);
4131 : }
4132 54830 : uel(y,3) = uel(x,3);
4133 54830 : uel(y,2) = uel(x,2);
4134 54830 : return y;
4135 : }
4136 :
4137 : GEN
4138 27566 : Flx_Laplace(GEN x, ulong p)
4139 : {
4140 27566 : long i, d = degpol(x);
4141 27566 : ulong t = 1;
4142 : GEN y;
4143 27566 : if (d <= 1) return Flx_copy(x);
4144 27566 : y = cgetg(d+3, t_VECSMALL);
4145 27566 : y[1] = x[1];
4146 27566 : uel(y,2) = uel(x,2);
4147 27566 : uel(y,3) = uel(x,3);
4148 766234 : for (i=2; i<=d; i++)
4149 : {
4150 738668 : t = Fl_mul(t, i%p, p);
4151 738668 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
4152 : }
4153 27566 : return y;
4154 : }
4155 :
4156 : GEN
4157 6416437 : Flxn_red(GEN a, long n)
4158 : {
4159 6416437 : long i, L, l = lg(a);
4160 : GEN b;
4161 6416437 : if (l == 2 || !n) return zero_Flx(a[1]);
4162 6018958 : L = n+2; if (L > l) L = l;
4163 6018958 : b = cgetg(L, t_VECSMALL); b[1] = a[1];
4164 62479605 : for (i=2; i<L; i++) b[i] = a[i];
4165 6018958 : return Flx_renormalize(b,L);
4166 : }
4167 :
4168 : GEN
4169 5222780 : Flxn_mul_pre(GEN a, GEN b, long n, ulong p, ulong pi)
4170 5222780 : { return Flxn_red(Flx_mul_pre(a, b, p, pi), n); }
4171 : GEN
4172 76054 : Flxn_mul(GEN a, GEN b, long n, ulong p)
4173 76054 : { return Flxn_mul_pre(a, b, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4174 :
4175 : GEN
4176 0 : Flxn_sqr_pre(GEN a, long n, ulong p, ulong pi)
4177 0 : { return Flxn_red(Flx_sqr_pre(a, p, pi), n); }
4178 : GEN
4179 0 : Flxn_sqr(GEN a, long n, ulong p)
4180 0 : { return Flxn_sqr_pre(a, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4181 :
4182 : /* (f*g) \/ x^n */
4183 : static GEN
4184 953825 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p, ulong pi)
4185 953825 : { return Flx_shift(Flx_mul_pre(f, g, p, pi),-n); }
4186 :
4187 : static GEN
4188 529859 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p, ulong pi)
4189 : {
4190 529859 : GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
4191 529859 : return Flx_add(Flx_mulhigh_i(fl, g, n2, p, pi),
4192 : Flxn_mul_pre(fh, g, n - n2, p, pi), p);
4193 : }
4194 :
4195 : /* g==NULL -> assume g==1 */
4196 : GEN
4197 56825 : Flxn_div_pre(GEN g, GEN f, long e, ulong p, ulong pi)
4198 : {
4199 56825 : pari_sp av = avma, av2;
4200 : ulong mask;
4201 : GEN W;
4202 56825 : long n = 1;
4203 56825 : if (lg(f) <= 2) pari_err_INV("Flxn_inv",f);
4204 56825 : W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
4205 56825 : mask = quadratic_prec_mask(e);
4206 56825 : av2 = avma;
4207 272333 : for (;mask>1;)
4208 : {
4209 : GEN u, fr;
4210 215508 : long n2 = n;
4211 215508 : n<<=1; if (mask & 1) n--;
4212 215508 : mask >>= 1;
4213 215508 : fr = Flxn_red(f, n);
4214 215508 : if (mask>1 || !g)
4215 : {
4216 159720 : u = Flxn_mul_pre(W, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
4217 159720 : W = Flx_sub(W, Flx_shift(u, n2), p);
4218 : } else
4219 : {
4220 55788 : GEN y = Flxn_mul_pre(g, W, n, p, pi), yt = Flxn_red(y, n-n2);
4221 55788 : u = Flxn_mul_pre(yt, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
4222 55788 : W = Flx_sub(y, Flx_shift(u, n2), p);
4223 : }
4224 215508 : if (gc_needed(av2,2))
4225 : {
4226 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_div, e = %ld", n);
4227 0 : W = gc_upto(av2, W);
4228 : }
4229 : }
4230 56825 : return gc_upto(av, W);
4231 : }
4232 : GEN
4233 56825 : Flxn_div(GEN g, GEN f, long e, ulong p)
4234 56825 : { return Flxn_div_pre(g, f, e, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4235 :
4236 : GEN
4237 1037 : Flxn_inv(GEN f, long e, ulong p)
4238 1037 : { return Flxn_div(NULL, f, e, p); }
4239 :
4240 : GEN
4241 109615 : Flxn_expint(GEN h, long e, ulong p)
4242 : {
4243 109615 : pari_sp av = avma, av2;
4244 109615 : long v = h[1], n=1;
4245 109615 : GEN f = pol1_Flx(v), g = pol1_Flx(v);
4246 109615 : ulong mask = quadratic_prec_mask(e), pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4247 109615 : av2 = avma;
4248 423966 : for (;mask>1;)
4249 : {
4250 : GEN u, w;
4251 423966 : long n2 = n;
4252 423966 : n<<=1; if (mask & 1) n--;
4253 423966 : mask >>= 1;
4254 423966 : u = Flxn_mul_pre(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p,pi), n-n2, p,pi);
4255 423966 : u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
4256 423966 : w = Flxn_mul_pre(f, Flx_integXn(u, n2-1, p), n-n2, p, pi);
4257 423966 : f = Flx_add(f, Flx_shift(w, n2), p);
4258 423966 : if (mask<=1) break;
4259 314351 : u = Flxn_mul_pre(g, Flxn_mulhigh(f, g, n2, n, p, pi), n-n2, p, pi);
4260 314351 : g = Flx_sub(g, Flx_shift(u, n2), p);
4261 314351 : if (gc_needed(av2,2))
4262 : {
4263 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
4264 0 : (void)gc_all(av2, 2, &f, &g);
4265 : }
4266 : }
4267 109615 : return gc_upto(av, f);
4268 : }
4269 :
4270 : GEN
4271 0 : Flxn_exp(GEN h, long e, ulong p)
4272 : {
4273 0 : if (degpol(h)<1 || uel(h,2)!=0)
4274 0 : pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
4275 0 : return Flxn_expint(Flx_deriv(h, p), e, p);
4276 : }
4277 :
4278 : INLINE GEN
4279 218671 : Flxn_recip(GEN x, long n)
4280 : {
4281 218671 : GEN z=Flx_recipspec(x+2,lgpol(x),n);
4282 218671 : z[1]=x[1];
4283 218671 : return z;
4284 : }
4285 :
4286 : GEN
4287 54528 : Flx_Newton(GEN P, long n, ulong p)
4288 : {
4289 54528 : pari_sp av = avma;
4290 54528 : long d = degpol(P);
4291 54528 : GEN dP = Flxn_recip(Flx_deriv(P, p), d);
4292 54528 : GEN Q = Flxn_div(dP, Flxn_recip(P, d+1), n, p);
4293 54528 : return gc_leaf(av, Q);
4294 : }
4295 :
4296 : GEN
4297 109615 : Flx_fromNewton(GEN P, ulong p)
4298 : {
4299 109615 : pari_sp av = avma;
4300 109615 : ulong n = Flx_constant(P)+1;
4301 109615 : GEN z = Flx_neg(Flx_shift(P, -1), p);
4302 109615 : GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
4303 109615 : return gc_leaf(av, Q);
4304 : }
4305 :
4306 : static void
4307 12521 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
4308 : {
4309 : long i;
4310 : ulong e;
4311 12521 : GEN P = cgetg(d+1, t_VECSMALL);
4312 12521 : GEN V = cgetg(d+1, t_VECSMALL);
4313 1400396 : for (i=1, e=1; i<=d; i++, e++)
4314 : {
4315 1387875 : if (e==p)
4316 : {
4317 459279 : e = 0;
4318 459279 : V[i] = u_lvalrem(i, p, &uel(P,i));
4319 : } else
4320 : {
4321 928596 : V[i] = 0; uel(P,i) = i;
4322 : }
4323 : }
4324 12521 : *pt_P = P; *pt_V = V;
4325 12521 : }
4326 :
4327 : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
4328 : * val large enough to compensate for the power of p in the factorials */
4329 :
4330 : static GEN
4331 504 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long sv)
4332 : {
4333 504 : pari_sp av = avma;
4334 504 : long i, d = n-1, w;
4335 : GEN y, W, E, t;
4336 504 : init_invlaplace(d, p, &E, &W);
4337 504 : t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
4338 504 : w = zv_sum(W);
4339 504 : if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
4340 504 : y = cgetg(d+3,t_POL);
4341 504 : y[1] = evalsigne(1) | sv;
4342 32697 : for (i=d; i>=1; i--)
4343 : {
4344 32193 : gel(y,i+2) = t;
4345 32193 : t = Fp_mulu(t, uel(E,i), q);
4346 32193 : if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
4347 : }
4348 504 : gel(y,2) = t;
4349 504 : return gc_GEN(av, ZX_renormalize(y, d+3));
4350 : }
4351 :
4352 : GEN
4353 27768 : Flx_composedsum(GEN P, GEN Q, ulong p)
4354 : {
4355 27768 : pari_sp av = avma;
4356 27768 : long n = 1 + degpol(P)*degpol(Q);
4357 27768 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4358 27768 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4359 : GEN R;
4360 27768 : if (p >= (ulong)n)
4361 : {
4362 27264 : GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
4363 27264 : GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
4364 27264 : GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
4365 27264 : R = Flx_fromNewton(L, p);
4366 : } else
4367 : {
4368 504 : long v = factorial_lval(n-1, p);
4369 504 : long w = 1 + ulogint(n-1, p);
4370 504 : GEN pv = powuu(p, v);
4371 504 : GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
4372 504 : GEN iL = ZpX_invLaplace_init(n, q, p, v, P[1]);
4373 504 : GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
4374 504 : GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
4375 504 : GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
4376 504 : GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
4377 504 : R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
4378 : }
4379 27768 : return gc_leaf(av, Flx_Fl_mul(R, lead, p));
4380 : }
4381 :
4382 : static GEN
4383 3896 : _Flx_composedsum(void *E, GEN a, GEN b)
4384 3896 : { return Flx_composedsum(a, b, (ulong)E); }
4385 :
4386 : GEN
4387 29009 : FlxV_composedsum(GEN V, ulong p)
4388 29009 : { return gen_product(V, (void *)p, &_Flx_composedsum); }
4389 :
4390 : GEN
4391 0 : Flx_composedprod(GEN P, GEN Q, ulong p)
4392 : {
4393 0 : pari_sp av = avma;
4394 0 : long n = 1+ degpol(P)*degpol(Q);
4395 0 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4396 0 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4397 : GEN R;
4398 0 : if (p >= (ulong)n)
4399 : {
4400 0 : GEN L = Flx_convol(Flx_Newton(P,n,p), Flx_Newton(Q,n,p), p);
4401 0 : R = Flx_fromNewton(L, p);
4402 : } else
4403 : {
4404 0 : long w = 1 + ulogint(n, p);
4405 0 : GEN qf = powuu(p, w);
4406 0 : GEN Pl = FpX_convol(FpX_Newton(Flx_to_ZX(P), n, qf), FpX_Newton(Flx_to_ZX(Q), n, qf), qf);
4407 0 : R = ZX_to_Flx(FpX_fromNewton(Pl, qf), p);
4408 : }
4409 0 : return gc_leaf(av, Flx_Fl_mul(R, lead, p));
4410 :
4411 : }
4412 :
4413 : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
4414 : static GEN
4415 0 : Fl_Xp1_powu(ulong n, ulong p, long v)
4416 : {
4417 0 : ulong k, d = (n + 1) >> 1;
4418 0 : GEN C, V = identity_zv(d);
4419 :
4420 0 : Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
4421 0 : C = cgetg(n+3, t_VECSMALL);
4422 0 : C[1] = v;
4423 0 : uel(C,2) = 1UL;
4424 0 : uel(C,3) = n%p;
4425 0 : uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
4426 : /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
4427 0 : if (SMALL_ULONG(p))
4428 0 : for (k = 3; k <= d; k++)
4429 0 : uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
4430 : else
4431 : {
4432 0 : ulong pi = get_Fl_red(p);
4433 0 : for (k = 3; k <= d; k++)
4434 0 : uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
4435 : }
4436 0 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4437 0 : return C; /* normalized */
4438 : }
4439 :
4440 : /* p arbitrary */
4441 : GEN
4442 675512 : Flx_translate1_basecase(GEN P, ulong p)
4443 : {
4444 675512 : GEN R = Flx_copy(P);
4445 675512 : long i, k, n = degpol(P);
4446 3504957 : for (i = 1; i <= n; i++)
4447 22817731 : for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
4448 675512 : return R;
4449 : }
4450 :
4451 : static int
4452 688677 : translate_basecase(long n, ulong p)
4453 : {
4454 : #ifdef LONG_IS_64BIT
4455 590910 : if (p <= 19) return n < 40;
4456 563070 : if (p < 1UL<<30) return n < 58;
4457 0 : if (p < 1UL<<59) return n < 100;
4458 0 : if (p < 1UL<<62) return n < 120;
4459 0 : if (p < 1UL<<63) return n < 240;
4460 0 : return n < 250;
4461 : #else
4462 97767 : if (p <= 13) return n < 18;
4463 94146 : if (p <= 17) return n < 22;
4464 93524 : if (p <= 29) return n < 39;
4465 91620 : if (p <= 67) return n < 69;
4466 86207 : if (p < 1UL<< 15) return n < 80;
4467 2047 : if (p < 1UL<< 16) return n < 100;
4468 0 : if (p < 1UL<< 28) return n < 300;
4469 0 : return n < 650;
4470 : #endif
4471 : }
4472 : /* assume p prime */
4473 : GEN
4474 663418 : Flx_translate1(GEN P, ulong p)
4475 : {
4476 663418 : long d, n = degpol(P);
4477 : GEN R, Q, S;
4478 663418 : if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
4479 : /* n > 0 */
4480 1148 : d = n >> 1;
4481 1148 : if ((ulong)n < p)
4482 : {
4483 0 : R = Flx_translate1(Flxn_red(P, d), p);
4484 0 : Q = Flx_translate1(Flx_shift(P, -d), p);
4485 0 : S = Fl_Xp1_powu(d, p, P[1]);
4486 0 : return Flx_add(Flx_mul(Q, S, p), R, p);
4487 : }
4488 : else
4489 : {
4490 : ulong q;
4491 1148 : if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
4492 1148 : R = Flx_translate1(Flxn_red(P, q), p);
4493 1148 : Q = Flx_translate1(Flx_shift(P, -q), p);
4494 1148 : S = Flx_add(Flx_shift(Q, q), Q, p);
4495 1148 : return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
4496 : }
4497 : }
4498 :
4499 : GEN
4500 0 : Flx_Fl_translate(GEN P, ulong c, ulong p)
4501 : {
4502 0 : pari_sp av = avma;
4503 : GEN Q;
4504 0 : if (c==0) return Flx_copy(P);
4505 0 : if (c==1) return Flx_translate1(P, p);
4506 0 : Q = Flx_unscale(Flx_translate1(Flx_unscale(P, c, p), p), Fl_inv(c, p), p);
4507 0 : return gc_leaf(av, Q);
4508 : }
4509 :
4510 : static GEN
4511 12017 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
4512 : {
4513 12017 : ulong k, d = n >> 1, c, v = 0;
4514 12017 : GEN C, V, W, U = upowers(p, e-1);
4515 12017 : init_invlaplace(d, p, &V, &W);
4516 12017 : Flv_inv_inplace(V, q);
4517 12017 : C = cgetg(n+3, t_VECSMALL);
4518 12017 : C[1] = vs;
4519 12017 : uel(C,2) = 1UL;
4520 12017 : uel(C,3) = n%q;
4521 12017 : v = u_lvalrem(n, p, &c);
4522 1355682 : for (k = 2; k <= d; k++)
4523 : {
4524 : ulong w;
4525 1343665 : v += u_lvalrem(n-k+1, p, &w) - W[k];
4526 1343665 : c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
4527 1343665 : uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
4528 : }
4529 1374521 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4530 12017 : return C; /* normalized */
4531 : }
4532 :
4533 : GEN
4534 25259 : zlx_translate1(GEN P, ulong p, long e)
4535 : {
4536 25259 : ulong d, q = upowuu(p,e), n = degpol(P);
4537 : GEN R, Q, S;
4538 25259 : if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
4539 : /* n > 0 */
4540 12017 : d = n >> 1;
4541 12017 : R = zlx_translate1(Flxn_red(P, d), p, e);
4542 12017 : Q = zlx_translate1(Flx_shift(P, -d), p, e);
4543 12017 : S = zl_Xp1_powu(d, p, q, e, P[1]);
4544 12017 : return Flx_add(Flx_mul(Q, S, q), R, q);
4545 : }
4546 :
4547 : /***********************************************************************/
4548 : /** Fl2 **/
4549 : /***********************************************************************/
4550 : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
4551 : * a nonsquare D. */
4552 :
4553 : INLINE GEN
4554 27825283 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
4555 :
4556 : /* allow pi = 0 */
4557 : GEN
4558 17131583 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4559 : {
4560 : ulong xaya, xbyb, Db2, mid, z1, z2;
4561 17131583 : ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
4562 17131583 : if (pi)
4563 : {
4564 17131583 : xaya = Fl_mul_pre(x1,y1,p,pi);
4565 17131583 : if (x2==0 && y2==0) return mkF2(xaya,0);
4566 17054423 : if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
4567 12249573 : if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
4568 12249353 : xbyb = Fl_mul_pre(x2,y2,p,pi);
4569 12249353 : mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
4570 12249353 : Db2 = Fl_mul_pre(D, xbyb, p,pi);
4571 : }
4572 0 : else if (p & HIGHMASK)
4573 : {
4574 0 : xaya = Fl_mul(x1,y1,p);
4575 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4576 0 : if (x2==0) return mkF2(xaya,Fl_mul(x1,y2,p));
4577 0 : if (y2==0) return mkF2(xaya,Fl_mul(x2,y1,p));
4578 0 : xbyb = Fl_mul(x2,y2,p);
4579 0 : mid = Fl_mul(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p);
4580 0 : Db2 = Fl_mul(D, xbyb, p);
4581 : }
4582 : else
4583 : {
4584 0 : xaya = (x1 * y1) % p;
4585 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4586 0 : if (x2==0) return mkF2(xaya, (x1 * y2) % p);
4587 0 : if (y2==0) return mkF2(xaya, (x2 * y1) % p);
4588 0 : xbyb = (x2 * y2) % p;
4589 0 : mid = (Fl_add(x1,x2,p) * Fl_add(y1,y2,p)) % p;
4590 0 : Db2 = (D * xbyb) % p;
4591 : }
4592 12249353 : z1 = Fl_add(xaya,Db2,p);
4593 12249353 : z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
4594 12249353 : return mkF2(z1,z2);
4595 : }
4596 :
4597 : /* allow pi = 0 */
4598 : GEN
4599 5022732 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
4600 : {
4601 5022732 : ulong a = x[1], b = x[2];
4602 : ulong a2, Db2, ab;
4603 5022732 : if (pi)
4604 : {
4605 5022732 : a2 = Fl_sqr_pre(a,p,pi);
4606 5022732 : if (b==0) return mkF2(a2,0);
4607 4788110 : Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
4608 4788110 : ab = Fl_mul_pre(a,b,p,pi);
4609 : }
4610 0 : else if (p & HIGHMASK)
4611 : {
4612 0 : a2 = Fl_sqr(a,p);
4613 0 : if (b==0) return mkF2(a2,0);
4614 0 : Db2= Fl_mul(D, Fl_sqr(b,p), p);
4615 0 : ab = Fl_mul(a,b,p);
4616 : }
4617 : else
4618 : {
4619 0 : a2 = (a * a) % p;
4620 0 : if (b==0) return mkF2(a2,0);
4621 0 : Db2= (D * ((b * b) % p)) % p;
4622 0 : ab = (a * b) % p;
4623 : }
4624 4788110 : return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
4625 : }
4626 :
4627 : /* allow pi = 0 */
4628 : ulong
4629 97119260 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
4630 : {
4631 97119260 : ulong a = x[1], b = x[2], a2;
4632 97119260 : if (pi)
4633 : {
4634 97067294 : a2 = Fl_sqr_pre(a,p,pi);
4635 97067294 : return b? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p): a2;
4636 : }
4637 51966 : else if (p & HIGHMASK)
4638 : {
4639 0 : a2 = Fl_sqr(a,p);
4640 0 : return b? Fl_sub(a2, Fl_mul(D, Fl_sqr(b, p), p), p): a2;
4641 : }
4642 : else
4643 : {
4644 51966 : a2 = (a * a) % p;
4645 51966 : return b? Fl_sub(a2, (D * ((b * b) % p)) % p, p): a2;
4646 : }
4647 : }
4648 :
4649 : /* allow pi = 0 */
4650 : GEN
4651 200599 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
4652 : {
4653 200599 : ulong a = x[1], b = x[2], n, ni;
4654 200599 : if (b == 0) return mkF2(Fl_inv(a,p), 0);
4655 166213 : b = Fl_neg(b, p);
4656 166213 : if (pi)
4657 : {
4658 166213 : n = Fl_sub(Fl_sqr_pre(a, p,pi),
4659 : Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p);
4660 166213 : ni = Fl_inv(n,p);
4661 166213 : return mkF2(Fl_mul_pre(a, ni, p,pi), Fl_mul_pre(b, ni, p,pi));
4662 : }
4663 0 : else if (p & HIGHMASK)
4664 : {
4665 0 : n = Fl_sub(Fl_sqr(a, p), Fl_mul(D, Fl_sqr(b, p), p), p);
4666 0 : ni = Fl_inv(n,p);
4667 0 : return mkF2(Fl_mul(a, ni, p), Fl_mul(b, ni, p));
4668 : }
4669 : else
4670 : {
4671 0 : n = Fl_sub((a * a) % p, (D * ((b * b) % p)) % p, p);
4672 0 : ni = Fl_inv(n,p);
4673 0 : return mkF2((a * ni) % p, (b * ni) % p);
4674 : }
4675 : }
4676 :
4677 : int
4678 457455 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
4679 :
4680 : struct _Fl2 {
4681 : ulong p, pi, D;
4682 : };
4683 :
4684 : static GEN
4685 5022732 : _Fl2_sqr(void *data, GEN x)
4686 : {
4687 5022732 : struct _Fl2 *D = (struct _Fl2*)data;
4688 5022732 : return Fl2_sqr_pre(x, D->D, D->p, D->pi);
4689 : }
4690 : static GEN
4691 1963220 : _Fl2_mul(void *data, GEN x, GEN y)
4692 : {
4693 1963220 : struct _Fl2 *D = (struct _Fl2*)data;
4694 1963220 : return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
4695 : }
4696 :
4697 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL; allow pi = 0 */
4698 : GEN
4699 682986 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
4700 : {
4701 682986 : pari_sp av = avma;
4702 : struct _Fl2 d;
4703 : GEN y;
4704 682986 : long s = signe(n);
4705 682986 : if (!s) return mkF2(1,0);
4706 605737 : if (s < 0)
4707 200599 : x = Fl2_inv_pre(x,D,p,pi);
4708 605737 : if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
4709 446225 : d.p = p; d.pi = pi; d.D=D;
4710 446225 : y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
4711 446225 : return gc_leaf(av, y);
4712 : }
4713 :
4714 : static GEN
4715 682986 : _Fl2_pow(void *data, GEN x, GEN n)
4716 : {
4717 682986 : struct _Fl2 *D = (struct _Fl2*)data;
4718 682986 : return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
4719 : }
4720 :
4721 : static GEN
4722 115645 : _Fl2_rand(void *data)
4723 : {
4724 115645 : struct _Fl2 *D = (struct _Fl2*)data;
4725 115645 : ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
4726 115645 : return mkF2(a,b);
4727 : }
4728 :
4729 : GEN
4730 65964 : Fl2_sqrt_pre(GEN z, ulong D, ulong p, ulong pi)
4731 : {
4732 65964 : ulong a = uel(z,1), b = uel(z,2), as2, u, v, s;
4733 65964 : ulong y = Fl_2gener_pre_i(D, p, pi);
4734 65964 : if (b == 0)
4735 19029 : return krouu(a, p)==1 ? mkF2(Fl_sqrt_pre_i(a, y, p, pi), 0)
4736 19029 : : mkF2(0, Fl_sqrt_pre_i(Fl_div(a, D, p), y, p, pi));
4737 52796 : s = Fl_sqrt_pre_i(Fl2_norm_pre(z, D, p, pi), y, p, pi);
4738 52796 : if (s==~0UL) return NULL;
4739 49601 : as2 = Fl_halve(Fl_add(a, s, p), p);
4740 49601 : if (krouu(as2, p)==-1) as2 = Fl_sub(as2, s, p);
4741 49601 : u = Fl_sqrt_pre_i(as2, y, p, pi);
4742 49601 : v = Fl_div(b, Fl_double(u, p), p);
4743 49601 : return mkF2(u,v);
4744 : }
4745 :
4746 : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
4747 : hash_GEN, zv_equal, Fl2_equal1, NULL};
4748 :
4749 : /* allow pi = 0 */
4750 : GEN
4751 77249 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
4752 : {
4753 : struct _Fl2 E;
4754 : GEN o;
4755 77249 : if (a[1]==0 && a[2]==0)
4756 : {
4757 0 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
4758 0 : if (zeta) *zeta=mkF2(1,0);
4759 0 : return zv_copy(a);
4760 : }
4761 77249 : E.p=p; E.pi = pi; E.D = D;
4762 77249 : o = subiu(powuu(p,2), 1);
4763 77249 : return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
4764 : }
4765 :
4766 : /* allow pi = 0 */
4767 : GEN
4768 5214706 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4769 : {
4770 : GEN p1;
4771 5214706 : long i = lg(x)-1;
4772 5214706 : if (i <= 2)
4773 771491 : return mkF2(i == 2? x[2]: 0, 0);
4774 4443215 : p1 = mkF2(x[i], 0);
4775 19611578 : for (i--; i>=2; i--)
4776 : {
4777 15168363 : p1 = Fl2_mul_pre(p1, y, D, p, pi);
4778 15168363 : uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
4779 : }
4780 4443215 : return p1;
4781 : }
4782 :
4783 : /***********************************************************************/
4784 : /** FlxV **/
4785 : /***********************************************************************/
4786 : /* FlxV are t_VEC with Flx coefficients. */
4787 :
4788 : GEN
4789 34482 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
4790 : {
4791 34482 : pari_sp ltop=avma;
4792 : long i;
4793 34482 : GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
4794 257068 : for(i=2;i<lg(V);i++)
4795 222586 : z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
4796 34482 : return gc_leaf(ltop,z);
4797 : }
4798 :
4799 : GEN
4800 0 : ZXV_to_FlxV(GEN x, ulong p)
4801 0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
4802 :
4803 : GEN
4804 3835129 : ZXT_to_FlxT(GEN x, ulong p)
4805 : {
4806 3835129 : if (typ(x) == t_POL)
4807 3774922 : return ZX_to_Flx(x, p);
4808 : else
4809 197303 : pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
4810 : }
4811 :
4812 : GEN
4813 172474 : FlxV_to_Flm(GEN x, long n)
4814 930524 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
4815 :
4816 : GEN
4817 0 : FlxV_red(GEN x, ulong p)
4818 0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4819 :
4820 : GEN
4821 301037 : FlxT_red(GEN x, ulong p)
4822 : {
4823 301037 : if (typ(x) == t_VECSMALL)
4824 202506 : return Flx_red(x, p);
4825 : else
4826 330340 : pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4827 : }
4828 :
4829 : GEN
4830 114702 : FlxqV_dotproduct_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
4831 : {
4832 114702 : long i, lx = lg(x);
4833 : pari_sp av;
4834 : GEN c;
4835 114702 : if (lx == 1) return pol0_Flx(get_Flx_var(T));
4836 114702 : av = avma; c = Flx_mul_pre(gel(x,1),gel(y,1), p, pi);
4837 489174 : for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4838 114702 : return gc_leaf(av, Flx_rem_pre(c,T,p,pi));
4839 : }
4840 : GEN
4841 0 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4842 0 : { return FlxqV_dotproduct_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4843 :
4844 : GEN
4845 1932 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4846 : {
4847 1932 : long i, l = minss(lg(x), lg(y));
4848 : ulong pi;
4849 : pari_sp av;
4850 : GEN c;
4851 1932 : if (l == 2) return pol0_Flx(get_Flx_var(T));
4852 1919 : av = avma; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4853 1919 : c = Flx_mul_pre(gel(x,2),gel(y,2), p, pi);
4854 6267 : for (i=3; i<l; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4855 1919 : return gc_leaf(av, Flx_rem_pre(c,T,p,pi));
4856 : }
4857 :
4858 : /* allow pi = 0 */
4859 : GEN
4860 325162 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4861 : {
4862 325162 : long i, l = lg(z);
4863 325162 : GEN y = cgetg(l, t_VECSMALL);
4864 15070763 : for (i=1; i<l; i++) uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4865 325162 : return y;
4866 : }
4867 :
4868 : /***********************************************************************/
4869 : /** FlxM **/
4870 : /***********************************************************************/
4871 : /* allow pi = 0 */
4872 : GEN
4873 22302 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4874 : {
4875 22302 : long i, l = lg(z);
4876 22302 : GEN y = cgetg(l, t_MAT);
4877 347464 : for (i=1; i<l; i++) gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4878 22302 : return y;
4879 : }
4880 :
4881 : GEN
4882 0 : zero_FlxC(long n, long sv)
4883 : {
4884 0 : GEN x = cgetg(n + 1, t_COL), z = zero_Flx(sv);
4885 : long i;
4886 0 : for (i = 1; i <= n; i++) gel(x, i) = z;
4887 0 : return x;
4888 : }
4889 :
4890 : GEN
4891 0 : FlxC_neg(GEN x, ulong p)
4892 0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4893 :
4894 : GEN
4895 0 : FlxC_sub(GEN x, GEN y, ulong p)
4896 0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4897 :
4898 : GEN
4899 0 : zero_FlxM(long r, long c, long sv)
4900 : {
4901 0 : GEN x = cgetg(c + 1, t_MAT), z = zero_FlxC(r, sv);
4902 : long j;
4903 0 : for (j = 1; j <= c; j++) gel(x, j) = z;
4904 0 : return x;
4905 : }
4906 :
4907 : GEN
4908 0 : zero_FlxM_copy(long r, long c, long sv)
4909 : {
4910 0 : GEN x = cgetg(c + 1, t_MAT);
4911 : long j;
4912 0 : for (j = 1; j <= c; j++) gel(x, j) = zero_FlxC(r, sv);
4913 0 : return x;
4914 : }
4915 :
4916 : GEN
4917 0 : FlxM_neg(GEN x, ulong p)
4918 0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4919 :
4920 : GEN
4921 0 : FlxM_sub(GEN x, GEN y, ulong p)
4922 0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4923 :
4924 : GEN
4925 0 : FlxC_Fl_translate(GEN x, ulong c, ulong p)
4926 0 : { pari_APPLY_type(t_COL, Flx_Fl_translate(gel(x,i), c, p)) }
4927 :
4928 : GEN
4929 0 : FlxM_Fl_translate(GEN x, ulong c, ulong p)
4930 0 : { pari_APPLY_same(FlxC_Fl_translate(gel(x,i), c, p)) }
4931 :
4932 : GEN
4933 247515 : FlxqC_red_pre(GEN x, GEN T, ulong p, ulong pi)
4934 4988454 : { pari_APPLY_type(t_COL, Flx_rem_pre(gel(x,i), T, p, pi)) }
4935 :
4936 : GEN
4937 82442 : FlxqM_red_pre(GEN x, GEN T, ulong p, ulong pi)
4938 329957 : { pari_APPLY_same(FlxqC_red_pre(gel(x,i), T, p, pi)) }
4939 :
4940 : GEN
4941 0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4942 0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4943 :
4944 : GEN
4945 0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4946 0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4947 :
4948 : static GEN
4949 47865 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4950 : long i, j, l, lc;
4951 47865 : GEN N = cgetg_copy(M, &l), x;
4952 47865 : if (l == 1)
4953 0 : return N;
4954 47865 : lc = lgcols(M);
4955 213857 : for (j = 1; j < l; j++) {
4956 165992 : gel(N, j) = cgetg(lc, t_COL);
4957 1146243 : for (i = 1; i < lc; i++) {
4958 980251 : x = gcoeff(M, i, j);
4959 980251 : gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4960 : }
4961 : }
4962 47865 : return N;
4963 : }
4964 :
4965 : static GEN
4966 891277 : kron_pack_Flx_spec_half(GEN x, long l) {
4967 891277 : if (l == 0) return gen_0;
4968 544951 : return Flx_to_int_halfspec(x, l);
4969 : }
4970 :
4971 : static GEN
4972 85585 : kron_pack_Flx_spec(GEN x, long l) {
4973 : long i;
4974 : GEN w, y;
4975 85585 : if (l == 0)
4976 28245 : return gen_0;
4977 57340 : y = cgetipos(l + 2);
4978 186172 : for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4979 128832 : *w = x[i];
4980 57340 : return y;
4981 : }
4982 :
4983 : static GEN
4984 3389 : kron_pack_Flx_spec_2(GEN x, long l) { return Flx_eval2BILspec(x, 2, l); }
4985 :
4986 : static GEN
4987 0 : kron_pack_Flx_spec_3(GEN x, long l) { return Flx_eval2BILspec(x, 3, l); }
4988 :
4989 : static GEN
4990 78984 : kron_unpack_Flx(GEN z, ulong p)
4991 : {
4992 78984 : long i, l = lgefint(z);
4993 78984 : GEN x = cgetg(l, t_VECSMALL), w;
4994 255639 : for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4995 176655 : x[i] = ((ulong) *w) % p;
4996 78984 : return Flx_renormalize(x, l);
4997 : }
4998 :
4999 : static GEN
5000 2930 : kron_unpack_Flx_2(GEN x, ulong p) {
5001 2930 : long d = (lgefint(x)-1)/2 - 1;
5002 2930 : return Z_mod2BIL_Flx_2(x, d, p);
5003 : }
5004 :
5005 : static GEN
5006 0 : kron_unpack_Flx_3(GEN x, ulong p) {
5007 0 : long d = lgefint(x)/3 - 1;
5008 0 : return Z_mod2BIL_Flx_3(x, d, p);
5009 : }
5010 :
5011 : static GEN
5012 116931 : FlxM_pack_ZM_bits(GEN M, long b)
5013 : {
5014 : long i, j, l, lc;
5015 116931 : GEN N = cgetg_copy(M, &l), x;
5016 116931 : if (l == 1)
5017 0 : return N;
5018 116931 : lc = lgcols(M);
5019 493598 : for (j = 1; j < l; j++) {
5020 376667 : gel(N, j) = cgetg(lc, t_COL);
5021 6501067 : for (i = 1; i < lc; i++) {
5022 6124400 : x = gcoeff(M, i, j);
5023 6124400 : gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
5024 : }
5025 : }
5026 116931 : return N;
5027 : }
5028 :
5029 : static GEN
5030 23936 : ZM_unpack_FlxM(GEN M, ulong p, ulong sv, GEN (*unpack)(GEN, ulong))
5031 : {
5032 : long i, j, l, lc;
5033 23936 : GEN N = cgetg_copy(M, &l), x;
5034 23936 : if (l == 1)
5035 0 : return N;
5036 23936 : lc = lgcols(M);
5037 116098 : for (j = 1; j < l; j++) {
5038 92162 : gel(N, j) = cgetg(lc, t_COL);
5039 902371 : for (i = 1; i < lc; i++) {
5040 810209 : x = unpack(gcoeff(M, i, j), p);
5041 810209 : x[1] = sv;
5042 810209 : gcoeff(N, i, j) = x;
5043 : }
5044 : }
5045 23936 : return N;
5046 : }
5047 :
5048 : static GEN
5049 58506 : ZM_unpack_FlxM_bits(GEN M, long b, ulong p, ulong pi, long sv)
5050 : {
5051 : long i, j, l, lc;
5052 58506 : GEN N = cgetg_copy(M, &l), x;
5053 58506 : if (l == 1)
5054 0 : return N;
5055 58506 : lc = lgcols(M);
5056 58506 : if (b < BITS_IN_LONG) {
5057 203927 : for (j = 1; j < l; j++) {
5058 147148 : gel(N, j) = cgetg(lc, t_COL);
5059 3910526 : for (i = 1; i < lc; i++) {
5060 3763378 : x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
5061 3763378 : x[1] = sv;
5062 3763378 : gcoeff(N, i, j) = x;
5063 : }
5064 : }
5065 : } else {
5066 1727 : if (!pi) pi = get_Fl_red(p); /* unset if !SMALL_ULONG(p) */
5067 9932 : for (j = 1; j < l; j++) {
5068 8205 : gel(N, j) = cgetg(lc, t_COL);
5069 175557 : for (i = 1; i < lc; i++) {
5070 167352 : x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
5071 167352 : x[1] = sv;
5072 167352 : gcoeff(N, i, j) = x;
5073 : }
5074 : }
5075 : }
5076 58506 : return N;
5077 : }
5078 :
5079 : static GEN
5080 82442 : FlxM_mul_Kronecker_i(GEN A, GEN B, ulong p, ulong pi, long d, long sv)
5081 : {
5082 82442 : long b, n = lg(A) - 1;
5083 : GEN C, z;
5084 : GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
5085 82442 : int is_sqr = A==B;
5086 :
5087 82442 : z = muliu(muliu(sqru(p - 1), d), n);
5088 82442 : b = expi(z) + 1;
5089 : /* only do expensive bit-packing if it saves at least 1 limb */
5090 82442 : if (b <= BITS_IN_HALFULONG)
5091 77887 : { if (nbits2nlong(d*b) == (d + 1)/2) b = BITS_IN_HALFULONG; }
5092 : else
5093 : {
5094 4555 : long l = lgefint(z) - 2;
5095 4555 : if (nbits2nlong(d*b) == d*l) b = l*BITS_IN_LONG;
5096 : }
5097 :
5098 82442 : switch (b) {
5099 22849 : case BITS_IN_HALFULONG:
5100 22849 : pack = kron_pack_Flx_spec_half;
5101 22849 : unpack = int_to_Flx_half;
5102 22849 : break;
5103 1038 : case BITS_IN_LONG:
5104 1038 : pack = kron_pack_Flx_spec;
5105 1038 : unpack = kron_unpack_Flx;
5106 1038 : break;
5107 49 : case 2*BITS_IN_LONG:
5108 49 : pack = kron_pack_Flx_spec_2;
5109 49 : unpack = kron_unpack_Flx_2;
5110 49 : break;
5111 0 : case 3*BITS_IN_LONG:
5112 0 : pack = kron_pack_Flx_spec_3;
5113 0 : unpack = kron_unpack_Flx_3;
5114 0 : break;
5115 58506 : default:
5116 58506 : A = FlxM_pack_ZM_bits(A, b);
5117 58506 : B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
5118 58506 : C = ZM_mul(A, B);
5119 58506 : return ZM_unpack_FlxM_bits(C, b, p, pi, sv);
5120 : }
5121 23936 : A = FlxM_pack_ZM(A, pack);
5122 23936 : B = is_sqr? A: FlxM_pack_ZM(B, pack);
5123 23936 : C = ZM_mul(A, B);
5124 23936 : return ZM_unpack_FlxM(C, p, sv, unpack);
5125 : }
5126 :
5127 : GEN
5128 82442 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
5129 : {
5130 82442 : pari_sp av = avma;
5131 82442 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
5132 82442 : long sv = get_Flx_var(T), d = get_Flx_degree(T);
5133 82442 : GEN C = FlxM_mul_Kronecker_i(A, B, p, pi, d, sv);
5134 82442 : C = FlxqM_red_pre(C, T, p, pi);
5135 82442 : return gc_upto(av, C);
5136 : }
5137 :
5138 : /* assume m > 1 */
5139 : static long
5140 0 : FlxV_max_degree_i(GEN x, long m)
5141 : {
5142 0 : long i, l = degpol(gel(x,1));
5143 0 : for (i = 2; i < m; i++) l = maxss(l, degpol(gel(x,i)));
5144 0 : return l;
5145 : }
5146 :
5147 : /* assume n > 1 and m > 1 */
5148 : static long
5149 0 : FlxM_max_degree_i(GEN x, long n, long m)
5150 : {
5151 0 : long j, l = FlxV_max_degree_i(gel(x,1), m);
5152 0 : for (j = 2; j < n; j++) l = maxss(l, FlxV_max_degree_i(gel(x,j), m));
5153 0 : return l;
5154 : }
5155 :
5156 : static long
5157 0 : FlxM_max_degree(GEN x)
5158 : {
5159 0 : long n = lg(x), m;
5160 0 : if (n == 1) return -1;
5161 0 : m = lgcols(x); return m == 1? -1: FlxM_max_degree_i(x, n, m);
5162 : }
5163 :
5164 : GEN
5165 0 : FlxM_mul(GEN x, GEN y, ulong p)
5166 : {
5167 0 : pari_sp av = avma;
5168 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
5169 : long sv, d;
5170 0 : if (lg(x) == 1) return cgetg(1,t_MAT);
5171 0 : if (lg(gel(x,1))==1) return FlxqM_mul(x, y, NULL, p);
5172 0 : sv = mael3(x,1,1,1);
5173 0 : d = maxss(FlxM_max_degree(x), FlxM_max_degree(y));
5174 0 : return gc_GEN(av, FlxM_mul_Kronecker_i(x, y, p, pi, d+1, sv));
5175 : }
|
