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 977099125 : get_Flx_red(GEN T, GEN *B)
22 : {
23 977099125 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 659608 : *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 37121917 : Flx_to_ZX(GEN z)
43 : {
44 37121917 : long i, l = lg(z);
45 37121917 : GEN x = cgetg(l,t_POL);
46 242448588 : for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
47 37108746 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
48 : }
49 :
50 : GEN
51 71360 : Flx_to_FlxX(GEN z, long sv)
52 : {
53 71360 : long i, l = lg(z);
54 71360 : GEN x = cgetg(l,t_POL);
55 278207 : for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
56 71360 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
57 : }
58 :
59 : /* same as Flx_to_ZX, in place */
60 : GEN
61 36445520 : Flx_to_ZX_inplace(GEN z)
62 : {
63 36445520 : long i, l = lg(z);
64 227402082 : for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
65 36437476 : 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 65823031 : Flx_to_Flv(GEN x, long N)
71 : {
72 65823031 : GEN z = cgetg(N+1,t_VECSMALL);
73 65817098 : long i, l = lg(x)-1;
74 65817098 : x++;
75 704903008 : for (i=1; i<l ; i++) z[i]=x[i];
76 328324868 : for ( ; i<=N; i++) z[i]=0;
77 65817098 : return z;
78 : }
79 :
80 : /*Flv_to_Flx=zv_to_zx*/
81 : GEN
82 25249869 : Flv_to_Flx(GEN x, long sv)
83 : {
84 25249869 : long i, l=lg(x)+1;
85 25249869 : GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
86 25245522 : x--;
87 278330621 : for (i=2; i<l ; i++) z[i]=x[i];
88 25245522 : return Flx_renormalize(z,l);
89 : }
90 :
91 : /*Flm_to_FlxV=zm_to_zxV*/
92 : GEN
93 2296 : Flm_to_FlxV(GEN x, long sv)
94 6272 : { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
95 :
96 : /*FlxC_to_ZXC=zxC_to_ZXC*/
97 : GEN
98 103966 : FlxC_to_ZXC(GEN x)
99 527193 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
100 :
101 : /*FlxC_to_ZXC=zxV_to_ZXV*/
102 : GEN
103 611267 : FlxV_to_ZXV(GEN x)
104 2472914 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
105 :
106 : void
107 2928567 : FlxV_to_ZXV_inplace(GEN v)
108 : {
109 : long i;
110 7778793 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
111 2928470 : }
112 :
113 : /*FlxM_to_ZXM=zxM_to_ZXM*/
114 : GEN
115 2399 : FlxM_to_ZXM(GEN x)
116 8123 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
117 :
118 : GEN
119 398050 : FlxV_to_FlxX(GEN x, long v)
120 : {
121 398050 : long i, l = lg(x)+1;
122 398050 : GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
123 398050 : x--;
124 4993804 : for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
125 398050 : 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 19844912 : Fl_to_Flx(ulong x, long sv) { return x? mkvecsmall2(sv, x): pol0_Flx(sv); }
156 :
157 : /* a X^d */
158 : GEN
159 913817 : monomial_Flx(ulong a, long d, long vs)
160 : {
161 : GEN P;
162 913817 : if (a==0) return pol0_Flx(vs);
163 913817 : P = const_vecsmall(d+2, 0);
164 913820 : P[1] = vs; P[d+2] = a; return P;
165 : }
166 :
167 : GEN
168 2595739 : Z_to_Flx(GEN x, ulong p, long sv)
169 : {
170 2595739 : long u = umodiu(x,p);
171 2595732 : 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 167455982 : ZX_to_Flx(GEN x, ulong p)
177 : {
178 167455982 : long i, lx = lg(x);
179 167455982 : GEN a = cgetg(lx, t_VECSMALL);
180 167414685 : a[1]=((ulong)x[1])&VARNBITS;
181 1111285461 : for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
182 167424346 : 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 6083880 : zx_to_Flx(GEN x, ulong p)
188 : {
189 6083880 : long i, lx = lg(x);
190 6083880 : GEN a = cgetg(lx, t_VECSMALL);
191 6077615 : a[1] = x[1];
192 18664959 : for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
193 6077512 : return Flx_renormalize(a,lx);
194 : }
195 :
196 : ulong
197 73024496 : Rg_to_Fl(GEN x, ulong p)
198 : {
199 73024496 : switch(typ(x))
200 : {
201 48030113 : case t_INT: return umodiu(x, p);
202 454302 : case t_FRAC: {
203 454302 : ulong z = umodiu(gel(x,1), p);
204 454302 : if (!z) return 0;
205 444608 : return Fl_div(z, umodiu(gel(x,2), p), p);
206 : }
207 205948 : case t_PADIC: return padic_to_Fl(x, p);
208 24334139 : case t_INTMOD: {
209 24334139 : GEN q = gel(x,1), a = gel(x,2);
210 24334139 : 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 1706940 : Rg_to_F2(GEN x)
221 : {
222 1706940 : switch(typ(x))
223 : {
224 274131 : 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(gel(x,2),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 1432809 : case t_INTMOD: {
233 1432809 : GEN q = gel(x,1), a = gel(x,2);
234 1432809 : if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
235 1432809 : 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 2354930 : RgX_to_Flx(GEN x, ulong p)
244 : {
245 2354930 : long i, lx = lg(x);
246 2354930 : GEN a = cgetg(lx, t_VECSMALL);
247 2354930 : a[1]=((ulong)x[1])&VARNBITS;
248 20434878 : for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
249 2354930 : return Flx_renormalize(a,lx);
250 : }
251 :
252 : GEN
253 7 : RgXV_to_FlxV(GEN x, ulong p)
254 175 : { 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 3566028 : Rg_to_Flxq(GEN x, GEN T, ulong p)
259 : {
260 3566028 : long ta, tx = typ(x), v = get_Flx_var(T);
261 : ulong pi;
262 : GEN a, b;
263 3566028 : if (is_const_t(tx))
264 : {
265 3315401 : if (tx == t_FFELT) return FF_to_Flxq(x);
266 2584393 : return Fl_to_Flx(Rg_to_Fl(x, p), v);
267 : }
268 250627 : 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 242051 : case t_POL:
280 242051 : x = RgX_to_Flx(x,p);
281 242051 : if (x[1] != v) break;
282 242051 : 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 2117928519 : Flx_renormalize(GEN /*in place*/ x, long lx)
298 : {
299 : long i;
300 2366584306 : for (i = lx-1; i>1; i--)
301 2272063167 : if (x[i]) break;
302 2117928519 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
303 2116639605 : setlg(x, i+1); return x;
304 : }
305 :
306 : GEN
307 1884139 : Flx_red(GEN z, ulong p)
308 : {
309 1884139 : long i, l = lg(z);
310 1884139 : GEN x = cgetg(l, t_VECSMALL);
311 1883972 : x[1] = z[1];
312 33979167 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
313 1883972 : return Flx_renormalize(x,l);
314 : }
315 :
316 : int
317 28890400 : Flx_equal(GEN V, GEN W)
318 : {
319 28890400 : long l = lg(V);
320 28890400 : if (lg(W) != l) return 0;
321 29917995 : while (--l > 1) /* do not compare variables, V[1] */
322 28784246 : if (V[l] != W[l]) return 0;
323 1133749 : return 1;
324 : }
325 :
326 : GEN
327 2589330 : random_Flx(long d1, long vs, ulong p)
328 : {
329 2589330 : long i, d = d1+2;
330 2589330 : GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
331 17931187 : for (i=2; i<d; i++) y[i] = random_Fl(p);
332 2589520 : return Flx_renormalize(y,d);
333 : }
334 :
335 : static GEN
336 7125885 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
337 : {
338 : long i,lz;
339 : GEN z;
340 :
341 7125885 : if (ly>lx) swapspec(x,y, lx,ly);
342 7125885 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
343 106030468 : for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
344 89789861 : for ( ; i<lx; i++) z[i+2] = x[i];
345 7125885 : z[1] = 0; return Flx_renormalize(z, lz);
346 : }
347 :
348 : GEN
349 62559028 : Flx_add(GEN x, GEN y, ulong p)
350 : {
351 : long i,lz;
352 : GEN z;
353 62559028 : long lx=lg(x);
354 62559028 : long ly=lg(y);
355 62559028 : if (ly>lx) swapspec(x,y, lx,ly);
356 62559028 : lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
357 572659061 : for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
358 128037228 : for ( ; i<lx; i++) z[i] = x[i];
359 62538569 : return Flx_renormalize(z, lz);
360 : }
361 :
362 : GEN
363 9902218 : Flx_Fl_add(GEN y, ulong x, ulong p)
364 : {
365 : GEN z;
366 : long lz, i;
367 9902218 : if (!lgpol(y))
368 228672 : return Fl_to_Flx(x,y[1]);
369 9674578 : lz=lg(y);
370 9674578 : z=cgetg(lz,t_VECSMALL);
371 9674213 : z[1]=y[1];
372 9674213 : z[2] = Fl_add(y[2],x,p);
373 46904993 : for(i=3;i<lz;i++)
374 37231084 : z[i] = y[i];
375 9673909 : if (lz==3) z = Flx_renormalize(z,lz);
376 9673678 : return z;
377 : }
378 :
379 : static GEN
380 898662 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
381 : {
382 : long i,lz;
383 : GEN z;
384 :
385 898662 : if (ly <= lx)
386 : {
387 898673 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
388 53986659 : for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
389 1449000 : 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 898339 : z[1] = 0; return Flx_renormalize(z, lz);
398 : }
399 :
400 : GEN
401 138071607 : Flx_sub(GEN x, GEN y, ulong p)
402 : {
403 138071607 : long i,lz,lx = lg(x), ly = lg(y);
404 : GEN z;
405 :
406 138071607 : if (ly <= lx)
407 : {
408 87904522 : lz = lx; z = cgetg(lz, t_VECSMALL);
409 456513344 : for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
410 175719307 : for ( ; i<lx; i++) z[i] = x[i];
411 : }
412 : else
413 : {
414 50167085 : lz = ly; z = cgetg(lz, t_VECSMALL);
415 259819454 : for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
416 232516988 : for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
417 : }
418 138064306 : z[1]=x[1]; return Flx_renormalize(z, lz);
419 : }
420 :
421 : GEN
422 151355 : Flx_Fl_sub(GEN y, ulong x, ulong p)
423 : {
424 : GEN z;
425 151355 : long lz = lg(y), i;
426 151355 : if (lz==2)
427 513 : return Fl_to_Flx(Fl_neg(x, p),y[1]);
428 150842 : z = cgetg(lz, t_VECSMALL);
429 150842 : z[1] = y[1];
430 150842 : z[2] = Fl_sub(uel(y,2), x, p);
431 751654 : for(i=3; i<lz; i++)
432 600812 : z[i] = y[i];
433 150842 : if (lz==3) z = Flx_renormalize(z,lz);
434 150842 : return z;
435 : }
436 :
437 : static GEN
438 3265163 : Flx_negspec(GEN x, ulong p, long l)
439 : {
440 : long i;
441 3265163 : GEN z = cgetg(l+2, t_VECSMALL) + 2;
442 21004845 : for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
443 3265086 : return z-2;
444 : }
445 :
446 : GEN
447 3265162 : Flx_neg(GEN x, ulong p)
448 : {
449 3265162 : GEN z = Flx_negspec(x+2, p, lgpol(x));
450 3265197 : z[1] = x[1];
451 3265197 : return z;
452 : }
453 :
454 : GEN
455 1746563 : Flx_neg_inplace(GEN x, ulong p)
456 : {
457 1746563 : long i, l = lg(x);
458 52077893 : for (i=2; i<l; i++)
459 50331330 : if (x[i]) x[i] = p - x[i];
460 1746563 : return x;
461 : }
462 :
463 : GEN
464 2444284 : Flx_double(GEN y, ulong p)
465 : {
466 : long i, l;
467 2444284 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
468 20390198 : for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
469 2444284 : return Flx_renormalize(z, l);
470 : }
471 : GEN
472 1049263 : Flx_triple(GEN y, ulong p)
473 : {
474 : long i, l;
475 1049263 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
476 8306859 : for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
477 1049263 : return Flx_renormalize(z, l);
478 : }
479 :
480 : GEN
481 18376576 : Flx_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
482 : {
483 : GEN z;
484 : long i, l;
485 18376576 : if (!x) return pol0_Flx(y[1]);
486 17594874 : z = cgetg_copy(y, &l); z[1] = y[1];
487 17594740 : if (pi==0)
488 : {
489 15414892 : if (HIGHWORD(x | p))
490 0 : for(i=2; i<l; i++) z[i] = Fl_mul(uel(y,i), x, p);
491 : else
492 92771755 : for(i=2; i<l; i++) z[i] = (uel(y,i) * x) % p;
493 : } else
494 17970208 : for(i=2; i<l; i++) z[i] = Fl_mul_pre(uel(y,i), x, p, pi);
495 17595518 : return Flx_renormalize(z, l);
496 : }
497 :
498 : GEN
499 7279171 : Flx_Fl_mul(GEN x, ulong y, ulong p)
500 7279171 : { 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 11965822 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
515 : {
516 : GEN z;
517 : long i, l;
518 11965822 : z = cgetg_copy(y, &l); z[1] = y[1];
519 11962585 : if (HIGHWORD(x | p))
520 5408429 : for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
521 : else
522 26839460 : for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
523 11962579 : 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 26807122 : Flx_shift(GEN a, long n)
529 : {
530 26807122 : long i, l = lg(a);
531 : GEN b;
532 26807122 : if (l==2 || !n) return Flx_copy(a);
533 26465134 : if (l+n<=2) return pol0_Flx(a[1]);
534 26250504 : b = cgetg(l+n, t_VECSMALL);
535 26248365 : b[1] = a[1];
536 26248365 : if (n < 0)
537 71716291 : for (i=2-n; i<l; i++) b[i+n] = a[i];
538 : else
539 : {
540 50889733 : for (i=0; i<n; i++) b[2+i] = 0;
541 148354841 : for (i=2; i<l; i++) b[i+n] = a[i];
542 : }
543 26248365 : return b;
544 : }
545 :
546 : GEN
547 62175119 : Flx_normalize(GEN z, ulong p)
548 : {
549 62175119 : long l = lg(z)-1;
550 62175119 : ulong p1 = z[l]; /* leading term */
551 62175119 : if (p1 == 1) return z;
552 11943610 : 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 3662141 : Flx_addshift(GEN x, GEN y, ulong p, long d)
558 : {
559 3662141 : GEN xd,yd,zd = (GEN)avma;
560 3662141 : long a,lz,ny = lgpol(y), nx = lgpol(x);
561 3662141 : long vs = x[1];
562 3662141 : if (nx == 0) return y;
563 3660289 : x += 2; y += 2; a = ny-d;
564 3660289 : if (a <= 0)
565 : {
566 85096 : lz = (a>nx)? ny+2: nx+d+2;
567 85096 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
568 1725373 : while (xd > x) *--zd = *--xd;
569 85096 : x = zd + a;
570 163436 : while (zd > x) *--zd = 0;
571 : }
572 : else
573 : {
574 3575193 : xd = new_chunk(d); yd = y+d;
575 3575193 : x = Flx_addspec(x,yd,p, nx,a);
576 3575193 : lz = (a>nx)? ny+2: lg(x)+d;
577 132058374 : x += 2; while (xd > x) *--zd = *--xd;
578 : }
579 60051102 : while (yd > y) *--zd = *--yd;
580 3660289 : *--zd = vs;
581 3660289 : *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
582 : }
583 :
584 : /* shift polynomial + gerepile */
585 : /* Do not set evalvarn*/
586 : static GEN
587 630003999 : Flx_shiftip(pari_sp av, GEN x, long v)
588 : {
589 630003999 : long i, lx = lg(x), ly;
590 : GEN y;
591 630003999 : if (!v || lx==2) return gerepileuptoleaf(av, x);
592 173933063 : ly = lx + v; /* result length */
593 173933063 : (void)new_chunk(ly); /* check that result fits */
594 173832027 : x += lx; y = (GEN)av;
595 1232175993 : for (i = 2; i<lx; i++) *--y = *--x;
596 700857925 : for (i = 0; i< v; i++) *--y = 0;
597 173832027 : y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
598 173978338 : return gc_const((pari_sp)y, y);
599 : }
600 :
601 : static long
602 2305331478 : get_Fl_threshold(ulong p, long mul, long mul2)
603 : {
604 2305331478 : return SMALL_ULONG(p) ? mul: mul2;
605 : }
606 :
607 : #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
608 : #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
609 : #define LLQUARTWORD(x) ((x) & QUARTMASK)
610 : #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
611 : #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
612 : #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
613 : INLINE long
614 8331200 : maxbitcoeffpol(ulong p, long n)
615 : {
616 8331200 : GEN z = muliu(sqru(p - 1), n);
617 8328271 : long b = expi(z) + 1;
618 : /* only do expensive bit-packing if it saves at least 1 limb */
619 8328986 : if (b <= BITS_IN_QUARTULONG)
620 : {
621 878237 : if (nbits2nlong(n*b) == (n + 3)>>2)
622 107333 : b = BITS_IN_QUARTULONG;
623 : }
624 7450749 : else if (b <= BITS_IN_HALFULONG)
625 : {
626 1544068 : if (nbits2nlong(n*b) == (n + 1)>>1)
627 5590 : b = BITS_IN_HALFULONG;
628 : }
629 : else
630 : {
631 5906681 : long l = lgefint(z) - 2;
632 5906681 : if (nbits2nlong(n*b) == n*l)
633 307363 : b = l*BITS_IN_LONG;
634 : }
635 8328807 : return b;
636 : }
637 :
638 : INLINE ulong
639 3373202348 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
640 : { /* Assume OK_ULONG*/
641 3373202348 : ulong p1 = 0;
642 : long i;
643 15971161894 : for (i=a; i<b; i++)
644 12597959546 : if (y[i])
645 : {
646 10585236904 : p1 += y[i] * x[-i];
647 10585236904 : if (p1 & HIGHBIT) p1 %= p;
648 : }
649 3373202348 : return p1 % p;
650 : }
651 :
652 : INLINE ulong
653 1151141937 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
654 : {
655 1151141937 : ulong p1 = 0;
656 : long i;
657 3616937000 : for (i=a; i<b; i++)
658 2466013540 : if (y[i])
659 2441825639 : p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
660 1150923460 : return p1;
661 : }
662 :
663 : /* assume nx >= ny > 0 */
664 : static GEN
665 341524159 : Flx_mulspec_basecase(GEN x, GEN y, ulong p, ulong pi, long nx, long ny)
666 : {
667 : long i,lz,nz;
668 : GEN z;
669 :
670 341524159 : lz = nx+ny+1; nz = lz-2;
671 341524159 : z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
672 341346855 : if (!pi)
673 : {
674 1139050841 : for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
675 727051038 : for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
676 887484473 : for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
677 : }
678 : else
679 : {
680 307250092 : for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
681 213880061 : for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
682 219050491 : for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
683 : }
684 341276722 : z -= 2; return Flx_renormalize(z, lz);
685 : }
686 :
687 : static GEN
688 12302 : int_to_Flx(GEN z, ulong p)
689 : {
690 12302 : long i, l = lgefint(z);
691 12302 : GEN x = cgetg(l, t_VECSMALL);
692 1060162 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
693 12299 : return Flx_renormalize(x, l);
694 : }
695 :
696 : INLINE GEN
697 10034 : Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
698 : {
699 10034 : GEN z=muliispec(a,b,na,nb);
700 10038 : return int_to_Flx(z,p);
701 : }
702 :
703 : static GEN
704 469541 : Flx_to_int_halfspec(GEN a, long na)
705 : {
706 : long j;
707 469541 : long n = (na+1)>>1UL;
708 469541 : GEN V = cgetipos(2+n);
709 : GEN w;
710 1378926 : for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
711 909385 : *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
712 469541 : if (j<na)
713 319481 : *w = a[j];
714 469541 : return V;
715 : }
716 :
717 : static GEN
718 506363 : int_to_Flx_half(GEN z, ulong p)
719 : {
720 : long i;
721 506363 : long lx = (lgefint(z)-2)*2+2;
722 506363 : GEN w, x = cgetg(lx, t_VECSMALL);
723 1909850 : for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
724 : {
725 1403487 : x[i] = LOWWORD((ulong)*w)%p;
726 1403487 : x[i+1] = HIGHWORD((ulong)*w)%p;
727 : }
728 506363 : return Flx_renormalize(x, lx);
729 : }
730 :
731 : static GEN
732 5454 : Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
733 : {
734 5454 : GEN A = Flx_to_int_halfspec(a,na);
735 5454 : GEN B = Flx_to_int_halfspec(b,nb);
736 5454 : GEN z = mulii(A,B);
737 5454 : return int_to_Flx_half(z,p);
738 : }
739 :
740 : static GEN
741 204445 : Flx_to_int_quartspec(GEN a, long na)
742 : {
743 : long j;
744 204445 : long n = (na+3)>>2UL;
745 204445 : GEN V = cgetipos(2+n);
746 : GEN w;
747 4377342 : for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
748 4172891 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
749 204451 : switch (na-j)
750 : {
751 116376 : case 3:
752 116376 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
753 116376 : break;
754 34469 : case 2:
755 34469 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
756 34469 : break;
757 27256 : case 1:
758 27256 : *w = a[j];
759 27256 : break;
760 26350 : case 0:
761 26350 : break;
762 : }
763 204451 : return V;
764 : }
765 :
766 : static GEN
767 107335 : int_to_Flx_quart(GEN z, ulong p)
768 : {
769 : long i;
770 107335 : long lx = (lgefint(z)-2)*4+2;
771 107335 : GEN w, x = cgetg(lx, t_VECSMALL);
772 4873496 : for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
773 : {
774 4766161 : x[i] = LLQUARTWORD((ulong)*w)%p;
775 4766161 : x[i+1] = HLQUARTWORD((ulong)*w)%p;
776 4766161 : x[i+2] = LHQUARTWORD((ulong)*w)%p;
777 4766161 : x[i+3] = HHQUARTWORD((ulong)*w)%p;
778 : }
779 107335 : return Flx_renormalize(x, lx);
780 : }
781 :
782 : static GEN
783 97112 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
784 : {
785 97112 : GEN A = Flx_to_int_quartspec(a,na);
786 97116 : GEN B = Flx_to_int_quartspec(b,nb);
787 97116 : GEN z = mulii(A,B);
788 97114 : return int_to_Flx_quart(z,p);
789 : }
790 :
791 : /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
792 : static GEN
793 582022 : Flx_eval2BILspec(GEN x, long k, long l)
794 : {
795 582022 : long i, lz = k*l, ki;
796 582022 : GEN pz = cgetipos(2+lz);
797 16364366 : for (i=0; i < lz; i++)
798 15782344 : *int_W(pz,i) = 0UL;
799 8473194 : for (i=0, ki=0; i<l; i++, ki+=k)
800 7891172 : *int_W(pz,ki) = x[i];
801 582022 : return int_normalize(pz,0);
802 : }
803 :
804 : static GEN
805 297997 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
806 : {
807 297997 : long i, offset, lm = lgefint(x)-2, l = d+3;
808 297997 : ulong pi = get_Fl_red(p);
809 297997 : GEN pol = cgetg(l, t_VECSMALL);
810 297997 : pol[1] = 0;
811 8007722 : for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
812 7709725 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
813 297997 : if (offset < lm)
814 225032 : pol[i+2] = (*int_W(x,offset)) % p;
815 297997 : return Flx_renormalize(pol,l);
816 : }
817 :
818 : static GEN
819 0 : Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
820 : {
821 0 : long i, offset, lm = lgefint(x)-2, l = d+3;
822 0 : ulong pi = get_Fl_red(p);
823 0 : GEN pol = cgetg(l, t_VECSMALL);
824 0 : pol[1] = 0;
825 0 : for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
826 0 : pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
827 0 : *int_W(x,offset), p, pi);
828 0 : if (offset+1 < lm)
829 0 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
830 0 : else if (offset < lm)
831 0 : pol[i+2] = (*int_W(x,offset)) % p;
832 0 : return Flx_renormalize(pol,l);
833 : }
834 :
835 : static GEN
836 295067 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
837 : {
838 295067 : return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
839 : }
840 :
841 : static GEN
842 283566 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
843 : {
844 283566 : pari_sp av = avma;
845 283566 : GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
846 283566 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
847 : }
848 :
849 : static GEN
850 20719281 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
851 : GEN y;
852 : long i;
853 20719281 : if (l == 0)
854 3427632 : return gen_0;
855 17291649 : y = cgetg(l + 1, t_VECSMALL);
856 813015279 : for(i = 1; i <= l; i++)
857 795727255 : y[i] = x[l - i];
858 17288024 : return nv_fromdigits_2k(y, b);
859 : }
860 :
861 : /* assume b < BITS_IN_LONG */
862 : static GEN
863 5645174 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
864 5645174 : GEN v = binary_2k_nv(z, b), x;
865 5645232 : long i, l = lg(v) + 1;
866 5645232 : x = cgetg(l, t_VECSMALL);
867 621255877 : for (i = 2; i < l; i++)
868 615610565 : x[i] = v[l - i] % p;
869 5645312 : return Flx_renormalize(x, l);
870 : }
871 :
872 : static GEN
873 5538843 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
874 5538843 : GEN v = binary_2k(z, b), x, y;
875 5537998 : long i, l = lg(v) + 1, ly;
876 5537998 : x = cgetg(l, t_VECSMALL);
877 232951589 : for (i = 2; i < l; i++) {
878 227415831 : y = gel(v, l - i);
879 227415831 : ly = lgefint(y);
880 227415831 : switch (ly) {
881 6286357 : case 2: x[i] = 0; break;
882 29285284 : case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
883 175973281 : case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
884 31741440 : case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
885 15870909 : *int_W_lg(y, 0, ly), p, pi); break;
886 0 : default: x[i] = umodiu(gel(v, l - i), p);
887 : }
888 : }
889 5535758 : return Flx_renormalize(x, l);
890 : }
891 :
892 : static GEN
893 7223628 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
894 : {
895 : GEN C, D;
896 7223628 : pari_sp av = avma;
897 7223628 : A = kron_pack_Flx_spec_bits(A, b, lA);
898 7229754 : B = kron_pack_Flx_spec_bits(B, b, lB);
899 7229804 : C = gerepileuptoint(av, mulii(A, B));
900 7228522 : if (b < BITS_IN_LONG)
901 2063408 : D = kron_unpack_Flx_bits_narrow(C, b, p);
902 : else
903 : {
904 5165114 : ulong pi = get_Fl_red(p);
905 5163925 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
906 : }
907 7226003 : return D;
908 : }
909 :
910 : static GEN
911 683859 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
912 : {
913 : GEN C, D;
914 683859 : A = kron_pack_Flx_spec_bits(A, b, lA);
915 683939 : C = sqri(A);
916 683965 : if (b < BITS_IN_LONG)
917 475779 : D = kron_unpack_Flx_bits_narrow(C, b, p);
918 : else
919 : {
920 208186 : ulong pi = get_Fl_red(p);
921 208181 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
922 : }
923 683922 : return D;
924 : }
925 :
926 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
927 : * b+2 were sent instead. na, nb = number of terms of a, b.
928 : * Only c, c0, c1, c2 are genuine GEN.
929 : */
930 : static GEN
931 378864027 : Flx_mulspec(GEN a, GEN b, ulong p, ulong pi, long na, long nb)
932 : {
933 : GEN a0,c,c0;
934 378864027 : long n0, n0a, i, v = 0;
935 : pari_sp av;
936 :
937 483178517 : while (na && !a[0]) { a++; na--; v++; }
938 563311277 : while (nb && !b[0]) { b++; nb--; v++; }
939 378864027 : if (na < nb) swapspec(a,b, na,nb);
940 378864027 : if (!nb) return pol0_Flx(0);
941 :
942 350703735 : av = avma;
943 350703735 : if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
944 : {
945 7623576 : long m = maxbitcoeffpol(p,nb);
946 7619515 : switch (m)
947 : {
948 97112 : case BITS_IN_QUARTULONG:
949 97112 : return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
950 5454 : case BITS_IN_HALFULONG:
951 5454 : return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
952 10034 : case BITS_IN_LONG:
953 10034 : return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
954 283566 : case 2*BITS_IN_LONG:
955 283566 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
956 0 : case 3*BITS_IN_LONG:
957 0 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
958 7223349 : default:
959 7223349 : return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
960 : }
961 : }
962 343248121 : if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
963 341462162 : return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,pi,na,nb), v);
964 1800427 : i=(na>>1); n0=na-i; na=i;
965 1800427 : a0=a+n0; n0a=n0;
966 2566113 : while (n0a && !a[n0a-1]) n0a--;
967 :
968 1800427 : if (nb > n0)
969 : {
970 : GEN b0,c1,c2;
971 : long n0b;
972 :
973 1746563 : nb -= n0; b0 = b+n0; n0b = n0;
974 2826052 : while (n0b && !b[n0b-1]) n0b--;
975 1746563 : c = Flx_mulspec(a,b,p,pi,n0a,n0b);
976 1746563 : c0 = Flx_mulspec(a0,b0,p,pi,na,nb);
977 :
978 1746563 : c2 = Flx_addspec(a0,a,p,na,n0a);
979 1746563 : c1 = Flx_addspec(b0,b,p,nb,n0b);
980 :
981 1746563 : c1 = Flx_mul_pre(c1,c2,p,pi);
982 1746563 : c2 = Flx_add(c0,c,p);
983 :
984 1746563 : c2 = Flx_neg_inplace(c2,p);
985 1746563 : c2 = Flx_add(c1,c2,p);
986 1746563 : c0 = Flx_addshift(c0,c2 ,p, n0);
987 : }
988 : else
989 : {
990 53864 : c = Flx_mulspec(a,b,p,pi,n0a,nb);
991 53864 : c0 = Flx_mulspec(a0,b,p,pi,na,nb);
992 : }
993 1800427 : c0 = Flx_addshift(c0,c,p,n0);
994 1800427 : return Flx_shiftip(av,c0, v);
995 : }
996 :
997 : GEN
998 373268010 : Flx_mul_pre(GEN x, GEN y, ulong p, ulong pi)
999 : {
1000 373268010 : GEN z = Flx_mulspec(x+2,y+2,p, pi, lgpol(x),lgpol(y));
1001 373427297 : z[1] = x[1]; return z;
1002 : }
1003 : GEN
1004 27690703 : Flx_mul(GEN x, GEN y, ulong p)
1005 27690703 : { return Flx_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1006 :
1007 : static GEN
1008 279705964 : Flx_sqrspec_basecase(GEN x, ulong p, ulong pi, long nx)
1009 : {
1010 : long i, lz, nz;
1011 : ulong p1;
1012 : GEN z;
1013 :
1014 279705964 : if (!nx) return pol0_Flx(0);
1015 279705964 : lz = (nx << 1) + 1, nz = lz-2;
1016 279705964 : z = cgetg(lz, t_VECSMALL) + 2;
1017 279187503 : if (!pi)
1018 : {
1019 214016457 : z[0] = x[0]*x[0]%p;
1020 915267765 : for (i=1; i<nx; i++)
1021 : {
1022 701358582 : p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1023 701251308 : p1 <<= 1;
1024 701251308 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1025 701251308 : z[i] = p1 % p;
1026 : }
1027 919368008 : for ( ; i<nz; i++)
1028 : {
1029 704814216 : p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1030 705458825 : p1 <<= 1;
1031 705458825 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1032 705458825 : z[i] = p1 % p;
1033 : }
1034 : }
1035 : else
1036 : {
1037 65171046 : z[0] = Fl_sqr_pre(x[0], p, pi);
1038 409164954 : for (i=1; i<nx; i++)
1039 : {
1040 344112274 : p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1041 344364832 : p1 = Fl_add(p1, p1, p);
1042 343917598 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1043 343925196 : z[i] = p1;
1044 : }
1045 409295672 : for ( ; i<nz; i++)
1046 : {
1047 344127880 : p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1048 344875282 : p1 = Fl_add(p1, p1, p);
1049 344507360 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1050 344242992 : z[i] = p1;
1051 : }
1052 : }
1053 279721584 : z -= 2; return Flx_renormalize(z, lz);
1054 : }
1055 :
1056 : static GEN
1057 2262 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
1058 : {
1059 2262 : GEN z=sqrispec(a,na);
1060 2265 : return int_to_Flx(z,p);
1061 : }
1062 :
1063 : static GEN
1064 136 : Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
1065 : {
1066 136 : GEN z = sqri(Flx_to_int_halfspec(a,na));
1067 136 : return int_to_Flx_half(z,p);
1068 : }
1069 :
1070 : static GEN
1071 10221 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1072 : {
1073 10221 : GEN z = sqri(Flx_to_int_quartspec(a,na));
1074 10221 : return int_to_Flx_quart(z,p);
1075 : }
1076 :
1077 : static GEN
1078 11501 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1079 : {
1080 11501 : pari_sp av = avma;
1081 11501 : GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1082 11501 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1083 : }
1084 :
1085 : static GEN
1086 279893416 : Flx_sqrspec(GEN a, ulong p, ulong pi, long na)
1087 : {
1088 : GEN a0, c, c0;
1089 279893416 : long n0, n0a, i, v = 0, m;
1090 : pari_sp av;
1091 :
1092 400892262 : while (na && !a[0]) { a++; na--; v += 2; }
1093 279893416 : if (!na) return pol0_Flx(0);
1094 :
1095 279647225 : av = avma;
1096 279647225 : if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1097 : {
1098 707939 : m = maxbitcoeffpol(p,na);
1099 707977 : switch(m)
1100 : {
1101 10221 : case BITS_IN_QUARTULONG:
1102 10221 : return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
1103 136 : case BITS_IN_HALFULONG:
1104 136 : return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
1105 2262 : case BITS_IN_LONG:
1106 2262 : return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1107 11501 : case 2*BITS_IN_LONG:
1108 11501 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
1109 0 : case 3*BITS_IN_LONG:
1110 0 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
1111 683857 : default:
1112 683857 : return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1113 : }
1114 : }
1115 279398473 : if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1116 279335488 : return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,pi,na), v);
1117 57585 : i=(na>>1); n0=na-i; na=i;
1118 57585 : a0=a+n0; n0a=n0;
1119 72314 : while (n0a && !a[n0a-1]) n0a--;
1120 :
1121 57585 : c = Flx_sqrspec(a,p,pi,n0a);
1122 57585 : c0= Flx_sqrspec(a0,p,pi,na);
1123 57585 : if (p == 2) n0 *= 2;
1124 : else
1125 : {
1126 57566 : GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1127 57566 : t = Flx_sqr_pre(t,p,pi);
1128 57566 : c1= Flx_add(c0,c, p);
1129 57566 : c1= Flx_sub(t, c1, p);
1130 57566 : c0 = Flx_addshift(c0,c1,p,n0);
1131 : }
1132 57585 : c0 = Flx_addshift(c0,c,p,n0);
1133 57585 : return Flx_shiftip(av,c0,v);
1134 : }
1135 :
1136 : GEN
1137 279658814 : Flx_sqr_pre(GEN x, ulong p, ulong pi)
1138 : {
1139 279658814 : GEN z = Flx_sqrspec(x+2,p, pi, lgpol(x));
1140 280908583 : z[1] = x[1]; return z;
1141 : }
1142 : GEN
1143 356367 : Flx_sqr(GEN x, ulong p)
1144 356367 : { return Flx_sqr_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1145 :
1146 : GEN
1147 8031 : Flx_powu_pre(GEN x, ulong n, ulong p, ulong pi)
1148 : {
1149 8031 : GEN y = pol1_Flx(x[1]), z;
1150 : ulong m;
1151 8030 : if (n == 0) return y;
1152 8030 : m = n; z = x;
1153 : for (;;)
1154 : {
1155 30973 : if (m&1UL) y = Flx_mul_pre(y,z, p, pi);
1156 30970 : m >>= 1; if (!m) return y;
1157 22941 : z = Flx_sqr_pre(z, p, pi);
1158 : }
1159 : }
1160 : GEN
1161 0 : Flx_powu(GEN x, ulong n, ulong p)
1162 : {
1163 0 : if (n == 0) return pol1_Flx(x[1]);
1164 0 : return Flx_powu_pre(x, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
1165 : }
1166 :
1167 : GEN
1168 14222 : Flx_halve(GEN y, ulong p)
1169 : {
1170 : GEN z;
1171 : long i, l;
1172 14222 : z = cgetg_copy(y, &l); z[1] = y[1];
1173 59732 : for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
1174 14222 : return z;
1175 : }
1176 :
1177 : static GEN
1178 7130920 : Flx_recipspec(GEN x, long l, long n)
1179 : {
1180 : long i;
1181 7130920 : GEN z=cgetg(n+2,t_VECSMALL)+2;
1182 115930915 : for(i=0; i<l; i++)
1183 108801516 : z[n-i-1] = x[i];
1184 15598279 : for( ; i<n; i++)
1185 8468880 : z[n-i-1] = 0;
1186 7129399 : return Flx_renormalize(z-2,n+2);
1187 : }
1188 :
1189 : GEN
1190 0 : Flx_recip(GEN x)
1191 : {
1192 0 : GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
1193 0 : z[1]=x[1];
1194 0 : return z;
1195 : }
1196 :
1197 : /* Return h^degpol(P) P(x / h) */
1198 : GEN
1199 1117 : Flx_rescale(GEN P, ulong h, ulong p)
1200 : {
1201 1117 : long i, l = lg(P);
1202 1117 : GEN Q = cgetg(l,t_VECSMALL);
1203 1117 : ulong hi = h;
1204 1117 : Q[l-1] = P[l-1];
1205 12538 : for (i=l-2; i>=2; i--)
1206 : {
1207 12538 : Q[i] = Fl_mul(P[i], hi, p);
1208 12538 : if (i == 2) break;
1209 11421 : hi = Fl_mul(hi,h, p);
1210 : }
1211 1117 : Q[1] = P[1]; return Q;
1212 : }
1213 :
1214 : /* x/polrecip(P)+O(x^n); allow pi = 0 */
1215 : static GEN
1216 134232 : Flx_invBarrett_basecase(GEN T, ulong p, ulong pi)
1217 : {
1218 134232 : long i, l=lg(T)-1, lr=l-1, k;
1219 134232 : GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1220 134232 : r[2] = 1;
1221 134232 : if (!pi)
1222 764017 : for (i=3;i<lr;i++)
1223 : {
1224 757026 : ulong u = uel(T, l-i+2);
1225 45370123 : for (k=3; k<i; k++)
1226 44613097 : { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1227 757026 : r[i] = Fl_neg(u % p, p);
1228 : }
1229 : else
1230 2109686 : for (i=3;i<lr;i++)
1231 : {
1232 1982449 : ulong u = Fl_neg(uel(T,l-i+2), p);
1233 59521945 : for (k=3; k<i; k++)
1234 : {
1235 57539499 : ulong t = Fl_neg(uel(T,l-i+k), p);
1236 57539499 : u = Fl_addmul_pre(u, t, uel(r,k), p, pi);
1237 : }
1238 1982446 : r[i] = u;
1239 : }
1240 134228 : return Flx_renormalize(r,lr);
1241 : }
1242 :
1243 : /* Return new lgpol */
1244 : static long
1245 2134146 : Flx_lgrenormalizespec(GEN x, long lx)
1246 : {
1247 : long i;
1248 7440078 : for (i = lx-1; i>=0; i--)
1249 7439246 : if (x[i]) break;
1250 2134146 : return i+1;
1251 : }
1252 : /* allow pi = 0 */
1253 : static GEN
1254 23115 : Flx_invBarrett_Newton(GEN T, ulong p, ulong pi)
1255 : {
1256 23115 : long nold, lx, lz, lq, l = degpol(T), lQ;
1257 23115 : GEN q, y, z, x = zero_zv(l+1) + 2;
1258 23116 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1259 : pari_sp av;
1260 :
1261 23116 : y = T+2;
1262 23116 : q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1263 23116 : av = avma;
1264 : /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1265 :
1266 : /* initialize */
1267 23116 : x[0] = Fl_inv(q[0], p);
1268 23115 : if (lQ>1 && q[1])
1269 5109 : {
1270 5109 : ulong u = q[1];
1271 5109 : if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1272 5109 : x[1] = p - u; lx = 2;
1273 : }
1274 : else
1275 18006 : lx = 1;
1276 23115 : nold = 1;
1277 158704 : for (; mask > 1; set_avma(av))
1278 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1279 135593 : long i, lnew, nnew = nold << 1;
1280 :
1281 135593 : if (mask & 1) nnew--;
1282 135593 : mask >>= 1;
1283 :
1284 135593 : lnew = nnew + 1;
1285 135593 : lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1286 135602 : z = Flx_mulspec(x, q, p, pi, lx, lq); /* FIXME: high product */
1287 135589 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1288 135588 : z += 2;
1289 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1290 290673 : for (i = nold; i < lz; i++) if (z[i]) break;
1291 135588 : nold = nnew;
1292 135588 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1293 :
1294 : /* z + i represents (x*q - 1) / t^i */
1295 100751 : lz = Flx_lgrenormalizespec (z+i, lz-i);
1296 100754 : z = Flx_mulspec(x, z+i, p, pi, lx, lz); /* FIXME: low product */
1297 100752 : lz = lgpol(z); z += 2;
1298 100752 : if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1299 :
1300 100752 : lx = lz+ i;
1301 100752 : y = x + i; /* x -= z * t^i, in place */
1302 915345 : for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1303 : }
1304 23116 : x -= 2; setlg(x, lx + 2); x[1] = T[1];
1305 23115 : return x;
1306 : }
1307 :
1308 : /* allow pi = 0 */
1309 : static GEN
1310 158659 : Flx_invBarrett_pre(GEN T, ulong p, ulong pi)
1311 : {
1312 158659 : pari_sp ltop = avma;
1313 158659 : long l = lgpol(T);
1314 : GEN r;
1315 158659 : if (l < 3) return pol0_Flx(T[1]);
1316 157347 : if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1317 : {
1318 134232 : ulong c = T[l+1];
1319 134232 : if (c != 1)
1320 : {
1321 98118 : ulong ci = Fl_inv(c,p);
1322 98118 : T = Flx_Fl_mul_pre(T, ci, p, pi);
1323 98118 : r = Flx_invBarrett_basecase(T, p, pi);
1324 98118 : r = Flx_Fl_mul_pre(r, ci, p, pi);
1325 : }
1326 : else
1327 36114 : r = Flx_invBarrett_basecase(T, p, pi);
1328 : }
1329 : else
1330 23115 : r = Flx_invBarrett_Newton(T, p, pi);
1331 157347 : return gerepileuptoleaf(ltop, r);
1332 : }
1333 : GEN
1334 0 : Flx_invBarrett(GEN T, ulong p)
1335 0 : { return Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1336 :
1337 : /* allow pi = 0 */
1338 : GEN
1339 98501199 : Flx_get_red_pre(GEN T, ulong p, ulong pi)
1340 : {
1341 98501199 : if (typ(T)!=t_VECSMALL
1342 98464805 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1343 : Flx_BARRETT2_LIMIT))
1344 98482137 : return T;
1345 7611 : retmkvec2(Flx_invBarrett_pre(T, p, pi),T);
1346 : }
1347 : GEN
1348 14256059 : Flx_get_red(GEN T, ulong p)
1349 : {
1350 14256059 : if (typ(T)!=t_VECSMALL
1351 14255960 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1352 : Flx_BARRETT2_LIMIT))
1353 14250331 : return T;
1354 5194 : retmkvec2(Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)),T);
1355 : }
1356 :
1357 : /* separate from Flx_divrem for maximal speed. */
1358 : static GEN
1359 790083783 : Flx_rem_basecase(GEN x, GEN y, ulong p, ulong pi)
1360 : {
1361 : pari_sp av;
1362 : GEN z, c;
1363 : long dx,dy,dy1,dz,i,j;
1364 : ulong p1,inv;
1365 790083783 : long vs=x[1];
1366 :
1367 790083783 : dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1368 754604901 : dx = degpol(x);
1369 754498125 : dz = dx-dy; if (dz < 0) return Flx_copy(x);
1370 754498125 : x += 2; y += 2;
1371 754498125 : inv = y[dy];
1372 754498125 : if (inv != 1UL) inv = Fl_inv(inv,p);
1373 909630150 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1374 :
1375 756052435 : c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1376 754800984 : z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1377 :
1378 753220581 : if (!pi)
1379 : {
1380 480548166 : z[dz] = (inv*x[dx]) % p;
1381 1802675882 : for (i=dx-1; i>=dy; --i)
1382 : {
1383 1322127716 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1384 10458859436 : for (j=i-dy1; j<=i && j<=dz; j++)
1385 : {
1386 9136731720 : p1 += z[j]*y[i-j];
1387 9136731720 : if (p1 & HIGHBIT) p1 %= p;
1388 : }
1389 1322127716 : p1 %= p;
1390 1322127716 : z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1391 : }
1392 3278379878 : for (i=0; i<dy; i++)
1393 : {
1394 2797868787 : p1 = z[0]*y[i];
1395 14449166507 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1396 : {
1397 11651297720 : p1 += z[j]*y[i-j];
1398 11651297720 : if (p1 & HIGHBIT) p1 %= p;
1399 : }
1400 2798129614 : c[i] = Fl_sub(x[i], p1%p, p);
1401 : }
1402 : }
1403 : else
1404 : {
1405 272672415 : z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1406 825386164 : for (i=dx-1; i>=dy; --i)
1407 : {
1408 552419279 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1409 2323922610 : for (j=i-dy1; j<=i && j<=dz; j++)
1410 1770756422 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1411 553166188 : z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1412 : }
1413 2002559716 : for (i=0; i<dy; i++)
1414 : {
1415 1730673772 : p1 = Fl_mul_pre(z[0],y[i],p,pi);
1416 4643840557 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1417 2904754574 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1418 1716608844 : c[i] = Fl_sub(x[i], p1, p);
1419 : }
1420 : }
1421 919449244 : i = dy-1; while (i>=0 && !c[i]) i--;
1422 752397035 : set_avma(av); return Flx_renormalize(c-2, i+3);
1423 : }
1424 :
1425 : /* as FpX_divrem but working only on ulong types.
1426 : * if relevant, *pr is the last object on stack */
1427 : static GEN
1428 61917578 : Flx_divrem_basecase(GEN x, GEN y, ulong p, ulong pi, GEN *pr)
1429 : {
1430 : GEN z,q,c;
1431 : long dx,dy,dy1,dz,i,j;
1432 : ulong p1,inv;
1433 61917578 : long sv=x[1];
1434 :
1435 61917578 : dy = degpol(y);
1436 61916405 : if (dy<0) pari_err_INV("Flx_divrem",y);
1437 61916519 : if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p, pi);
1438 61916121 : if (!dy)
1439 : {
1440 7223565 : if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1441 7223544 : if (y[2] == 1UL) return Flx_copy(x);
1442 5210562 : return Flx_Fl_mul_pre(x, Fl_inv(y[2], p), p, pi);
1443 : }
1444 54692556 : dx = degpol(x);
1445 54695250 : dz = dx-dy;
1446 54695250 : if (dz < 0)
1447 : {
1448 1028271 : q = pol0_Flx(sv);
1449 1028265 : if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1450 1028265 : return q;
1451 : }
1452 53666979 : x += 2;
1453 53666979 : y += 2;
1454 53666979 : z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1455 53665707 : inv = uel(y, dy);
1456 53665707 : if (inv != 1UL) inv = Fl_inv(inv,p);
1457 79016069 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1458 :
1459 53668635 : if (SMALL_ULONG(p))
1460 : {
1461 51792091 : z[dz] = (inv*x[dx]) % p;
1462 131472381 : for (i=dx-1; i>=dy; --i)
1463 : {
1464 79680290 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1465 257736550 : for (j=i-dy1; j<=i && j<=dz; j++)
1466 : {
1467 178056260 : p1 += z[j]*y[i-j];
1468 178056260 : if (p1 & HIGHBIT) p1 %= p;
1469 : }
1470 79680290 : p1 %= p;
1471 79680290 : z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1472 : }
1473 : }
1474 : else
1475 : {
1476 1876544 : z[dz] = Fl_mul(inv, x[dx], p);
1477 9245165 : for (i=dx-1; i>=dy; --i)
1478 : { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1479 7368525 : p1 = p - uel(x,i);
1480 26361706 : for (j=i-dy1; j<=i && j<=dz; j++)
1481 18993178 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1482 7368528 : z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1483 : }
1484 : }
1485 53668731 : q = Flx_renormalize(z-2, dz+3);
1486 53667625 : if (!pr) return q;
1487 :
1488 26526674 : c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1489 26528710 : if (SMALL_ULONG(p))
1490 : {
1491 225707414 : for (i=0; i<dy; i++)
1492 : {
1493 200815118 : p1 = (ulong)z[0]*y[i];
1494 470788592 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1495 : {
1496 269973474 : p1 += (ulong)z[j]*y[i-j];
1497 269973474 : if (p1 & HIGHBIT) p1 %= p;
1498 : }
1499 200814760 : c[i] = Fl_sub(x[i], p1%p, p);
1500 : }
1501 : }
1502 : else
1503 : {
1504 16026511 : for (i=0; i<dy; i++)
1505 : {
1506 14390785 : p1 = Fl_mul(z[0],y[i],p);
1507 50219813 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1508 35829026 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1509 14390785 : c[i] = Fl_sub(x[i], p1, p);
1510 : }
1511 : }
1512 35665812 : i=dy-1; while (i>=0 && !c[i]) i--;
1513 26528022 : c = Flx_renormalize(c-2, i+3);
1514 26528665 : if (pr == ONLY_DIVIDES)
1515 427 : { if (lg(c) != 2) return NULL; }
1516 : else
1517 26528238 : *pr = c;
1518 26528525 : return q;
1519 : }
1520 :
1521 : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1522 : * and mg is the Barrett inverse of T. */
1523 : static GEN
1524 906243 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1525 : {
1526 : GEN q, r;
1527 906243 : long lt = degpol(T); /*We discard the leading term*/
1528 : long ld, lm, lT, lmg;
1529 906212 : ld = l-lt;
1530 906212 : lm = minss(ld, lgpol(mg));
1531 906580 : lT = Flx_lgrenormalizespec(T+2,lt);
1532 906714 : lmg = Flx_lgrenormalizespec(mg+2,lm);
1533 906537 : q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1534 905760 : q = Flx_mulspec(q+2,mg+2,p,pi,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1535 906488 : q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1536 905671 : if (!pr) return q;
1537 897980 : r = Flx_mulspec(q+2,T+2,p,pi,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1538 898793 : r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1539 898589 : if (pr == ONLY_REM) return r;
1540 427909 : *pr = r; return q;
1541 : }
1542 :
1543 : static GEN
1544 605791 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1545 : {
1546 605791 : GEN q = NULL, r = Flx_copy(x);
1547 605818 : long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1548 : long i;
1549 605815 : if (l <= lt)
1550 : {
1551 0 : if (pr == ONLY_REM) return Flx_copy(x);
1552 0 : if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1553 0 : if (pr) *pr = Flx_copy(x);
1554 0 : return pol0_Flx(v);
1555 : }
1556 605815 : if (lt <= 1)
1557 1312 : return Flx_divrem_basecase(x,T,p,pi,pr);
1558 604503 : if (pr != ONLY_REM && l>lm)
1559 28917 : { q = zero_zv(l-lt+1); q[1] = T[1]; }
1560 907830 : while (l>lm)
1561 : {
1562 303382 : GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,pi,&zr);
1563 303392 : long lz = lgpol(zr);
1564 303327 : if (pr != ONLY_REM)
1565 : {
1566 58007 : long lq = lgpol(zq);
1567 872592 : for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1568 : }
1569 4392088 : for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1570 303327 : l = l-lm+lz;
1571 : }
1572 604448 : if (pr == ONLY_REM)
1573 : {
1574 470743 : if (l > lt)
1575 470701 : r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi,ONLY_REM);
1576 : else
1577 42 : r = Flx_renormalize(r, l+2);
1578 470722 : r[1] = v; return r;
1579 : }
1580 133705 : if (l > lt)
1581 : {
1582 132187 : GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi, pr ? &r: NULL);
1583 132187 : if (!q) q = zq;
1584 : else
1585 : {
1586 27343 : long lq = lgpol(zq);
1587 158719 : for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1588 : }
1589 : }
1590 1518 : else if (pr)
1591 1535 : r = Flx_renormalize(r, l+2);
1592 133705 : q[1] = v; q = Flx_renormalize(q, lg(q));
1593 133761 : if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1594 133761 : if (pr) { r[1] = v; *pr = r; }
1595 133761 : return q;
1596 : }
1597 :
1598 : /* allow pi = 0 (SMALL_ULONG) */
1599 : GEN
1600 79362585 : Flx_divrem_pre(GEN x, GEN T, ulong p, ulong pi, GEN *pr)
1601 : {
1602 : GEN B, y;
1603 : long dy, dx, d;
1604 79362585 : if (pr==ONLY_REM) return Flx_rem_pre(x, T, p, pi);
1605 62042896 : y = get_Flx_red(T, &B);
1606 62053618 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1607 62050000 : if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1608 61915982 : return Flx_divrem_basecase(x,y,p,pi,pr);
1609 : else
1610 : {
1611 134675 : pari_sp av = avma;
1612 134675 : GEN mg = B? B: Flx_invBarrett_pre(y, p, pi);
1613 134675 : GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pi,pr);
1614 134675 : if (!q1) return gc_NULL(av);
1615 134675 : if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
1616 126378 : return gc_all(av, 2, &q1, pr);
1617 : }
1618 : }
1619 : GEN
1620 30293976 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1621 30293976 : { return Flx_divrem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p), pr); }
1622 :
1623 : GEN
1624 913520350 : Flx_rem_pre(GEN x, GEN T, ulong p, ulong pi)
1625 : {
1626 913520350 : GEN B, y = get_Flx_red(T, &B);
1627 913401938 : long d = degpol(x) - degpol(y);
1628 913092527 : if (d < 0) return Flx_copy(x);
1629 790677287 : if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1630 790105191 : return Flx_rem_basecase(x,y,p, pi);
1631 : else
1632 : {
1633 471117 : pari_sp av=avma;
1634 471117 : GEN mg = B ? B: Flx_invBarrett_pre(y, p, pi);
1635 471117 : GEN r = Flx_divrem_Barrett(x, mg, y, p, pi, ONLY_REM);
1636 471120 : return gerepileuptoleaf(av, r);
1637 : }
1638 : }
1639 : GEN
1640 41844239 : Flx_rem(GEN x, GEN T, ulong p)
1641 41844239 : { return Flx_rem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1642 :
1643 : /* reduce T mod (X^n - 1, p). Shallow function */
1644 : GEN
1645 5107324 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1646 : {
1647 5107324 : long i, j, L = lg(T), l = n+2;
1648 : GEN S;
1649 5107324 : if (L <= l || n & ~LGBITS) return T;
1650 3450 : S = cgetg(l, t_VECSMALL);
1651 3450 : S[1] = T[1];
1652 14013 : for (i = 2; i < l; i++) S[i] = T[i];
1653 9420 : for (j = 2; i < L; i++) {
1654 5970 : S[j] = Fl_add(S[j], T[i], p);
1655 5970 : if (++j == l) j = 2;
1656 : }
1657 3450 : return Flx_renormalize(S, l);
1658 : }
1659 : /* reduce T mod (X^n + 1, p). Shallow function */
1660 : GEN
1661 30303 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
1662 : {
1663 30303 : long i, j, L = lg(T), l = n+2;
1664 : GEN S;
1665 30303 : if (L <= l || n & ~LGBITS) return T;
1666 2682 : S = cgetg(l, t_VECSMALL);
1667 2682 : S[1] = T[1];
1668 11347 : for (i = 2; i < l; i++) S[i] = T[i];
1669 6974 : for (j = 2; i < L; i++) {
1670 4292 : S[j] = Fl_sub(S[j], T[i], p);
1671 4292 : if (++j == l) j = 2;
1672 : }
1673 2682 : return Flx_renormalize(S, l);
1674 : }
1675 :
1676 : struct _Flxq {
1677 : GEN aut, T;
1678 : ulong p, pi;
1679 : };
1680 : /* allow pi = 0 */
1681 : static void
1682 71100339 : set_Flxq_pre(struct _Flxq *D, GEN T, ulong p, ulong pi)
1683 : {
1684 71100339 : D->p = p;
1685 71100339 : D->pi = pi;
1686 71100339 : D->T = Flx_get_red_pre(T, p, pi);
1687 71095659 : }
1688 : static void
1689 68965 : set_Flxq(struct _Flxq *D, GEN T, ulong p)
1690 68965 : { set_Flxq_pre(D, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1691 :
1692 : static GEN
1693 0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1694 : {
1695 0 : struct _Flxq *D = (struct _Flxq*) E;
1696 0 : return Flx_divrem_pre(x, y, D->p, D->pi, r);
1697 : }
1698 : static GEN
1699 389822 : _Flx_add(void * E, GEN x, GEN y) {
1700 389822 : struct _Flxq *D = (struct _Flxq*) E;
1701 389822 : return Flx_add(x, y, D->p);
1702 : }
1703 : static GEN
1704 10481032 : _Flx_mul(void *E, GEN x, GEN y) {
1705 10481032 : struct _Flxq *D = (struct _Flxq*) E;
1706 10481032 : return Flx_mul_pre(x, y, D->p, D->pi);
1707 : }
1708 : static GEN
1709 0 : _Flx_sqr(void *E, GEN x) {
1710 0 : struct _Flxq *D = (struct _Flxq*) E;
1711 0 : return Flx_sqr_pre(x, D->p, D->pi);
1712 : }
1713 :
1714 : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1715 :
1716 : GEN
1717 0 : Flx_digits(GEN x, GEN T, ulong p)
1718 : {
1719 : struct _Flxq D;
1720 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
1721 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1722 0 : return gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1723 : }
1724 :
1725 : GEN
1726 0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1727 : {
1728 : struct _Flxq D;
1729 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1730 0 : return gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1731 : }
1732 :
1733 : long
1734 4160349 : Flx_val(GEN x)
1735 : {
1736 4160349 : long i, l=lg(x);
1737 4160349 : if (l==2) return LONG_MAX;
1738 4169279 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1739 4160349 : return i-2;
1740 : }
1741 : long
1742 26307797 : Flx_valrem(GEN x, GEN *Z)
1743 : {
1744 26307797 : long v, i, l=lg(x);
1745 : GEN y;
1746 26307797 : if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1747 28490327 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1748 26307797 : v = i-2;
1749 26307797 : if (v == 0) { *Z = x; return 0; }
1750 1028129 : l -= v;
1751 1028129 : y = cgetg(l, t_VECSMALL); y[1] = x[1];
1752 2636895 : for (i=2; i<l; i++) y[i] = x[i+v];
1753 1029110 : *Z = y; return v;
1754 : }
1755 :
1756 : GEN
1757 21160810 : Flx_deriv(GEN z, ulong p)
1758 : {
1759 21160810 : long i,l = lg(z)-1;
1760 : GEN x;
1761 21160810 : if (l < 2) l = 2;
1762 21160810 : x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1763 21159355 : if (HIGHWORD(l | p))
1764 57457923 : for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1765 : else
1766 85425428 : for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1767 21160459 : return Flx_renormalize(x,l);
1768 : }
1769 :
1770 : static GEN
1771 422528 : Flx_integXn(GEN x, long n, ulong p)
1772 : {
1773 422528 : long i, lx = lg(x);
1774 : GEN y;
1775 422528 : if (lx == 2) return Flx_copy(x);
1776 412714 : y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1777 2096167 : for (i=2; i<lx; i++)
1778 : {
1779 1682978 : ulong xi = uel(x,i);
1780 1682978 : if (xi == 0)
1781 13345 : uel(y,i) = 0;
1782 : else
1783 : {
1784 1669633 : ulong j = n+i-1;
1785 1669633 : ulong d = ugcd(j, xi);
1786 1669587 : if (d==1)
1787 1018121 : uel(y,i) = Fl_div(xi, j, p);
1788 : else
1789 651466 : uel(y,i) = Fl_div(xi/d, j/d, p);
1790 : }
1791 : }
1792 413189 : return Flx_renormalize(y, lx);;
1793 : }
1794 :
1795 : GEN
1796 0 : Flx_integ(GEN x, ulong p)
1797 : {
1798 0 : long i, lx = lg(x);
1799 : GEN y;
1800 0 : if (lx == 2) return Flx_copy(x);
1801 0 : y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1802 0 : uel(y,2) = 0;
1803 0 : for (i=3; i<=lx; i++)
1804 0 : uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1805 0 : return Flx_renormalize(y, lx+1);;
1806 : }
1807 :
1808 : /* assume p prime */
1809 : GEN
1810 13482 : Flx_diff1(GEN P, ulong p)
1811 : {
1812 13482 : return Flx_sub(Flx_translate1(P, p), P, p);
1813 : }
1814 :
1815 : GEN
1816 420111 : Flx_deflate(GEN x0, long d)
1817 : {
1818 : GEN z, y, x;
1819 420111 : long i,id, dy, dx = degpol(x0);
1820 420110 : if (d == 1 || dx <= 0) return Flx_copy(x0);
1821 356612 : dy = dx/d;
1822 356612 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1823 356610 : z = y + 2;
1824 356610 : x = x0+ 2;
1825 1159807 : for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1826 356610 : return y;
1827 : }
1828 :
1829 : GEN
1830 158464 : Flx_inflate(GEN x0, long d)
1831 : {
1832 158464 : long i, id, dy, dx = degpol(x0);
1833 158462 : GEN x = x0 + 2, z, y;
1834 158462 : if (dx <= 0) return Flx_copy(x0);
1835 157390 : dy = dx*d;
1836 157390 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1837 157388 : z = y + 2;
1838 8935243 : for (i=0; i<=dy; i++) z[i] = 0;
1839 4353571 : for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1840 157388 : return y;
1841 : }
1842 :
1843 : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1844 : GEN
1845 148009 : Flx_splitting(GEN p, long k)
1846 : {
1847 148009 : long n = degpol(p), v = p[1], m, i, j, l;
1848 : GEN r;
1849 :
1850 148009 : m = n/k;
1851 148009 : r = cgetg(k+1,t_VEC);
1852 681339 : for(i=1; i<=k; i++)
1853 : {
1854 533333 : gel(r,i) = cgetg(m+3, t_VECSMALL);
1855 533329 : mael(r,i,1) = v;
1856 : }
1857 4546790 : for (j=1, i=0, l=2; i<=n; i++)
1858 : {
1859 4398784 : mael(r,j,l) = p[2+i];
1860 4398784 : if (j==k) { j=1; l++; } else j++;
1861 : }
1862 681344 : for(i=1; i<=k; i++)
1863 533340 : gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1864 148004 : return r;
1865 : }
1866 :
1867 : /* ux + vy */
1868 : static GEN
1869 416900 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p, ulong pi)
1870 416900 : { return Flx_add(Flx_mul_pre(u,x, p,pi), Flx_mul_pre(v,y, p,pi), p); }
1871 :
1872 : static GEN
1873 24752 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p, ulong pi)
1874 : {
1875 24752 : GEN res = cgetg(3, t_COL);
1876 24752 : gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p, pi);
1877 24752 : gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p, pi);
1878 24752 : return res;
1879 : }
1880 :
1881 : #if 0
1882 : static GEN
1883 : FlxM_mul2_old(GEN M, GEN N, ulong p)
1884 : {
1885 : GEN res = cgetg(3, t_MAT);
1886 : gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1887 : gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1888 : return res;
1889 : }
1890 : #endif
1891 : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1892 : static GEN
1893 6517 : FlxM_mul2(GEN A, GEN B, ulong p, ulong pi)
1894 : {
1895 6517 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1896 6517 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1897 6517 : GEN M1 = Flx_mul_pre(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p, pi);
1898 6517 : GEN M2 = Flx_mul_pre(Flx_add(A21,A22, p), B11, p, pi);
1899 6517 : GEN M3 = Flx_mul_pre(A11, Flx_sub(B12,B22, p), p, pi);
1900 6517 : GEN M4 = Flx_mul_pre(A22, Flx_sub(B21,B11, p), p, pi);
1901 6517 : GEN M5 = Flx_mul_pre(Flx_add(A11,A12, p), B22, p, pi);
1902 6517 : GEN M6 = Flx_mul_pre(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p, pi);
1903 6517 : GEN M7 = Flx_mul_pre(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p, pi);
1904 6517 : GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1905 6517 : GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1906 6517 : retmkmat22(Flx_add(T1,T2, p), Flx_add(M3,M5, p),
1907 : Flx_add(M2,M4, p), Flx_add(T3,T4, p));
1908 : }
1909 :
1910 : /* Return [0,1;1,-q]*M */
1911 : static GEN
1912 6345 : Flx_FlxM_qmul(GEN q, GEN M, ulong p, ulong pi)
1913 : {
1914 6345 : GEN u = Flx_mul_pre(gcoeff(M,2,1), q, p, pi);
1915 6345 : GEN v = Flx_mul_pre(gcoeff(M,2,2), q, p, pi);
1916 6345 : retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
1917 : Flx_sub(gcoeff(M,1,1), u, p), Flx_sub(gcoeff(M,1,2), v, p));
1918 : }
1919 :
1920 : static GEN
1921 895 : matid2_FlxM(long v)
1922 895 : { retmkmat22(pol1_Flx(v),pol0_Flx(v),pol0_Flx(v),pol1_Flx(v)); }
1923 :
1924 : static GEN
1925 13 : matJ2_FlxM(long v)
1926 13 : { retmkmat22(pol0_Flx(v),pol1_Flx(v),pol1_Flx(v),pol0_Flx(v)); }
1927 :
1928 : struct Flx_res
1929 : {
1930 : ulong res, lc;
1931 : long deg0, deg1, off;
1932 : };
1933 :
1934 : INLINE void
1935 9405 : Flx_halfres_update_pre(long da, long db, long dr, ulong p, ulong pi, struct Flx_res *res)
1936 : {
1937 9405 : if (dr >= 0)
1938 : {
1939 9405 : if (res->lc != 1)
1940 : {
1941 7596 : if (pi)
1942 : {
1943 3127 : res->lc = Fl_powu_pre(res->lc, da - dr, p, pi);
1944 3127 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1945 : } else
1946 : {
1947 4469 : res->lc = Fl_powu(res->lc, da - dr, p);
1948 4469 : res->res = Fl_mul(res->res, res->lc, p);
1949 : }
1950 : }
1951 9405 : if (both_odd(da + res->off, db + res->off))
1952 63 : res->res = Fl_neg(res->res, p);
1953 : } else
1954 : {
1955 0 : if (db == 0)
1956 : {
1957 0 : if (res->lc != 1)
1958 : {
1959 0 : if (pi)
1960 : {
1961 0 : res->lc = Fl_powu_pre(res->lc, da, p, pi);
1962 0 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1963 : } else
1964 : {
1965 0 : res->lc = Fl_powu(res->lc, da, p);
1966 0 : res->res = Fl_mul(res->res, res->lc, p);
1967 : }
1968 : }
1969 : } else
1970 0 : res->res = 0;
1971 : }
1972 9405 : }
1973 :
1974 : static GEN
1975 1109658 : Flx_halfres_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *pa, GEN *pb, struct Flx_res *res)
1976 : {
1977 1109658 : pari_sp av = avma;
1978 : GEN u, u1, v, v1, M;
1979 1109658 : long vx = a[1], n = lgpol(a)>>1;
1980 1109656 : u1 = v = pol0_Flx(vx);
1981 1109650 : u = v1 = pol1_Flx(vx);
1982 6851042 : while (lgpol(b)>n)
1983 : {
1984 : GEN r, q;
1985 5741397 : q = Flx_divrem_pre(a,b,p,pi, &r);
1986 5741415 : if (res)
1987 : {
1988 8362 : long da = degpol(a), db=degpol(b), dr = degpol(r);
1989 8362 : res->lc = b[db+2];
1990 8362 : if (dr >= n)
1991 7133 : Flx_halfres_update_pre(da, db, dr, p, pi, res);
1992 : else
1993 : {
1994 1229 : res->deg0 = da;
1995 1229 : res->deg1 = db;
1996 : }
1997 : }
1998 5741415 : a = b; b = r; swap(u,u1); swap(v,v1);
1999 5741415 : u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
2000 5741358 : v1 = Flx_sub(v1, Flx_mul(v, q, p), p);
2001 5741398 : if (gc_needed(av,2))
2002 : {
2003 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
2004 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
2005 : }
2006 : }
2007 1109504 : M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
2008 1109632 : return gc_all(av,3, &M, pa, pb);
2009 : }
2010 :
2011 : static GEN Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res);
2012 :
2013 : static GEN
2014 19282 : Flx_halfres_split(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2015 : {
2016 19282 : pari_sp av = avma;
2017 : GEN R, S, T, V1, V2;
2018 : GEN x1, y1, r, q;
2019 19282 : long l = lgpol(x), n = l>>1, k;
2020 19282 : if (lgpol(y) <= n)
2021 855 : { *a = Flx_copy(x); *b = Flx_copy(y); return matid2_FlxM(x[1]); }
2022 18427 : if (res)
2023 : {
2024 3263 : res->lc = Flx_lead(y);
2025 3263 : res->deg0 -= n;
2026 3263 : res->deg1 -= n;
2027 3263 : res->off += n;
2028 : }
2029 18427 : R = Flx_halfres_i(Flx_shift(x,-n),Flx_shift(y,-n),p,pi,a,b,res);
2030 18427 : if (res)
2031 : {
2032 3263 : res->off -= n;
2033 3263 : res->deg0 += n;
2034 3263 : res->deg1 += n;
2035 : }
2036 18427 : V1 = FlxM_Flx_mul2(R, Flxn_red(x,n), Flxn_red(y,n), p, pi);
2037 18427 : x1 = Flx_add(Flx_shift(*a,n), gel(V1,1), p);
2038 18427 : y1 = Flx_add(Flx_shift(*b,n), gel(V1,2), p);
2039 18426 : if (lgpol(y1) <= n)
2040 12101 : { *a = x1; *b = y1; return gc_all(av, 3, &R, a, b); }
2041 6325 : k = 2*n-degpol(y1);
2042 6325 : q = Flx_divrem_pre(x1, y1, p, pi, &r);
2043 6325 : if (res)
2044 : {
2045 1043 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
2046 1043 : if (dy1 < degpol(y))
2047 185 : Flx_halfres_update_pre(res->deg0, res->deg1, dy1, p, pi, res);
2048 1043 : res->lc = uel(y1, dy1+2);
2049 1043 : res->deg0 = dx1;
2050 1043 : res->deg1 = dy1;
2051 1043 : if (dr >= n)
2052 : {
2053 1043 : Flx_halfres_update_pre(dx1, dy1, dr, p, pi, res);
2054 1043 : res->deg0 = dy1;
2055 1043 : res->deg1 = dr;
2056 : }
2057 1043 : res->deg0 -= k;
2058 1043 : res->deg1 -= k;
2059 1043 : res->off += k;
2060 : }
2061 6325 : S = Flx_halfres_i(Flx_shift(y1,-k), Flx_shift(r,-k), p, pi, a, b, res);
2062 6325 : if (res)
2063 : {
2064 1043 : res->deg0 += k;
2065 1043 : res->deg1 += k;
2066 1043 : res->off -= k;
2067 : }
2068 6325 : T = FlxM_mul2(S, Flx_FlxM_qmul(q, R, p,pi), p, pi);
2069 6325 : V2 = FlxM_Flx_mul2(S, Flxn_red(y1,k), Flxn_red(r,k), p, pi);
2070 6325 : *a = Flx_add(Flx_shift(*a,k), gel(V2,1), p);
2071 6325 : *b = Flx_add(Flx_shift(*b,k), gel(V2,2), p);
2072 6325 : return gc_all(av, 3, &T, a, b);
2073 : }
2074 :
2075 : static GEN
2076 1128944 : Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2077 : {
2078 1128944 : if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
2079 1109659 : return Flx_halfres_basecase(x, y, p, pi, a, b, res);
2080 19282 : return Flx_halfres_split(x, y, p, pi, a, b, res);
2081 : }
2082 :
2083 : static GEN
2084 1103149 : Flx_halfgcd_all_i(GEN x, GEN y, ulong p, ulong pi, GEN *pa, GEN *pb)
2085 : {
2086 : GEN a, b, R;
2087 1103149 : R = Flx_halfres_i(x, y, p, pi, &a, &b, NULL);
2088 1103151 : if (pa) *pa = a;
2089 1103151 : if (pb) *pb = b;
2090 1103151 : return R;
2091 : }
2092 :
2093 : /* Return M in GL_2(Fl[X]) such that:
2094 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
2095 : */
2096 :
2097 : GEN
2098 1103150 : Flx_halfgcd_all_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b)
2099 : {
2100 : pari_sp av;
2101 : GEN R, q ,r;
2102 1103150 : long lx = lgpol(x), ly = lgpol(y);
2103 1103149 : if (!lx)
2104 : {
2105 0 : if (a) *a = Flx_copy(y);
2106 0 : if (b) *b = Flx_copy(x);
2107 0 : return matJ2_FlxM(x[1]);
2108 : }
2109 1103149 : if (ly < lx) return Flx_halfgcd_all_i(x, y, p, pi, a, b);
2110 8627 : av = avma;
2111 8627 : q = Flx_divrem(y,x,p,&r);
2112 8627 : R = Flx_halfgcd_all_i(x, r, p, pi, a, b);
2113 8627 : gcoeff(R,1,1) = Flx_sub(gcoeff(R,1,1), Flx_mul_pre(q,gcoeff(R,1,2), p,pi), p);
2114 8627 : gcoeff(R,2,1) = Flx_sub(gcoeff(R,2,1), Flx_mul_pre(q,gcoeff(R,2,2), p,pi), p);
2115 8627 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
2116 : }
2117 :
2118 : GEN
2119 154 : Flx_halfgcd_all(GEN x, GEN y, ulong p, GEN *a, GEN *b)
2120 154 : { return Flx_halfgcd_all_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), a, b); }
2121 :
2122 : GEN
2123 846966 : Flx_halfgcd_pre(GEN x, GEN y, ulong p, ulong pi)
2124 846966 : { return Flx_halfgcd_all_pre(x, y, p, pi, NULL, NULL); }
2125 :
2126 : GEN
2127 0 : Flx_halfgcd(GEN x, GEN y, ulong p)
2128 0 : { return Flx_halfgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2129 :
2130 : /*Do not garbage collect*/
2131 : static GEN
2132 82959941 : Flx_gcd_basecase(GEN a, GEN b, ulong p, ulong pi)
2133 : {
2134 82959941 : pari_sp av = avma;
2135 82959941 : ulong iter = 0;
2136 82959941 : if (lg(b) > lg(a)) swap(a, b);
2137 286368933 : while (lgpol(b))
2138 : {
2139 202994851 : GEN c = Flx_rem_pre(a,b,p,pi);
2140 203408992 : iter++; a = b; b = c;
2141 203408992 : if (gc_needed(av,2))
2142 : {
2143 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
2144 0 : gerepileall(av,2, &a,&b);
2145 : }
2146 : }
2147 82908208 : return iter < 2 ? Flx_copy(a) : a;
2148 : }
2149 :
2150 : GEN
2151 84601303 : Flx_gcd_pre(GEN x, GEN y, ulong p, ulong pi)
2152 : {
2153 84601303 : pari_sp av = avma;
2154 : long lim;
2155 84601303 : if (!lgpol(x)) return Flx_copy(y);
2156 82961716 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2157 82967011 : while (lgpol(y) >= lim)
2158 : {
2159 150 : if (lgpol(y)<=(lgpol(x)>>1))
2160 : {
2161 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2162 0 : x = y; y = r;
2163 : }
2164 150 : (void) Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2165 150 : if (gc_needed(av,2))
2166 : {
2167 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
2168 0 : gerepileall(av,2,&x,&y);
2169 : }
2170 : }
2171 82956433 : return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p,pi));
2172 : }
2173 : GEN
2174 32433494 : Flx_gcd(GEN x, GEN y, ulong p)
2175 32433494 : { return Flx_gcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2176 :
2177 : int
2178 8532297 : Flx_is_squarefree(GEN z, ulong p)
2179 : {
2180 8532297 : pari_sp av = avma;
2181 8532297 : GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
2182 8532166 : return gc_bool(av, degpol(d) == 0);
2183 : }
2184 :
2185 : static long
2186 126443 : Flx_is_smooth_squarefree(GEN f, long r, ulong p, ulong pi)
2187 : {
2188 126443 : pari_sp av = avma;
2189 : long i;
2190 126443 : GEN sx = polx_Flx(f[1]), a = sx;
2191 530284 : for(i=1;;i++)
2192 : {
2193 530284 : if (degpol(f)<=r) return gc_long(av,1);
2194 508368 : a = Flxq_powu_pre(Flx_rem_pre(a,f,p,pi), p, f, p, pi);
2195 508160 : if (Flx_equal(a, sx)) return gc_long(av,1);
2196 504447 : if (i==r) return gc_long(av,0);
2197 403494 : f = Flx_div_pre(f, Flx_gcd_pre(Flx_sub(a,sx,p),f,p,pi),p,pi);
2198 : }
2199 : }
2200 :
2201 : static long
2202 8358 : Flx_is_l_pow(GEN x, ulong p)
2203 : {
2204 8358 : ulong i, lx = lgpol(x);
2205 16591 : for (i=1; i<lx; i++)
2206 14898 : if (x[i+2] && i%p) return 0;
2207 1693 : return 1;
2208 : }
2209 :
2210 : int
2211 126427 : Flx_is_smooth_pre(GEN g, long r, ulong p, ulong pi)
2212 : {
2213 : while (1)
2214 8357 : {
2215 126427 : GEN f = Flx_gcd_pre(g, Flx_deriv(g, p), p, pi);
2216 126312 : if (!Flx_is_smooth_squarefree(Flx_div_pre(g, f, p, pi), r, p, pi))
2217 100954 : return 0;
2218 25510 : if (degpol(f)==0) return 1;
2219 8347 : g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2220 : }
2221 : }
2222 : int
2223 74256 : Flx_is_smooth(GEN g, long r, ulong p)
2224 74256 : { return Flx_is_smooth_pre(g, r, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2225 :
2226 : static GEN
2227 6353364 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2228 : {
2229 6353364 : pari_sp av=avma;
2230 : GEN u,v,u1,v1;
2231 6353364 : long vx = a[1];
2232 6353364 : v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2233 6353207 : if (ptu) { u = pol1_Flx(vx); u1 = pol0_Flx(vx); }
2234 28234526 : while (lgpol(b))
2235 : {
2236 21880208 : GEN r, q = Flx_divrem_pre(a,b,p,pi, &r);
2237 21881548 : a = b; b = r;
2238 21881548 : if (ptu)
2239 : {
2240 2435356 : swap(u,u1);
2241 2435356 : u1 = Flx_sub(u1, Flx_mul_pre(u, q, p, pi), p);
2242 : }
2243 21881557 : swap(v,v1);
2244 21881557 : v1 = Flx_sub(v1, Flx_mul_pre(v, q, p, pi), p);
2245 21881294 : if (gc_needed(av,2))
2246 : {
2247 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(a));
2248 0 : gerepileall(av,ptu ? 6: 4, &a,&b,&v,&v1,&u,&u1);
2249 : }
2250 : }
2251 6353233 : if (ptu) *ptu = u;
2252 6353233 : *ptv = v;
2253 6353233 : return a;
2254 : }
2255 :
2256 : static GEN
2257 147553 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2258 : {
2259 : GEN u, v;
2260 147553 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2261 147553 : GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
2262 147553 : long i, n = 0, vs = x[1];
2263 401703 : while (lgpol(y) >= lim)
2264 : {
2265 254150 : if (lgpol(y)<=(lgpol(x)>>1))
2266 : {
2267 26 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2268 26 : x = y; y = r;
2269 26 : gel(V,++n) = mkmat22(pol0_Flx(vs),pol1_Flx(vs),pol1_Flx(vs),Flx_neg(q,p));
2270 : } else
2271 254124 : gel(V,++n) = Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2272 : }
2273 147553 : y = Flx_extgcd_basecase(x,y,p,pi,&u,&v);
2274 254150 : for (i = n; i>1; i--)
2275 : {
2276 106597 : GEN R = gel(V,i);
2277 106597 : GEN u1 = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2278 106597 : GEN v1 = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2279 106597 : u = u1; v = v1;
2280 : }
2281 : {
2282 147553 : GEN R = gel(V,1);
2283 147553 : if (ptu)
2284 6543 : *ptu = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2285 147553 : *ptv = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2286 : }
2287 147553 : return y;
2288 : }
2289 :
2290 : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2291 : * ux + vy = gcd (mod p) */
2292 : GEN
2293 6353361 : Flx_extgcd_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2294 : {
2295 6353361 : pari_sp av = avma;
2296 : GEN d;
2297 6353361 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2298 6353364 : if (lgpol(y) >= lim)
2299 147553 : d = Flx_extgcd_halfgcd(x, y, p, pi, ptu, ptv);
2300 : else
2301 6205798 : d = Flx_extgcd_basecase(x, y, p, pi, ptu, ptv);
2302 6353238 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
2303 : }
2304 : GEN
2305 854714 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2306 854714 : { return Flx_extgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptu, ptv); }
2307 :
2308 : static GEN
2309 1044 : Flx_halfres_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, ulong *r)
2310 : {
2311 : struct Flx_res res;
2312 : GEN R;
2313 : long dB;
2314 :
2315 1044 : res.res = *r;
2316 1044 : res.lc = Flx_lead(y);
2317 1044 : res.deg0 = degpol(x);
2318 1044 : res.deg1 = degpol(y);
2319 1044 : res.off = 0;
2320 1044 : R = Flx_halfres_i(x, y, p, pi, a, b, &res);
2321 1044 : dB = degpol(*b);
2322 1044 : if (dB < degpol(y))
2323 1044 : Flx_halfres_update_pre(res.deg0, res.deg1, dB, p, pi, &res);
2324 1044 : *r = res.res;
2325 1044 : return R;
2326 : }
2327 :
2328 : static ulong
2329 10266174 : Flx_resultant_basecase_pre(GEN a, GEN b, ulong p, ulong pi)
2330 : {
2331 : pari_sp av;
2332 : long da,db,dc;
2333 10266174 : ulong lb, res = 1UL;
2334 : GEN c;
2335 :
2336 10266174 : da = degpol(a);
2337 10266046 : db = degpol(b);
2338 10266195 : if (db > da)
2339 : {
2340 0 : swapspec(a,b, da,db);
2341 0 : if (both_odd(da,db)) res = p-res;
2342 : }
2343 10266195 : else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2344 10266195 : av = avma;
2345 107279914 : while (db)
2346 : {
2347 97037577 : lb = b[db+2];
2348 97037577 : c = Flx_rem_pre(a,b, p,pi);
2349 96708460 : a = b; b = c; dc = degpol(c);
2350 96684350 : if (dc < 0) return gc_long(av,0);
2351 :
2352 96678872 : if (both_odd(da,db)) res = p - res;
2353 96671970 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, da - dc, p, pi), p);
2354 97008856 : if (gc_needed(av,2))
2355 : {
2356 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant (da = %ld)",da);
2357 0 : gerepileall(av,2, &a,&b);
2358 : }
2359 97013719 : da = db; /* = degpol(a) */
2360 97013719 : db = dc; /* = degpol(b) */
2361 : }
2362 10242337 : return gc_ulong(av, Fl_mul(res, Fl_powu_pre(b[2], da, p, pi), p));
2363 : }
2364 :
2365 : ulong
2366 10268029 : Flx_resultant_pre(GEN x, GEN y, ulong p, ulong pi)
2367 : {
2368 10268029 : pari_sp av = avma;
2369 : long lim;
2370 10268029 : ulong res = 1;
2371 10268029 : long dx = degpol(x), dy = degpol(y);
2372 10267705 : if (dx < 0 || dy < 0) return 0;
2373 10266263 : if (dx < dy)
2374 : {
2375 1064948 : swap(x,y);
2376 1064948 : if (both_odd(dx, dy))
2377 1906 : res = Fl_neg(res, p);
2378 : }
2379 10266264 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2380 10267143 : while (lgpol(y) >= lim)
2381 : {
2382 852 : if (lgpol(y)<=(lgpol(x)>>1))
2383 : {
2384 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2385 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2386 0 : ulong ly = y[dy+2];
2387 0 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2388 0 : if (both_odd(dx, dy))
2389 0 : res = Fl_neg(res, p);
2390 0 : x = y; y = r;
2391 : }
2392 852 : (void) Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2393 852 : if (gc_needed(av,2))
2394 : {
2395 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_res (y = %ld)",degpol(y));
2396 0 : gerepileall(av,2,&x,&y);
2397 : }
2398 : }
2399 10266191 : return gc_ulong(av, Fl_mul(res, Flx_resultant_basecase_pre(x, y, p, pi), p));
2400 : }
2401 :
2402 : ulong
2403 4731470 : Flx_resultant(GEN a, GEN b, ulong p)
2404 4731470 : { return Flx_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2405 :
2406 : /* If resultant is 0, *ptU and *ptV are not set */
2407 : static ulong
2408 53 : Flx_extresultant_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2409 : {
2410 53 : GEN z,q,u,v, x = a, y = b;
2411 53 : ulong lb, res = 1UL;
2412 53 : pari_sp av = avma;
2413 : long dx, dy, dz;
2414 53 : long vs = a[1];
2415 :
2416 53 : u = pol0_Flx(vs);
2417 53 : v = pol1_Flx(vs); /* v = 1 */
2418 53 : dx = degpol(x);
2419 53 : dy = degpol(y);
2420 764 : while (dy)
2421 : { /* b u = x (a), b v = y (a) */
2422 711 : lb = y[dy+2];
2423 711 : q = Flx_divrem_pre(x,y, p, pi, &z);
2424 711 : x = y; y = z; /* (x,y) = (y, x - q y) */
2425 711 : dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2426 711 : z = Flx_sub(u, Flx_mul_pre(q,v, p, pi), p);
2427 711 : u = v; v = z; /* (u,v) = (v, u - q v) */
2428 :
2429 711 : if (both_odd(dx,dy)) res = p - res;
2430 711 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, dx-dz, p, pi), p);
2431 711 : dx = dy; /* = degpol(x) */
2432 711 : dy = dz; /* = degpol(y) */
2433 : }
2434 53 : res = Fl_mul(res, Fl_powu_pre(y[2], dx, p, pi), p);
2435 53 : lb = Fl_mul(res, Fl_inv(y[2],p), p);
2436 53 : v = gerepileuptoleaf(av, Flx_Fl_mul_pre(v, lb, p, pi));
2437 53 : av = avma;
2438 53 : u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul_pre(b,v,p,pi), p);
2439 53 : u = gerepileuptoleaf(av, Flx_div_pre(u,a,p,pi)); /* = (res - b v) / a */
2440 53 : *ptU = u;
2441 53 : *ptV = v; return res;
2442 : }
2443 :
2444 : ulong
2445 53 : Flx_extresultant_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2446 : {
2447 53 : pari_sp av=avma;
2448 : GEN u, v, R;
2449 53 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2450 53 : ulong res = 1, res1;
2451 53 : long dx = degpol(x), dy = degpol(y);
2452 53 : if (dy > dx)
2453 : {
2454 13 : swap(x,y); lswap(dx,dy);
2455 13 : if (both_odd(dx,dy)) res = p-res;
2456 13 : R = matJ2_FlxM(x[1]);
2457 40 : } else R = matid2_FlxM(x[1]);
2458 53 : if (dy < 0) return 0;
2459 245 : while (lgpol(y) >= lim)
2460 : {
2461 : GEN M;
2462 192 : if (lgpol(y)<=(lgpol(x)>>1))
2463 : {
2464 20 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2465 20 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2466 20 : ulong ly = y[dy+2];
2467 20 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2468 20 : if (both_odd(dx, dy))
2469 0 : res = Fl_neg(res, p);
2470 20 : x = y; y = r;
2471 20 : R = Flx_FlxM_qmul(q, R, p,pi);
2472 : }
2473 192 : M = Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2474 192 : if (!res) return gc_ulong(av, 0);
2475 192 : R = FlxM_mul2(M, R, p, pi);
2476 192 : gerepileall(av,3,&x,&y,&R);
2477 : }
2478 53 : res1 = Flx_extresultant_basecase(x,y,p,pi,&u,&v);
2479 53 : if (!res1) return gc_ulong(av, 0);
2480 53 : *ptU = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi), res, p, pi);
2481 53 : *ptV = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi), res, p, pi);
2482 53 : gerepileall(av, 2, ptU, ptV);
2483 53 : return Fl_mul(res1,res,p);
2484 : }
2485 :
2486 : ulong
2487 53 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2488 53 : { return Flx_extresultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptU, ptV); }
2489 :
2490 : /* allow pi = 0 (SMALL_ULONG) */
2491 : ulong
2492 43773566 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2493 : {
2494 43773566 : ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2495 :
2496 43773566 : if (lx == 1) return 0;
2497 41016072 : x++;
2498 41016072 : if (pi)
2499 : {
2500 : LOCAL_OVERFLOW;
2501 : LOCAL_HIREMAINDER;
2502 40953000 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2503 97814515 : while (++i < lx)
2504 : {
2505 56861515 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2506 56861515 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2507 : }
2508 81118 : return v1? remlll_pre(v1, h1, l1, p, pi)
2509 41034118 : : remll_pre(h1, l1, p, pi);
2510 : }
2511 : else
2512 : {
2513 63072 : l1 = x[i] * y[i];
2514 30920402 : while (++i < lx) { l1 += x[i] * y[i]; if (l1 & HIGHBIT) l1 %= p; }
2515 63072 : return l1 % p;
2516 : }
2517 : }
2518 :
2519 : /* allow pi = 0 (SMALL_ULONG) */
2520 : ulong
2521 100639847 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2522 : {
2523 100639847 : long i, n = degpol(x);
2524 : ulong t;
2525 100641135 : if (n <= 0) return n? 0: x[2];
2526 32923214 : if (n > 15)
2527 : {
2528 180428 : pari_sp av = avma;
2529 180428 : GEN v = Fl_powers_pre(y, n, p, pi);
2530 180437 : return gc_ulong(av, Flx_eval_powers_pre(x, v, p, pi));
2531 : }
2532 32742786 : i = n+2; t = x[i];
2533 32742786 : if (pi)
2534 : {
2535 123046021 : for (i--; i>=2; i--) t = Fl_addmul_pre(uel(x, i), t, y, p, pi);
2536 31636697 : return t;
2537 : }
2538 2670981 : for (i--; i>=2; i--) t = (t * y + x[i]) % p;
2539 1115237 : return t %= p;
2540 : }
2541 : ulong
2542 20382292 : Flx_eval(GEN x, ulong y, ulong p)
2543 20382292 : { return Flx_eval_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2544 :
2545 : ulong
2546 3255 : Flv_prod_pre(GEN x, ulong p, ulong pi)
2547 : {
2548 3255 : pari_sp ltop = avma;
2549 : GEN v;
2550 3255 : long i,k,lx = lg(x);
2551 3255 : if (lx == 1) return 1UL;
2552 3255 : if (lx == 2) return uel(x,1);
2553 3003 : v = cgetg(1+(lx << 1), t_VECSMALL);
2554 3003 : k = 1;
2555 28593 : for (i=1; i<lx-1; i+=2)
2556 25590 : uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2557 3003 : if (i < lx) uel(v,k++) = uel(x,i);
2558 13529 : while (k > 2)
2559 : {
2560 10526 : lx = k; k = 1;
2561 36116 : for (i=1; i<lx-1; i+=2)
2562 25590 : uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2563 10526 : if (i < lx) uel(v,k++) = uel(v,i);
2564 : }
2565 3003 : return gc_ulong(ltop, uel(v,1));
2566 : }
2567 :
2568 : ulong
2569 0 : Flv_prod(GEN v, ulong p)
2570 : {
2571 0 : return Flv_prod_pre(v, p, get_Fl_red(p));
2572 : }
2573 :
2574 : GEN
2575 0 : FlxV_prod(GEN V, ulong p)
2576 : {
2577 : struct _Flxq D;
2578 0 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2579 0 : return gen_product(V, (void *)&D, &_Flx_mul);
2580 : }
2581 :
2582 : /* compute prod (x - a[i]) */
2583 : GEN
2584 737469 : Flv_roots_to_pol(GEN a, ulong p, long vs)
2585 : {
2586 : struct _Flxq D;
2587 737469 : long i,k,lx = lg(a);
2588 : GEN p1;
2589 737469 : if (lx == 1) return pol1_Flx(vs);
2590 737469 : p1 = cgetg(lx, t_VEC);
2591 11903929 : for (k=1,i=1; i<lx-1; i+=2)
2592 11164839 : gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2593 11166731 : Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2594 737198 : if (i < lx)
2595 58112 : gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2596 737196 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2597 737194 : setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2598 : }
2599 :
2600 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2601 : INLINE void
2602 18911486 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2603 : {
2604 18911486 : pari_sp av = avma;
2605 18911486 : long n = lg(w), i;
2606 : ulong u;
2607 : GEN c;
2608 :
2609 18911486 : if (n == 1) return;
2610 18911486 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2611 80282240 : for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2612 19029518 : i = n-1; u = Fl_inv(c[i], p);
2613 80600843 : for ( ; i > 1; --i)
2614 : {
2615 61498302 : ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2616 61472258 : u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2617 : }
2618 19102541 : v[1] = u; set_avma(av);
2619 : }
2620 :
2621 : void
2622 18304726 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2623 :
2624 : GEN
2625 10850 : Flv_inv_pre(GEN w, ulong p, ulong pi)
2626 10850 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2627 :
2628 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2629 : INLINE void
2630 49734 : Flv_inv_indir(GEN w, GEN v, ulong p)
2631 : {
2632 49734 : pari_sp av = avma;
2633 49734 : long n = lg(w), i;
2634 : ulong u;
2635 : GEN c;
2636 :
2637 49734 : if (n == 1) return;
2638 49734 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2639 1718534 : for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2640 49753 : i = n-1; u = Fl_inv(c[i], p);
2641 1718586 : for ( ; i > 1; --i)
2642 : {
2643 1668852 : ulong t = Fl_mul(u, c[i-1], p);
2644 1668852 : u = Fl_mul(u, w[i], p); v[i] = t;
2645 : }
2646 49734 : v[1] = u; set_avma(av);
2647 : }
2648 : static void
2649 635683 : Flv_inv_i(GEN v, GEN w, ulong p)
2650 : {
2651 635683 : if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2652 585949 : else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2653 635687 : }
2654 : void
2655 12017 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2656 : GEN
2657 623668 : Flv_inv(GEN w, ulong p)
2658 623668 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2659 :
2660 : GEN
2661 33024662 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2662 : {
2663 33024662 : long l = lg(a), i;
2664 : GEN a0, z0, z;
2665 33024662 : if (l <= 3)
2666 : {
2667 0 : if (rem) *rem = l == 2? 0: a[2];
2668 0 : return zero_Flx(a[1]);
2669 : }
2670 33024662 : z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2671 32893761 : a0 = a + l-1;
2672 32893761 : z0 = z + l-2; *z0 = *a0--;
2673 32893761 : if (SMALL_ULONG(p))
2674 : {
2675 79748584 : for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2676 : {
2677 59085238 : ulong t = (*a0-- + x * *z0--) % p;
2678 59085238 : *z0 = (long)t;
2679 : }
2680 20663346 : if (rem) *rem = (*a0 + x * *z0) % p;
2681 : }
2682 : else
2683 : {
2684 48246755 : for (i=l-3; i>1; i--)
2685 : {
2686 35996319 : ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2687 36016340 : *z0 = (long)t;
2688 : }
2689 12250436 : if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2690 : }
2691 32912190 : return z;
2692 : }
2693 :
2694 : /* xa, ya = t_VECSMALL */
2695 : static GEN
2696 624863 : Flv_producttree(GEN xa, GEN s, ulong p, ulong pi, long vs)
2697 : {
2698 624863 : long n = lg(xa)-1;
2699 624863 : long m = n==1 ? 1: expu(n-1)+1;
2700 624863 : long i, j, k, ls = lg(s);
2701 624863 : GEN T = cgetg(m+1, t_VEC);
2702 624862 : GEN t = cgetg(ls, t_VEC);
2703 7832868 : for (j=1, k=1; j<ls; k+=s[j++])
2704 7207865 : gel(t, j) = s[j] == 1 ?
2705 7208007 : mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2706 1516060 : mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2707 1516064 : Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2708 624861 : gel(T,1) = t;
2709 2356054 : for (i=2; i<=m; i++)
2710 : {
2711 1731252 : GEN u = gel(T, i-1);
2712 1731252 : long n = lg(u)-1;
2713 1731252 : GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2714 8313608 : for (j=1, k=1; k<n; j++, k+=2)
2715 6582415 : gel(t, j) = Flx_mul_pre(gel(u, k), gel(u, k+1), p, pi);
2716 1731193 : gel(T, i) = t;
2717 : }
2718 624802 : return T;
2719 : }
2720 :
2721 : static GEN
2722 665173 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p, ulong pi)
2723 : {
2724 : long i,j,k;
2725 665173 : long m = lg(T)-1;
2726 665173 : GEN R = cgetg(lg(xa), t_VECSMALL);
2727 665170 : GEN Tp = cgetg(m+1, t_VEC), t;
2728 665166 : gel(Tp, m) = mkvec(P);
2729 2581738 : for (i=m-1; i>=1; i--)
2730 : {
2731 1916575 : GEN u = gel(T, i), v = gel(Tp, i+1);
2732 1916575 : long n = lg(u)-1;
2733 1916575 : t = cgetg(n+1, t_VEC);
2734 9529416 : for (j=1, k=1; k<n; j++, k+=2)
2735 : {
2736 7612850 : gel(t, k) = Flx_rem_pre(gel(v, j), gel(u, k), p, pi);
2737 7612811 : gel(t, k+1) = Flx_rem_pre(gel(v, j), gel(u, k+1), p, pi);
2738 : }
2739 1916566 : gel(Tp, i) = t;
2740 : }
2741 : {
2742 665163 : GEN u = gel(T, i+1), v = gel(Tp, i+1);
2743 665163 : long n = lg(u)-1;
2744 8945504 : for (j=1, k=1; j<=n; j++)
2745 : {
2746 8280309 : long c, d = degpol(gel(u,j));
2747 18326042 : for (c=1; c<=d; c++, k++) R[k] = Flx_eval_pre(gel(v, j), xa[k], p, pi);
2748 : }
2749 665195 : return gc_const((pari_sp)R, R);
2750 : }
2751 : }
2752 :
2753 : static GEN
2754 1386232 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, ulong pi, long vs)
2755 : {
2756 1386232 : pari_sp av = avma;
2757 1386232 : long m = lg(T)-1;
2758 1386232 : long i, j, k, ls = lg(s);
2759 1386232 : GEN Tp = cgetg(m+1, t_VEC);
2760 1385879 : GEN t = cgetg(ls, t_VEC);
2761 24935826 : for (j=1, k=1; j<ls; k+=s[j++])
2762 23550088 : if (s[j]==2)
2763 : {
2764 6941194 : ulong a = Fl_mul(ya[k], R[k], p);
2765 6940858 : ulong b = Fl_mul(ya[k+1], R[k+1], p);
2766 6946582 : gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2767 6941166 : Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2768 6944561 : gel(t, j) = Flx_renormalize(gel(t, j), 4);
2769 : }
2770 : else
2771 16608894 : gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2772 1385738 : gel(Tp, 1) = t;
2773 6387133 : for (i=2; i<=m; i++)
2774 : {
2775 5001435 : GEN u = gel(T, i-1);
2776 5001435 : GEN t = cgetg(lg(gel(T,i)), t_VEC);
2777 4998702 : GEN v = gel(Tp, i-1);
2778 4998702 : long n = lg(v)-1;
2779 27110326 : for (j=1, k=1; k<n; j++, k+=2)
2780 22091177 : gel(t, j) = Flx_add(Flx_mul_pre(gel(u, k), gel(v, k+1), p, pi),
2781 22108931 : Flx_mul_pre(gel(u, k+1), gel(v, k), p, pi), p);
2782 5001395 : gel(Tp, i) = t;
2783 : }
2784 1385698 : return gerepileuptoleaf(av, gmael(Tp,m,1));
2785 : }
2786 :
2787 : GEN
2788 0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2789 : {
2790 0 : pari_sp av = avma;
2791 0 : GEN s = producttree_scheme(lg(xa)-1);
2792 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2793 0 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2794 0 : return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p, pi));
2795 : }
2796 :
2797 : static GEN
2798 2471 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p, ulong pi)
2799 45248 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p, pi)) }
2800 :
2801 : GEN
2802 2471 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2803 : {
2804 2471 : pari_sp av = avma;
2805 2471 : GEN s = producttree_scheme(lg(xa)-1);
2806 2471 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2807 2471 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2808 2471 : return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p, pi));
2809 : }
2810 :
2811 : GEN
2812 368326 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2813 : {
2814 368326 : pari_sp av = avma;
2815 368326 : GEN s = producttree_scheme(lg(xa)-1);
2816 368327 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2817 368326 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2818 368326 : long m = lg(T)-1;
2819 368326 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2820 368328 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2821 368329 : return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, pi, vs));
2822 : }
2823 :
2824 : GEN
2825 101074 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2826 : {
2827 101074 : pari_sp av = avma;
2828 101074 : GEN s = producttree_scheme(lg(xa)-1);
2829 101072 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2830 101072 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2831 101074 : long i, m = lg(T)-1, l = lg(ya)-1;
2832 101074 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2833 101074 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2834 101074 : GEN M = cgetg(l+1, t_VEC);
2835 1118768 : for (i=1; i<=l; i++)
2836 1017697 : gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, pi, vs);
2837 101071 : return gerepileupto(av, M);
2838 : }
2839 :
2840 : GEN
2841 152995 : Flv_invVandermonde(GEN L, ulong den, ulong p)
2842 : {
2843 152995 : pari_sp av = avma;
2844 152995 : long i, n = lg(L);
2845 : GEN M, R;
2846 152995 : GEN s = producttree_scheme(n-1);
2847 152995 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2848 152995 : GEN tree = Flv_producttree(L, s, p, pi, 0);
2849 152995 : long m = lg(tree)-1;
2850 152995 : GEN T = gmael(tree, m, 1);
2851 152995 : R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p, pi), p);
2852 152995 : if (den!=1) R = Flv_Fl_mul(R, den, p);
2853 152995 : M = cgetg(n, t_MAT);
2854 600537 : for (i = 1; i < n; i++)
2855 : {
2856 447542 : GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2857 447542 : gel(M,i) = Flx_to_Flv(P, n-1);
2858 : }
2859 152995 : return gerepilecopy(av, M);
2860 : }
2861 :
2862 : /***********************************************************************/
2863 : /** Flxq **/
2864 : /***********************************************************************/
2865 : /* Flxq objects are Flx modulo another Flx called q. */
2866 :
2867 : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2868 : GEN
2869 192135925 : Flxq_mul_pre(GEN x,GEN y,GEN T,ulong p,ulong pi)
2870 192135925 : { return Flx_rem_pre(Flx_mul_pre(x,y,p,pi),T,p,pi); }
2871 : GEN
2872 13191915 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2873 13191915 : { return Flxq_mul_pre(x,y,T,p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2874 :
2875 : GEN
2876 278693904 : Flxq_sqr_pre(GEN x,GEN T,ulong p,ulong pi)
2877 278693904 : { return Flx_rem_pre(Flx_sqr_pre(x, p,pi), T, p,pi); }
2878 : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2879 : GEN
2880 2757885 : Flxq_sqr(GEN x,GEN T,ulong p)
2881 2757885 : { return Flxq_sqr_pre(x,T,p,SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2882 :
2883 : static GEN
2884 1551415 : _Flxq_red(void *E, GEN x)
2885 1551415 : { struct _Flxq *s = (struct _Flxq *)E;
2886 1551415 : return Flx_rem_pre(x, s->T, s->p, s->pi); }
2887 : #if 0
2888 : static GEN
2889 : _Flx_sub(void *E, GEN x, GEN y)
2890 : { struct _Flxq *s = (struct _Flxq *)E;
2891 : return Flx_sub(x,y,s->p); }
2892 : #endif
2893 : static GEN
2894 270802241 : _Flxq_sqr(void *data, GEN x)
2895 : {
2896 270802241 : struct _Flxq *D = (struct _Flxq*)data;
2897 270802241 : return Flxq_sqr_pre(x, D->T, D->p, D->pi);
2898 : }
2899 : static GEN
2900 151212694 : _Flxq_mul(void *data, GEN x, GEN y)
2901 : {
2902 151212694 : struct _Flxq *D = (struct _Flxq*)data;
2903 151212694 : return Flxq_mul_pre(x,y, D->T, D->p, D->pi);
2904 : }
2905 : static GEN
2906 22219121 : _Flxq_one(void *data)
2907 : {
2908 22219121 : struct _Flxq *D = (struct _Flxq*)data;
2909 22219121 : return pol1_Flx(get_Flx_var(D->T));
2910 : }
2911 :
2912 : static GEN
2913 22885091 : _Flxq_powu_i(struct _Flxq *D, GEN x, ulong n)
2914 22885091 : { return gen_powu_i(x, n, (void*)D, &_Flxq_sqr, &_Flxq_mul); }
2915 : static GEN
2916 68 : _Flxq_powu(struct _Flxq *D, GEN x, ulong n)
2917 68 : { pari_sp av = avma; return gerepileuptoleaf(av, _Flxq_powu_i(D, x, n)); }
2918 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2919 : GEN
2920 24134162 : Flxq_powu_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
2921 : {
2922 : pari_sp av;
2923 : struct _Flxq D;
2924 24134162 : switch(n)
2925 : {
2926 0 : case 0: return pol1_Flx(get_Flx_var(T));
2927 275251 : case 1: return Flx_copy(x);
2928 972981 : case 2: return Flxq_sqr_pre(x, T, p, pi);
2929 : }
2930 22885930 : av = avma; set_Flxq_pre(&D, T, p, pi);
2931 22884881 : return gerepileuptoleaf(av, _Flxq_powu_i(&D, x, n));
2932 : }
2933 : GEN
2934 488187 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2935 488187 : { return Flxq_powu_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2936 :
2937 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2938 : GEN
2939 25597882 : Flxq_pow_pre(GEN x, GEN n, GEN T, ulong p, ulong pi)
2940 : {
2941 25597882 : pari_sp av = avma;
2942 : struct _Flxq D;
2943 : GEN y;
2944 25597882 : long s = signe(n);
2945 25597882 : if (!s) return pol1_Flx(get_Flx_var(T));
2946 25506881 : if (s < 0) x = Flxq_inv_pre(x,T,p,pi);
2947 25506877 : if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2948 24868279 : set_Flxq_pre(&D, T, p, pi);
2949 24868298 : y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2950 24868228 : return gerepileuptoleaf(av, y);
2951 : }
2952 : GEN
2953 931126 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
2954 931126 : { return Flxq_pow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2955 :
2956 : GEN
2957 28 : Flxq_pow_init_pre(GEN x, GEN n, long k, GEN T, ulong p, ulong pi)
2958 : {
2959 28 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2960 28 : return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2961 : }
2962 : GEN
2963 0 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
2964 0 : { return Flxq_pow_init_pre(x, n, k, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2965 :
2966 : GEN
2967 4392 : Flxq_pow_table_pre(GEN R, GEN n, GEN T, ulong p, ulong pi)
2968 : {
2969 4392 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2970 4392 : return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
2971 : }
2972 : GEN
2973 0 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
2974 0 : { return Flxq_pow_table_pre(R, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2975 :
2976 : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
2977 : * not stack clean. */
2978 : GEN
2979 5498654 : Flxq_invsafe_pre(GEN x, GEN T, ulong p, ulong pi)
2980 : {
2981 5498654 : GEN V, z = Flx_extgcd_pre(get_Flx_mod(T), x, p, pi, NULL, &V);
2982 : ulong iz;
2983 5498768 : if (degpol(z)) return NULL;
2984 5498108 : iz = Fl_inv(uel(z,2), p);
2985 5498112 : return Flx_Fl_mul_pre(V, iz, p, pi);
2986 : }
2987 : GEN
2988 668885 : Flxq_invsafe(GEN x, GEN T, ulong p)
2989 668885 : { return Flxq_invsafe_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2990 :
2991 : GEN
2992 4371808 : Flxq_inv_pre(GEN x, GEN T, ulong p, ulong pi)
2993 : {
2994 4371808 : pari_sp av=avma;
2995 4371808 : GEN U = Flxq_invsafe_pre(x, T, p, pi);
2996 4371791 : if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
2997 4371784 : return gerepileuptoleaf(av, U);
2998 : }
2999 : GEN
3000 335768 : Flxq_inv(GEN x, GEN T, ulong p)
3001 335768 : { return Flxq_inv_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3002 :
3003 : GEN
3004 2416590 : Flxq_div_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
3005 : {
3006 2416590 : pari_sp av = avma;
3007 2416590 : return gerepileuptoleaf(av, Flxq_mul_pre(x,Flxq_inv_pre(y,T,p,pi),T,p,pi));
3008 : }
3009 : GEN
3010 237703 : Flxq_div(GEN x, GEN y, GEN T, ulong p)
3011 237703 : { return Flxq_div_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3012 :
3013 : GEN
3014 22219597 : Flxq_powers_pre(GEN x, long l, GEN T, ulong p, ulong pi)
3015 : {
3016 22219597 : int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
3017 22216867 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3018 22214861 : return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
3019 : }
3020 : GEN
3021 232049 : Flxq_powers(GEN x, long l, GEN T, ulong p)
3022 232049 : { return Flxq_powers_pre(x, l, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3023 :
3024 : GEN
3025 170729 : Flxq_matrix_pow_pre(GEN y, long n, long m, GEN P, ulong l, ulong li)
3026 170729 : { return FlxV_to_Flm(Flxq_powers_pre(y,m-1,P,l,li),n); }
3027 : GEN
3028 399 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
3029 399 : { return Flxq_matrix_pow_pre(y, n, m, P, l, SMALL_ULONG(l)? 0: get_Fl_red(l)); }
3030 :
3031 : GEN
3032 13674799 : Flx_Frobenius_pre(GEN T, ulong p, ulong pi)
3033 13674799 : { return Flxq_powu_pre(polx_Flx(get_Flx_var(T)), p, T, p, pi); }
3034 : GEN
3035 86498 : Flx_Frobenius(GEN T, ulong p)
3036 86498 : { return Flx_Frobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3037 :
3038 : GEN
3039 86614 : Flx_matFrobenius_pre(GEN T, ulong p, ulong pi)
3040 : {
3041 86614 : long n = get_Flx_degree(T);
3042 86614 : return Flxq_matrix_pow_pre(Flx_Frobenius_pre(T, p, pi), n, n, T, p, pi);
3043 : }
3044 : GEN
3045 0 : Flx_matFrobenius(GEN T, ulong p)
3046 0 : { return Flx_matFrobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3047 :
3048 : static GEN
3049 12804790 : Flx_blocks_Flm(GEN P, long n, long m)
3050 : {
3051 12804790 : GEN z = cgetg(m+1,t_MAT);
3052 12804513 : long i,j, k=2, l = lg(P);
3053 36689322 : for(i=1; i<=m; i++)
3054 : {
3055 23888540 : GEN zi = cgetg(n+1,t_VECSMALL);
3056 23884809 : gel(z,i) = zi;
3057 110865362 : for(j=1; j<=n; j++)
3058 86980553 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3059 : }
3060 12800782 : return z;
3061 : }
3062 :
3063 : GEN
3064 516137 : Flx_blocks(GEN P, long n, long m)
3065 : {
3066 516137 : GEN z = cgetg(m+1,t_VEC);
3067 515842 : long i,j, k=2, l = lg(P);
3068 1545415 : for(i=1; i<=m; i++)
3069 : {
3070 1030003 : GEN zi = cgetg(n+2,t_VECSMALL);
3071 1029168 : zi[1] = P[1];
3072 1029168 : gel(z,i) = zi;
3073 6456481 : for(j=2; j<n+2; j++)
3074 5427313 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3075 1029168 : zi = Flx_renormalize(zi, n+2);
3076 : }
3077 515412 : return z;
3078 : }
3079 :
3080 : static GEN
3081 12805897 : FlxV_to_Flm_lg(GEN x, long m, long n)
3082 : {
3083 : long i;
3084 12805897 : GEN y = cgetg(n+1, t_MAT);
3085 60838911 : for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
3086 12802857 : return y;
3087 : }
3088 :
3089 : /* allow pi = 0 (SMALL_ULONG) */
3090 : GEN
3091 13004292 : Flx_FlxqV_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3092 : {
3093 13004292 : pari_sp btop, av = avma;
3094 13004292 : long sv = get_Flx_var(T), m = get_Flx_degree(T);
3095 13004548 : long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
3096 : GEN A, B, C, S, g;
3097 13005361 : if (lQ == 0) return pol0_Flx(sv);
3098 12806613 : if (lQ <= l)
3099 : {
3100 6346973 : n = l;
3101 6346973 : d = 1;
3102 : }
3103 : else
3104 : {
3105 6459640 : n = l-1;
3106 6459640 : d = (lQ+n-1)/n;
3107 : }
3108 12806613 : A = FlxV_to_Flm_lg(x, m, n);
3109 12804676 : B = Flx_blocks_Flm(Q, n, d);
3110 12803848 : C = gerepileupto(av, Flm_mul(A, B, p));
3111 12807113 : g = gel(x, l);
3112 12807113 : if (pi && SMALL_ULONG(p)) pi = 0;
3113 12807113 : T = Flx_get_red_pre(T, p, pi);
3114 12806683 : btop = avma;
3115 12806683 : S = Flv_to_Flx(gel(C, d), sv);
3116 23892867 : for (i = d-1; i>0; i--)
3117 : {
3118 11087778 : S = Flx_add(Flxq_mul_pre(S, g, T, p, pi), Flv_to_Flx(gel(C,i), sv), p);
3119 11087385 : if (gc_needed(btop,1))
3120 0 : S = gerepileuptoleaf(btop, S);
3121 : }
3122 12805089 : return gerepileuptoleaf(av, S);
3123 : }
3124 : GEN
3125 5082 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
3126 5082 : { return Flx_FlxqV_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3127 :
3128 : /* allow pi = 0 (SMALL_ULONG) */
3129 : GEN
3130 2404554 : Flx_Flxq_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3131 : {
3132 2404554 : pari_sp av = avma;
3133 : GEN z, V;
3134 2404554 : long d = degpol(Q), rtd;
3135 2404551 : if (d < 0) return pol0_Flx(get_Flx_var(T));
3136 2404460 : rtd = (long) sqrt((double)d);
3137 2404460 : T = Flx_get_red_pre(T, p, pi);
3138 2404461 : V = Flxq_powers_pre(x, rtd, T, p, pi);
3139 2404536 : z = Flx_FlxqV_eval_pre(Q, V, T, p, pi);
3140 2404487 : return gerepileupto(av, z);
3141 : }
3142 : GEN
3143 789349 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
3144 789349 : { return Flx_Flxq_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3145 :
3146 : /* allow pi = 0 (SMALL_ULONG) */
3147 : GEN
3148 0 : FlxC_FlxqV_eval_pre(GEN x, GEN v, GEN T, ulong p, ulong pi)
3149 0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval_pre(gel(x,i), v, T, p, pi)) }
3150 : GEN
3151 0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
3152 0 : { return FlxC_FlxqV_eval_pre(x, v, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3153 :
3154 : /* allow pi = 0 (SMALL_ULONG) */
3155 : GEN
3156 0 : FlxC_Flxq_eval_pre(GEN x, GEN F, GEN T, ulong p, ulong pi)
3157 : {
3158 0 : long d = brent_kung_optpow(get_Flx_degree(T)-1,lg(x)-1,1);
3159 0 : GEN Fp = Flxq_powers_pre(F, d, T, p, pi);
3160 0 : return FlxC_FlxqV_eval_pre(x, Fp, T, p, pi);
3161 : }
3162 : GEN
3163 0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
3164 0 : { return FlxC_Flxq_eval_pre(x, F, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3165 :
3166 : #if 0
3167 : static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
3168 : _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
3169 : #endif
3170 :
3171 : static GEN
3172 46398 : Flxq_autpow_sqr(void *E, GEN x)
3173 : {
3174 46398 : struct _Flxq *D = (struct _Flxq*)E;
3175 46398 : return Flx_Flxq_eval_pre(x, x, D->T, D->p, D->pi);
3176 : }
3177 : static GEN
3178 20772 : Flxq_autpow_msqr(void *E, GEN x)
3179 : {
3180 20772 : struct _Flxq *D = (struct _Flxq*)E;
3181 20772 : return Flx_FlxqV_eval_pre(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p, D->pi);
3182 : }
3183 :
3184 : GEN
3185 31980 : Flxq_autpow_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3186 : {
3187 31980 : pari_sp av = avma;
3188 : struct _Flxq D;
3189 : long d;
3190 31980 : if (n==0) return Flx_rem_pre(polx_Flx(x[1]), T, p, pi);
3191 31973 : if (n==1) return Flx_rem_pre(x, T, p, pi);
3192 31490 : set_Flxq_pre(&D, T, p, pi);
3193 31490 : d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
3194 31490 : D.aut = Flxq_powers_pre(x, d, T, p, D.pi);
3195 31490 : x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
3196 31490 : return gerepilecopy(av, x);
3197 : }
3198 : GEN
3199 7 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
3200 7 : { return Flxq_autpow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3201 :
3202 : GEN
3203 1679 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
3204 : {
3205 1679 : long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
3206 : ulong i, pi;
3207 1679 : pari_sp av = avma;
3208 1679 : GEN xp, V = cgetg(l+2,t_VEC);
3209 1679 : gel(V,1) = polx_Flx(vT); if (l==0) return V;
3210 1679 : gel(V,2) = gcopy(x); if (l==1) return V;
3211 1679 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3212 1679 : T = Flx_get_red_pre(T, p, pi);
3213 1679 : d = brent_kung_optpow(dT-1, l-1, 1);
3214 1679 : xp = Flxq_powers_pre(x, d, T, p, pi);
3215 7082 : for(i = 3; i < l+2; i++)
3216 5403 : gel(V,i) = Flx_FlxqV_eval_pre(gel(V,i-1), xp, T, p, pi);
3217 1679 : return gerepilecopy(av, V);
3218 : }
3219 :
3220 : static GEN
3221 113384 : Flxq_autsum_mul(void *E, GEN x, GEN y)
3222 : {
3223 113384 : struct _Flxq *D = (struct _Flxq*)E;
3224 113384 : GEN T = D->T;
3225 113384 : ulong p = D->p, pi = D->pi;
3226 113384 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3227 113384 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3228 113384 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3229 113384 : GEN V2 = Flxq_powers_pre(phi2, d, T, p, pi);
3230 113384 : GEN phi3 = Flx_FlxqV_eval_pre(phi1, V2, T, p, pi);
3231 113384 : GEN aphi = Flx_FlxqV_eval_pre(a1, V2, T, p, pi);
3232 113384 : GEN a3 = Flxq_mul_pre(aphi, a2, T, p, pi);
3233 113384 : return mkvec2(phi3, a3);
3234 : }
3235 : static GEN
3236 105664 : Flxq_autsum_sqr(void *E, GEN x)
3237 105664 : { return Flxq_autsum_mul(E, x, x); }
3238 :
3239 : static GEN
3240 99245 : Flxq_autsum_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3241 : {
3242 99245 : pari_sp av = avma;
3243 99245 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3244 99245 : x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
3245 99245 : return gerepilecopy(av, x);
3246 : }
3247 : GEN
3248 0 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
3249 0 : { return Flxq_autsum_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3250 :
3251 : static GEN
3252 763555 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
3253 : {
3254 763555 : struct _Flxq *D = (struct _Flxq*)E;
3255 763555 : GEN T = D->T;
3256 763555 : ulong p = D->p, pi = D->pi;
3257 763555 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3258 763555 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3259 763555 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3260 763568 : GEN V1 = Flxq_powers_pre(phi1, d, T, p, pi);
3261 763520 : GEN phi3 = Flx_FlxqV_eval_pre(phi2, V1, T, p, pi);
3262 763514 : GEN aphi = Flx_FlxqV_eval_pre(a2, V1, T, p, pi);
3263 763528 : GEN a3 = Flx_add(a1, aphi, p);
3264 763546 : return mkvec2(phi3, a3);
3265 : }
3266 :
3267 : static GEN
3268 636104 : Flxq_auttrace_sqr(void *E, GEN x)
3269 636104 : { return Flxq_auttrace_mul(E, x, x); }
3270 :
3271 : GEN
3272 935909 : Flxq_auttrace_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3273 : {
3274 935909 : pari_sp av = avma;
3275 : struct _Flxq D;
3276 935909 : set_Flxq_pre(&D, T, p, pi);
3277 935910 : x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
3278 935887 : return gerepilecopy(av, x);
3279 : }
3280 : GEN
3281 0 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
3282 0 : { return Flxq_auttrace_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3283 :
3284 : static long
3285 422755 : bounded_order(ulong p, GEN b, long k)
3286 : {
3287 422755 : GEN a = modii(utoipos(p), b);
3288 : long i;
3289 870660 : for(i = 1; i < k; i++)
3290 : {
3291 547785 : if (equali1(a)) return i;
3292 447905 : a = modii(muliu(a,p),b);
3293 : }
3294 322875 : return 0;
3295 : }
3296 :
3297 : /* n = (p^d-a)\b
3298 : * b = bb*p^vb
3299 : * p^k = 1 [bb]
3300 : * d = m*k+r+vb
3301 : * u = (p^k-1)/bb;
3302 : * v = (p^(r+vb)-a)/b;
3303 : * w = (p^(m*k)-1)/(p^k-1)
3304 : * n = p^r*w*u+v
3305 : * w*u = p^vb*(p^(m*k)-1)/b
3306 : * n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b */
3307 : static GEN
3308 24672172 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p, ulong pi)
3309 : {
3310 24672172 : pari_sp av=avma;
3311 24672172 : long d = get_Flx_degree(T);
3312 24672171 : GEN an = absi_shallow(n), z, q;
3313 24672180 : if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow_pre(x, n, T, p, pi);
3314 423149 : q = powuu(p, d);
3315 423148 : if (dvdii(q, n))
3316 : {
3317 332 : long vn = logint(an, utoipos(p));
3318 332 : GEN autvn = vn==1 ? aut: Flxq_autpow_pre(aut,vn,T,p,pi);
3319 332 : z = Flx_Flxq_eval_pre(x,autvn,T,p,pi);
3320 : } else
3321 : {
3322 422814 : GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
3323 : GEN bb, u, v, autk;
3324 422816 : long vb = Z_lvalrem(b,p,&bb);
3325 422817 : long m, r, k = is_pm1(bb)? 1: bounded_order(p,bb,d);
3326 422817 : if (!k || d-vb < k) return Flxq_pow_pre(x,n, T,p,pi);
3327 99935 : m = (d-vb)/k; r = (d-vb)%k;
3328 99935 : u = diviiexact(subiu(powuu(p,k),1),bb);
3329 99935 : v = diviiexact(subii(powuu(p,r+vb),a),b);
3330 99935 : autk = k==1 ? aut: Flxq_autpow_pre(aut,k,T,p,pi);
3331 99935 : if (r)
3332 : {
3333 606 : GEN autr = r==1 ? aut: Flxq_autpow_pre(aut,r,T,p,pi);
3334 606 : z = Flx_Flxq_eval_pre(x,autr,T,p,pi);
3335 99329 : } else z = x;
3336 99935 : if (m > 1) z = gel(Flxq_autsum_pre(mkvec2(autk, z), m, T, p, pi), 2);
3337 99935 : if (!is_pm1(u)) z = Flxq_pow_pre(z, u, T, p, pi);
3338 99935 : if (signe(v)) z = Flxq_mul_pre(z, Flxq_pow_pre(x, v, T, p, pi), T, p, pi);
3339 : }
3340 100267 : return gerepileupto(av,signe(n)>0 ? z : Flxq_inv_pre(z,T,p,pi));
3341 : }
3342 :
3343 : static GEN
3344 24664784 : _Flxq_pow(void *data, GEN x, GEN n)
3345 : {
3346 24664784 : struct _Flxq *D = (struct _Flxq*)data;
3347 24664784 : return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p, D->pi);
3348 : }
3349 :
3350 : static GEN
3351 41074 : _Flxq_rand(void *data)
3352 : {
3353 41074 : pari_sp av=avma;
3354 41074 : struct _Flxq *D = (struct _Flxq*)data;
3355 : GEN z;
3356 : do
3357 : {
3358 41470 : set_avma(av);
3359 41470 : z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3360 41470 : } while (lgpol(z)==0);
3361 41074 : return z;
3362 : }
3363 :
3364 : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3365 : static GEN
3366 35479 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3367 : {
3368 35479 : pari_sp av = avma;
3369 : GEN q,n_q,ord,ordp, op;
3370 :
3371 35479 : if (a == 1UL) return gen_0;
3372 : /* p > 2 */
3373 :
3374 35479 : ordp = utoi(p - 1);
3375 35479 : ord = get_arith_Z(o);
3376 35479 : if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3377 35479 : if (a == p - 1) /* -1 */
3378 7739 : return gerepileuptoint(av, shifti(ord,-1));
3379 27740 : ordp = gcdii(ordp, ord);
3380 27740 : op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3381 :
3382 27740 : q = NULL;
3383 27740 : if (T)
3384 : { /* we want < g > = Fp^* */
3385 27740 : if (!equalii(ord,ordp)) {
3386 11906 : q = diviiexact(ord,ordp);
3387 11906 : g = Flxq_pow(g,q,T,p);
3388 : }
3389 : }
3390 27740 : n_q = Fp_log(utoi(a), utoipos(uel(g,2)), op, utoipos(p));
3391 27740 : if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
3392 27740 : if (q) n_q = mulii(q, n_q);
3393 27740 : return gerepileuptoint(av, n_q);
3394 : }
3395 :
3396 : static GEN
3397 548304 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3398 : {
3399 548304 : struct _Flxq *f = (struct _Flxq *)E;
3400 548304 : GEN T = f->T;
3401 548304 : ulong p = f->p;
3402 548304 : long d = get_Flx_degree(T);
3403 548304 : if (Flx_equal1(a)) return gen_0;
3404 388500 : if (Flx_equal(a,g)) return gen_1;
3405 174286 : if (!degpol(a))
3406 35479 : return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3407 138807 : if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3408 138779 : return NULL;
3409 28 : return Flxq_log_index(a, g, ord, T, p);
3410 : }
3411 :
3412 : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3413 :
3414 : const struct bb_group *
3415 281868 : get_Flxq_star(void **E, GEN T, ulong p)
3416 : {
3417 281868 : struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3418 281868 : e->T = T; e->p = p; e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3419 281868 : e->aut = Flx_Frobenius_pre(T, p, e->pi);
3420 281868 : *E = (void*)e; return &Flxq_star;
3421 : }
3422 :
3423 : GEN
3424 97253 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3425 : {
3426 : void *E;
3427 97253 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3428 97253 : return gen_order(a,ord,E,S);
3429 : }
3430 :
3431 : GEN
3432 164497 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3433 : {
3434 : void *E;
3435 164497 : pari_sp av = avma;
3436 164497 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3437 164497 : GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3438 164497 : if (lg(F) > 1 && Flxq_log_use_index(veclast(F), T, p))
3439 24381 : v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3440 164497 : return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
3441 : }
3442 :
3443 : GEN
3444 20125 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3445 : {
3446 20125 : if (!lgpol(a))
3447 : {
3448 7 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3449 0 : if (zeta)
3450 0 : *zeta=pol1_Flx(get_Flx_var(T));
3451 0 : return pol0_Flx(get_Flx_var(T));
3452 : }
3453 : else
3454 : {
3455 : void *E;
3456 20118 : pari_sp av = avma;
3457 20118 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3458 20118 : GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3459 20118 : GEN s = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
3460 20118 : if (!s) return gc_NULL(av);
3461 20076 : return gc_all(av, zeta?2:1, &s, zeta);
3462 : }
3463 : }
3464 :
3465 : static GEN
3466 291114 : Flxq_sumautsum_sqr(void *E, GEN xzd)
3467 : {
3468 291114 : struct _Flxq *D = (struct _Flxq*)E;
3469 291114 : pari_sp av = avma;
3470 : GEN xi, zeta, delta, xi2, zeta2, delta2, temp, xipow;
3471 291114 : GEN T = D->T;
3472 291114 : ulong d, p = D-> p, pi = D->pi;
3473 291114 : xi = gel(xzd, 1); zeta = gel(xzd, 2); delta = gel(xzd, 3);
3474 :
3475 291114 : d = brent_kung_optpow(get_Flx_degree(T)-1,3,1);
3476 291114 : xipow = Flxq_powers_pre(xi, d, T, p, pi);
3477 :
3478 291114 : xi2 = Flx_FlxqV_eval_pre(xi, xipow, T, p, pi);
3479 291114 : zeta2 = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(zeta, xipow, T, p, pi), T, p, pi);
3480 291114 : temp = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(delta, xipow, T, p, pi), T, p, pi);
3481 291114 : delta2 = Flx_add(delta, temp, p);
3482 291114 : return gerepilecopy(av, mkvec3(xi2, zeta2, delta2));
3483 : }
3484 :
3485 : static GEN
3486 40480 : Flxq_sumautsum_msqr(void *E, GEN xzd)
3487 : {
3488 40480 : struct _Flxq *D = (struct _Flxq*)E;
3489 40480 : pari_sp av = avma;
3490 : GEN xii, zetai, deltai, xzd2;
3491 40480 : GEN T = D->T, xi0pow = gel(D->aut, 1), zeta0 = gel(D->aut, 2);
3492 40480 : ulong p = D-> p, pi = D->pi;
3493 40480 : xzd2 = Flxq_sumautsum_sqr(E, xzd);
3494 40480 : xii = Flx_FlxqV_eval_pre(gel(xzd2, 1), xi0pow, T, p, pi);
3495 40480 : zetai = Flxq_mul_pre(zeta0, Flx_FlxqV_eval_pre(gel(xzd2, 2), xi0pow, T, p, pi), T, p, pi);
3496 40480 : deltai = Flx_add(gel(xzd2, 3), zetai, p);
3497 :
3498 40480 : return gerepilecopy(av, mkvec3(xii, zetai, deltai));
3499 : }
3500 :
3501 : /*returns a + a^(1+s) + a^(1+s+2s) + ... + a^(1+s+...+is)
3502 : where ax = [a,s] with s an automorphism */
3503 : static GEN
3504 207592 : Flxq_sumautsum_pre(GEN ax, long i, GEN T, ulong p, ulong pi) {
3505 207592 : pari_sp av = avma;
3506 : GEN a, xi, zeta, vec, res;
3507 : struct _Flxq D;
3508 : ulong d;
3509 207592 : D.T = Flx_get_red(T, p); D.p = p; D.pi = pi;
3510 207592 : a = gel(ax, 1); xi = gel(ax,2);
3511 207592 : d = brent_kung_optpow(get_Flx_degree(T)-1,2*(hammingl(i)-1),1);
3512 207592 : zeta = Flx_Flxq_eval_pre(a, xi, T, p, pi);
3513 207592 : D.aut = mkvec2(Flxq_powers_pre(xi, d, T, p, pi), zeta);
3514 :
3515 207592 : vec = gen_powu_fold(mkvec3(xi, zeta, zeta), i, (void *)&D, Flxq_sumautsum_sqr, Flxq_sumautsum_msqr);
3516 207592 : res = Flxq_mul_pre(a, Flx_add(pol1_Flx(get_Flx_var(T)), gel(vec, 3), p), T, p, pi);
3517 :
3518 207592 : return gerepilecopy(av, res);
3519 : }
3520 :
3521 : GEN
3522 232376 : Flxq_sqrt_pre(GEN z, GEN T, ulong p, ulong pi)
3523 : {
3524 232376 : pari_sp av = avma;
3525 : long d;
3526 232376 : if (p==2)
3527 : {
3528 0 : GEN r = F2xq_sqrt(Flx_to_F2x(z), Flx_to_F2x(get_Flx_mod(T)));
3529 0 : return gerepileupto(av, F2x_to_Flx(r));
3530 : }
3531 232376 : d = get_Flx_degree(T);
3532 232376 : if (d==2)
3533 : {
3534 67676 : GEN P = get_Flx_mod(T), s;
3535 67676 : ulong c = uel(P,2), b = uel(P,3), a = uel(P,4);
3536 67676 : ulong y = degpol(z)<1 ? 0: uel(z,3);
3537 67676 : if (a==1 && b==0)
3538 14890 : {
3539 15670 : ulong x = degpol(z)<1 ? Flx_constant(z): uel(z,2);
3540 15670 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), Fl_neg(c, p), p, pi);
3541 15670 : if (!r) return gc_NULL(av);
3542 14890 : s = mkvecsmall3(P[1], uel(r,1), uel(r,2));
3543 : }
3544 : else
3545 : {
3546 52006 : ulong b2 = Fl_halve(b, p), t = Fl_div(b2, a, p);
3547 52006 : ulong D = Fl_sub(Fl_sqr(b2, p), Fl_mul(a, c, p), p);
3548 52006 : ulong x = degpol(z)<1 ? Flx_constant(z): Fl_sub(uel(z,2), Fl_mul(uel(z,3), t, p), p);
3549 52006 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), D, p, pi);
3550 52006 : if (!r) return gc_NULL(av);
3551 49612 : s = mkvecsmall3(P[1], Fl_add(uel(r,1), Fl_mul(uel(r,2),t,p), p), uel(r,2));
3552 : }
3553 64502 : return gerepileuptoleaf(av, Flx_renormalize(s, 4));
3554 : }
3555 164700 : if (lgpol(z)<=1 && odd(d))
3556 : {
3557 11612 : pari_sp av = avma;
3558 11612 : ulong s = Fl_sqrt(Flx_constant(z), p);
3559 11612 : if (s==~0UL) return gc_NULL(av);
3560 11598 : return gerepilecopy(av, Fl_to_Flx(s, get_Flx_var(T)));
3561 : } else
3562 : {
3563 : GEN c, b, new_z, x, y, w, ax;
3564 : ulong p2, beta;
3565 153088 : long v = get_Flx_var(T);
3566 153088 : if (!lgpol(z)) return pol0_Flx(v);
3567 152423 : T = Flx_get_red_pre(T, p, pi);
3568 152423 : ax = mkvec2(NULL, Flx_Frobenius_pre(T, p, pi));
3569 152423 : p2 = p >> 1; /* (p-1) / 2 */
3570 : do {
3571 208264 : do c = random_Flx(d, v, p); while (!lgpol(c));
3572 :
3573 207592 : new_z = Flxq_mul_pre(z, Flxq_sqr_pre(c, T, p, pi), T, p, pi);
3574 207592 : gel(ax, 1) = Flxq_powu_pre(new_z, p2, T, p, pi);
3575 207592 : y = Flxq_sumautsum_pre(ax, d-2, T, p, pi); /* d > 2 */
3576 207592 : b = Flx_Fl_add(y, 1UL, p);
3577 207592 : } while (!lgpol(b));
3578 :
3579 152423 : x = Flxq_mul_pre(new_z, Flxq_sqr_pre(b, T, p, pi), T, p, pi);
3580 152423 : if (degpol(x) > 0) return gc_NULL(av);
3581 145381 : beta = Fl_sqrt_pre(Flx_constant(x), p, pi);
3582 145381 : if (beta==~0UL) return gc_NULL(av);
3583 145381 : w = Flx_Fl_mul(Flxq_inv_pre(Flxq_mul_pre(b, c, T,p,pi), T,p,pi), beta, p);
3584 145381 : return gerepilecopy(av, w);
3585 : }
3586 : }
3587 :
3588 : GEN
3589 232376 : Flxq_sqrt(GEN a, GEN T, ulong p)
3590 232376 : { return Flxq_sqrt_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3591 :
3592 : /* assume T irreducible mod p */
3593 : int
3594 404318 : Flxq_issquare(GEN x, GEN T, ulong p)
3595 : {
3596 404318 : if (lgpol(x) == 0 || p == 2) return 1;
3597 398011 : return krouu(Flxq_norm(x,T,p), p) == 1;
3598 : }
3599 :
3600 : /* assume T irreducible mod p */
3601 : int
3602 0 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3603 : {
3604 : pari_sp av;
3605 : GEN m;
3606 0 : if (n==1) return Flxq_issquare(x, T, p);
3607 0 : if (lgpol(x) == 0 || p == 2) return 1;
3608 0 : av = avma;
3609 0 : m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3610 0 : return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3611 : }
3612 :
3613 : GEN
3614 113589 : Flxq_lroot_fast_pre(GEN a, GEN sqx, GEN T, long p, ulong pi)
3615 : {
3616 113589 : pari_sp av=avma;
3617 113589 : GEN A = Flx_splitting(a,p);
3618 113589 : return gerepileuptoleaf(av, FlxqV_dotproduct_pre(A,sqx,T,p,pi));
3619 : }
3620 : GEN
3621 0 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
3622 0 : { return Flxq_lroot_fast_pre(a, sqx, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3623 :
3624 : GEN
3625 25053 : Flxq_lroot_pre(GEN a, GEN T, long p, ulong pi)
3626 : {
3627 25053 : pari_sp av=avma;
3628 25053 : long n = get_Flx_degree(T), d = degpol(a);
3629 : GEN sqx, V;
3630 25053 : if (n==1) return leafcopy(a);
3631 25053 : if (n==2) return Flxq_powu_pre(a, p, T, p, pi);
3632 25053 : sqx = Flxq_autpow_pre(Flx_Frobenius_pre(T, p, pi), n-1, T, p, pi);
3633 25053 : if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
3634 0 : if (d>=p)
3635 : {
3636 0 : V = Flxq_powers_pre(sqx,p-1,T,p,pi);
3637 0 : return gerepileuptoleaf(av, Flxq_lroot_fast_pre(a,V,T,p,pi));
3638 : } else
3639 0 : return gerepileuptoleaf(av, Flx_Flxq_eval_pre(a,sqx,T,p,pi));
3640 : }
3641 : GEN
3642 0 : Flxq_lroot(GEN a, GEN T, long p)
3643 0 : { return Flxq_lroot_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3644 :
3645 : ulong
3646 443551 : Flxq_norm(GEN x, GEN TB, ulong p)
3647 : {
3648 443551 : GEN T = get_Flx_mod(TB);
3649 443551 : ulong y = Flx_resultant(T, x, p), L = Flx_lead(T);
3650 443551 : if (L==1 || lgpol(x)==0) return y;
3651 0 : return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3652 : }
3653 :
3654 : ulong
3655 4696 : Flxq_trace(GEN x, GEN TB, ulong p)
3656 : {
3657 4696 : pari_sp av = avma;
3658 : ulong t;
3659 4696 : GEN T = get_Flx_mod(TB);
3660 4696 : long n = degpol(T)-1;
3661 4696 : GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3662 4696 : t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3663 4696 : return gc_ulong(av, t);
3664 : }
3665 :
3666 : /*x must be reduced*/
3667 : GEN
3668 3632 : Flxq_charpoly(GEN x, GEN TB, ulong p)
3669 : {
3670 3632 : pari_sp ltop=avma;
3671 3632 : GEN T = get_Flx_mod(TB);
3672 3632 : long vs = evalvarn(fetch_var());
3673 3632 : GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3674 3632 : GEN r = Flx_FlxY_resultant(T, xm1, p);
3675 3632 : r[1] = x[1];
3676 3632 : (void)delete_var(); return gerepileupto(ltop, r);
3677 : }
3678 :
3679 : /* Computing minimal polynomial : */
3680 : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3681 : /* in Algebraic Extensions of Finite Fields' */
3682 :
3683 : /* Let v a linear form, return the linear form z->v(tau*z)
3684 : that is, v*(M_tau) */
3685 :
3686 : static GEN
3687 1693919 : Flxq_transmul_init(GEN tau, GEN T, ulong p, ulong pi)
3688 : {
3689 : GEN bht;
3690 1693919 : GEN h, Tp = get_Flx_red(T, &h);
3691 1693914 : long n = degpol(Tp), vT = Tp[1];
3692 1693902 : GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3693 1693898 : GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3694 1693901 : ft[1] = vT; bt[1] = vT;
3695 1693901 : if (h)
3696 2688 : bht = Flxn_mul_pre(bt, h, n-1, p, pi);
3697 : else
3698 : {
3699 1691213 : GEN bh = Flx_div_pre(Flx_shift(tau, n-1), T, p, pi);
3700 1691208 : bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3701 1691213 : bht[1] = vT;
3702 : }
3703 1693901 : return mkvec3(bt, bht, ft);
3704 : }
3705 :
3706 : static GEN
3707 4087831 : Flxq_transmul(GEN tau, GEN a, long n, ulong p, ulong pi)
3708 : {
3709 4087831 : pari_sp ltop = avma;
3710 : GEN t1, t2, t3, vec;
3711 4087831 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3712 4087831 : if (lgpol(a)==0) return pol0_Flx(a[1]);
3713 4057525 : t2 = Flx_shift(Flx_mul_pre(bt, a, p, pi),1-n);
3714 4057172 : if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
3715 3062893 : t1 = Flx_shift(Flx_mul_pre(ft, a, p, pi),-n);
3716 3062891 : t3 = Flxn_mul_pre(t1, bht, n-1, p, pi);
3717 3062894 : vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3718 3062970 : return gerepileuptoleaf(ltop, vec);
3719 : }
3720 :
3721 : GEN
3722 784733 : Flxq_minpoly_pre(GEN x, GEN T, ulong p, ulong pi)
3723 : {
3724 784733 : pari_sp ltop = avma;
3725 784733 : long vT = get_Flx_var(T), n = get_Flx_degree(T);
3726 : GEN v_x;
3727 784726 : GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3728 784700 : T = Flx_get_red_pre(T, p, pi);
3729 784701 : v_x = Flxq_powers_pre(x, usqrt(2*n), T, p, pi);
3730 1631646 : while (lgpol(tau) != 0)
3731 : {
3732 : long i, j, m, k1;
3733 : GEN M, v, tr, g_prime, c;
3734 846929 : if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3735 846925 : v = random_Flx(n, vT, p);
3736 846961 : tr = Flxq_transmul_init(tau, T, p, pi);
3737 846941 : v = Flxq_transmul(tr, v, n, p, pi);
3738 846962 : m = 2*(n-degpol(g));
3739 846961 : k1 = usqrt(m);
3740 846962 : tr = Flxq_transmul_init(gel(v_x,k1+1), T, p, pi);
3741 846949 : c = cgetg(m+2,t_VECSMALL);
3742 846917 : c[1] = vT;
3743 4087659 : for (i=0; i<m; i+=k1)
3744 : {
3745 3240700 : long mj = minss(m-i, k1);
3746 12662452 : for (j=0; j<mj; j++)
3747 9421445 : uel(c,m+1-(i+j)) = Flx_dotproduct_pre(v, gel(v_x,j+1), p, pi);
3748 3241007 : v = Flxq_transmul(tr, v, n, p, pi);
3749 : }
3750 846959 : c = Flx_renormalize(c, m+2);
3751 : /* now c contains <v,x^i> , i = 0..m-1 */
3752 846961 : M = Flx_halfgcd_pre(monomial_Flx(1, m, vT), c, p, pi);
3753 846967 : g_prime = gmael(M, 2, 2);
3754 846967 : if (degpol(g_prime) < 1) continue;
3755 835093 : g = Flx_mul_pre(g, g_prime, p, pi);
3756 835080 : tau = Flxq_mul_pre(tau, Flx_FlxqV_eval_pre(g_prime, v_x, T,p,pi), T,p,pi);
3757 : }
3758 784677 : g = Flx_normalize(g,p);
3759 784726 : return gerepileuptoleaf(ltop,g);
3760 : }
3761 : GEN
3762 44454 : Flxq_minpoly(GEN x, GEN T, ulong p)
3763 44454 : { return Flxq_minpoly_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3764 :
3765 : GEN
3766 20 : Flxq_conjvec(GEN x, GEN T, ulong p)
3767 : {
3768 20 : long i, l = 1+get_Flx_degree(T);
3769 20 : GEN z = cgetg(l,t_COL);
3770 20 : struct _Flxq D; set_Flxq(&D, T, p);
3771 20 : gel(z,1) = Flx_copy(x);
3772 88 : for (i=2; i<l; i++) gel(z,i) = _Flxq_powu(&D, gel(z,i-1), p);
3773 20 : return z;
3774 : }
3775 :
3776 : GEN
3777 7201 : gener_Flxq(GEN T, ulong p, GEN *po)
3778 : {
3779 7201 : long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
3780 : ulong p_1, pi;
3781 : GEN g, L, L2, o, q, F;
3782 : pari_sp av0, av;
3783 :
3784 7201 : if (f == 1) {
3785 : GEN fa;
3786 28 : o = utoipos(p-1);
3787 28 : fa = Z_factor(o);
3788 28 : L = gel(fa,1);
3789 28 : L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
3790 28 : g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
3791 28 : if (po) *po = mkvec2(o, fa);
3792 28 : return g;
3793 : }
3794 :
3795 7173 : av0 = avma; p_1 = p - 1;
3796 7173 : q = diviuexact(subiu(powuu(p,f), 1), p_1);
3797 :
3798 7173 : L = cgetg(1, t_VECSMALL);
3799 7173 : if (p > 3)
3800 : {
3801 2371 : ulong t = p_1 >> vals(p_1);
3802 2371 : GEN P = gel(factoru(t), 1);
3803 2371 : L = cgetg_copy(P, &i);
3804 3787 : while (--i) L[i] = p_1 / P[i];
3805 : }
3806 7173 : o = factor_pn_1(utoipos(p),f);
3807 7173 : L2 = leafcopy( gel(o, 1) );
3808 19212 : for (i = j = 1; i < lg(L2); i++)
3809 : {
3810 12039 : if (umodui(p_1, gel(L2,i)) == 0) continue;
3811 6488 : gel(L2,j++) = diviiexact(q, gel(L2,i));
3812 : }
3813 7173 : setlg(L2, j); pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3814 7173 : F = Flx_Frobenius_pre(T, p, pi);
3815 17647 : for (av = avma;; set_avma(av))
3816 10474 : {
3817 : GEN tt;
3818 17647 : g = random_Flx(f, vT, p);
3819 17647 : if (degpol(g) < 1) continue;
3820 12082 : if (p == 2) tt = g;
3821 : else
3822 : {
3823 8883 : ulong t = Flxq_norm(g, T, p);
3824 8883 : if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
3825 4764 : tt = Flxq_powu_pre(g, p_1>>1, T, p, pi);
3826 : }
3827 14556 : for (i = 1; i < j; i++)
3828 : {
3829 7383 : GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p, pi);
3830 7383 : if (!degpol(a) && uel(a,2) == p_1) break;
3831 : }
3832 7963 : if (i == j) break;
3833 : }
3834 7173 : if (!po)
3835 : {
3836 187 : set_avma((pari_sp)g);
3837 187 : g = gerepileuptoleaf(av0, g);
3838 : }
3839 : else {
3840 6986 : *po = mkvec2(subiu(powuu(p,f), 1), o);
3841 6986 : gerepileall(av0, 2, &g, po);
3842 : }
3843 7173 : return g;
3844 : }
3845 :
3846 : static GEN
3847 366530 : _Flxq_neg(void *E, GEN x)
3848 366530 : { struct _Flxq *s = (struct _Flxq *)E;
3849 366530 : return Flx_neg(x,s->p); }
3850 :
3851 : static GEN
3852 1462529 : _Flxq_rmul(void *E, GEN x, GEN y)
3853 1462529 : { struct _Flxq *s = (struct _Flxq *)E;
3854 1462529 : return Flx_mul_pre(x,y,s->p,s->pi); }
3855 :
3856 : static GEN
3857 9418 : _Flxq_inv(void *E, GEN x)
3858 9418 : { struct _Flxq *s = (struct _Flxq *)E;
3859 9418 : return Flxq_inv(x,s->T,s->p); }
3860 :
3861 : static int
3862 68965 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
3863 :
3864 : static GEN
3865 6453 : _Flxq_s(void *E, long x)
3866 6453 : { struct _Flxq *s = (struct _Flxq *)E;
3867 6453 : ulong u = x<0 ? s->p+x: (ulong)x;
3868 6453 : return Fl_to_Flx(u, get_Flx_var(s->T));
3869 : }
3870 :
3871 : static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
3872 : _Flxq_inv,_Flxq_equal0,_Flxq_s};
3873 :
3874 68945 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
3875 : {
3876 68945 : GEN z = new_chunk(sizeof(struct _Flxq));
3877 68945 : set_Flxq((struct _Flxq *)z, T, p); *E = (void*)z; return &Flxq_field;
3878 : }
3879 :
3880 : /***********************************************************************/
3881 : /** Flxn **/
3882 : /***********************************************************************/
3883 :
3884 : GEN
3885 54262 : Flx_invLaplace(GEN x, ulong p)
3886 : {
3887 54262 : long i, d = degpol(x);
3888 : ulong t;
3889 : GEN y;
3890 54258 : if (d <= 1) return Flx_copy(x);
3891 54258 : t = Fl_inv(factorial_Fl(d, p), p);
3892 54297 : y = cgetg(d+3, t_VECSMALL);
3893 54251 : y[1] = x[1];
3894 1325558 : for (i=d; i>=2; i--)
3895 : {
3896 1271273 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3897 1271253 : t = Fl_mul(t, i, p);
3898 : }
3899 54285 : uel(y,3) = uel(x,3);
3900 54285 : uel(y,2) = uel(x,2);
3901 54285 : return y;
3902 : }
3903 :
3904 : GEN
3905 27285 : Flx_Laplace(GEN x, ulong p)
3906 : {
3907 27285 : long i, d = degpol(x);
3908 27285 : ulong t = 1;
3909 : GEN y;
3910 27285 : if (d <= 1) return Flx_copy(x);
3911 27285 : y = cgetg(d+3, t_VECSMALL);
3912 27270 : y[1] = x[1];
3913 27270 : uel(y,2) = uel(x,2);
3914 27270 : uel(y,3) = uel(x,3);
3915 755299 : for (i=2; i<=d; i++)
3916 : {
3917 728006 : t = Fl_mul(t, i%p, p);
3918 728010 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3919 : }
3920 27293 : return y;
3921 : }
3922 :
3923 : GEN
3924 6231179 : Flxn_red(GEN a, long n)
3925 : {
3926 6231179 : long i, L, l = lg(a);
3927 : GEN b;
3928 6231179 : if (l == 2 || !n) return zero_Flx(a[1]);
3929 5841581 : L = n+2; if (L > l) L = l;
3930 5841581 : b = cgetg(L, t_VECSMALL); b[1] = a[1];
3931 58607137 : for (i=2; i<L; i++) b[i] = a[i];
3932 5839078 : return Flx_renormalize(b,L);
3933 : }
3934 :
3935 : GEN
3936 5064275 : Flxn_mul_pre(GEN a, GEN b, long n, ulong p, ulong pi)
3937 5064275 : { return Flxn_red(Flx_mul_pre(a, b, p, pi), n); }
3938 : GEN
3939 75308 : Flxn_mul(GEN a, GEN b, long n, ulong p)
3940 75308 : { return Flxn_mul_pre(a, b, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3941 :
3942 : GEN
3943 0 : Flxn_sqr_pre(GEN a, long n, ulong p, ulong pi)
3944 0 : { return Flxn_red(Flx_sqr_pre(a, p, pi), n); }
3945 : GEN
3946 0 : Flxn_sqr(GEN a, long n, ulong p)
3947 0 : { return Flxn_sqr_pre(a, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3948 :
3949 : /* (f*g) \/ x^n */
3950 : static GEN
3951 937507 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p, ulong pi)
3952 937507 : { return Flx_shift(Flx_mul_pre(f, g, p, pi),-n); }
3953 :
3954 : static GEN
3955 516040 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p, ulong pi)
3956 : {
3957 516040 : GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
3958 515599 : return Flx_add(Flx_mulhigh_i(fl, g, n2, p, pi),
3959 : Flxn_mul_pre(fh, g, n - n2, p, pi), p);
3960 : }
3961 :
3962 : /* g==NULL -> assume g==1 */
3963 : GEN
3964 55061 : Flxn_div_pre(GEN g, GEN f, long e, ulong p, ulong pi)
3965 : {
3966 55061 : pari_sp av = avma, av2;
3967 : ulong mask;
3968 : GEN W;
3969 55061 : long n = 1;
3970 55061 : if (lg(f) <= 2) pari_err_INV("Flxn_inv",f);
3971 55061 : W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
3972 55076 : mask = quadratic_prec_mask(e);
3973 55074 : av2 = avma;
3974 258334 : for (;mask>1;)
3975 : {
3976 : GEN u, fr;
3977 203245 : long n2 = n;
3978 203245 : n<<=1; if (mask & 1) n--;
3979 203245 : mask >>= 1;
3980 203245 : fr = Flxn_red(f, n);
3981 203043 : if (mask>1 || !g)
3982 : {
3983 149043 : u = Flxn_mul_pre(W, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
3984 149505 : W = Flx_sub(W, Flx_shift(u, n2), p);
3985 : } else
3986 : {
3987 54000 : GEN y = Flxn_mul_pre(g, W, n, p, pi), yt = Flxn_red(y, n-n2);
3988 53979 : u = Flxn_mul_pre(yt, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
3989 53998 : W = Flx_sub(y, Flx_shift(u, n2), p);
3990 : }
3991 203259 : if (gc_needed(av2,2))
3992 : {
3993 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_div, e = %ld", n);
3994 0 : W = gerepileupto(av2, W);
3995 : }
3996 : }
3997 55089 : return gerepileupto(av, W);
3998 : }
3999 : GEN
4000 55037 : Flxn_div(GEN g, GEN f, long e, ulong p)
4001 55037 : { return Flxn_div_pre(g, f, e, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4002 :
4003 : GEN
4004 1030 : Flxn_inv(GEN f, long e, ulong p)
4005 1030 : { return Flxn_div(NULL, f, e, p); }
4006 :
4007 : GEN
4008 109346 : Flxn_expint(GEN h, long e, ulong p)
4009 : {
4010 109346 : pari_sp av = avma, av2;
4011 109346 : long v = h[1], n=1;
4012 109346 : GEN f = pol1_Flx(v), g = pol1_Flx(v);
4013 109325 : ulong mask = quadratic_prec_mask(e), pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4014 109330 : av2 = avma;
4015 422597 : for (;mask>1;)
4016 : {
4017 : GEN u, w;
4018 422519 : long n2 = n;
4019 422519 : n<<=1; if (mask & 1) n--;
4020 422519 : mask >>= 1;
4021 422519 : u = Flxn_mul_pre(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p,pi), n-n2, p,pi);
4022 422489 : u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
4023 422513 : w = Flxn_mul_pre(f, Flx_integXn(u, n2-1, p), n-n2, p, pi);
4024 422472 : f = Flx_add(f, Flx_shift(w, n2), p);
4025 422610 : if (mask<=1) break;
4026 313279 : u = Flxn_mul_pre(g, Flxn_mulhigh(f, g, n2, n, p, pi), n-n2, p, pi);
4027 313260 : g = Flx_sub(g, Flx_shift(u, n2), p);
4028 313267 : if (gc_needed(av2,2))
4029 : {
4030 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
4031 0 : gerepileall(av2, 2, &f, &g);
4032 : }
4033 : }
4034 109409 : return gerepileupto(av, f);
4035 : }
4036 :
4037 : GEN
4038 0 : Flxn_exp(GEN h, long e, ulong p)
4039 : {
4040 0 : if (degpol(h)<1 || uel(h,2)!=0)
4041 0 : pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
4042 0 : return Flxn_expint(Flx_deriv(h, p), e, p);
4043 : }
4044 :
4045 : INLINE GEN
4046 217044 : Flxn_recip(GEN x, long n)
4047 : {
4048 217044 : GEN z=Flx_recipspec(x+2,lgpol(x),n);
4049 216869 : z[1]=x[1];
4050 216869 : return z;
4051 : }
4052 :
4053 : GEN
4054 53999 : Flx_Newton(GEN P, long n, ulong p)
4055 : {
4056 53999 : pari_sp av = avma;
4057 53999 : long d = degpol(P);
4058 53991 : GEN dP = Flxn_recip(Flx_deriv(P, p), d);
4059 53904 : GEN Q = Flxn_div(dP, Flxn_recip(P, d+1), n, p);
4060 53960 : return gerepileuptoleaf(av, Q);
4061 : }
4062 :
4063 : GEN
4064 109342 : Flx_fromNewton(GEN P, ulong p)
4065 : {
4066 109342 : pari_sp av = avma;
4067 109342 : ulong n = Flx_constant(P)+1;
4068 109341 : GEN z = Flx_neg(Flx_shift(P, -1), p);
4069 109346 : GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
4070 109327 : return gerepileuptoleaf(av, Q);
4071 : }
4072 :
4073 : static void
4074 12514 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
4075 : {
4076 : long i;
4077 : ulong e;
4078 12514 : GEN P = cgetg(d+1, t_VECSMALL);
4079 12514 : GEN V = cgetg(d+1, t_VECSMALL);
4080 1396581 : for (i=1, e=1; i<=d; i++, e++)
4081 : {
4082 1384067 : if (e==p)
4083 : {
4084 459153 : e = 0;
4085 459153 : V[i] = u_lvalrem(i, p, &uel(P,i));
4086 : } else
4087 : {
4088 924914 : V[i] = 0; uel(P,i) = i;
4089 : }
4090 : }
4091 12514 : *pt_P = P; *pt_V = V;
4092 12514 : }
4093 :
4094 : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
4095 : * val large enough to compensate for the power of p in the factorials */
4096 :
4097 : static GEN
4098 497 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long sv)
4099 : {
4100 497 : pari_sp av = avma;
4101 497 : long i, d = n-1, w;
4102 : GEN y, W, E, t;
4103 497 : init_invlaplace(d, p, &E, &W);
4104 497 : t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
4105 497 : w = zv_sum(W);
4106 497 : if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
4107 497 : y = cgetg(d+3,t_POL);
4108 497 : y[1] = evalsigne(1) | sv;
4109 28882 : for (i=d; i>=1; i--)
4110 : {
4111 28385 : gel(y,i+2) = t;
4112 28385 : t = Fp_mulu(t, uel(E,i), q);
4113 28385 : if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
4114 : }
4115 497 : gel(y,2) = t;
4116 497 : return gerepilecopy(av, ZX_renormalize(y, d+3));
4117 : }
4118 :
4119 : GEN
4120 27493 : Flx_composedsum(GEN P, GEN Q, ulong p)
4121 : {
4122 27493 : pari_sp av = avma;
4123 27493 : long n = 1 + degpol(P)*degpol(Q);
4124 27489 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4125 27490 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4126 : GEN R;
4127 27494 : if (p >= (ulong)n)
4128 : {
4129 26997 : GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
4130 27006 : GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
4131 26994 : GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
4132 26998 : R = Flx_fromNewton(L, p);
4133 : } else
4134 : {
4135 497 : long v = factorial_lval(n-1, p);
4136 497 : long w = 1 + ulogint(n-1, p);
4137 497 : GEN pv = powuu(p, v);
4138 497 : GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
4139 497 : GEN iL = ZpX_invLaplace_init(n, q, p, v, P[1]);
4140 497 : GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
4141 497 : GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
4142 497 : GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
4143 497 : GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
4144 497 : R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
4145 : }
4146 27485 : return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
4147 : }
4148 :
4149 : static GEN
4150 3826 : _Flx_composedsum(void *E, GEN a, GEN b)
4151 3826 : { return Flx_composedsum(a, b, (ulong)E); }
4152 :
4153 : GEN
4154 28917 : FlxV_composedsum(GEN V, ulong p)
4155 28917 : { return gen_product(V, (void *)p, &_Flx_composedsum); }
4156 :
4157 : GEN
4158 0 : Flx_composedprod(GEN P, GEN Q, ulong p)
4159 : {
4160 0 : pari_sp av = avma;
4161 0 : long n = 1+ degpol(P)*degpol(Q);
4162 0 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4163 0 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4164 : GEN R;
4165 0 : if (p >= (ulong)n)
4166 : {
4167 0 : GEN L = Flx_convol(Flx_Newton(P,n,p), Flx_Newton(Q,n,p), p);
4168 0 : R = Flx_fromNewton(L, p);
4169 : } else
4170 : {
4171 0 : long w = 1 + ulogint(n, p);
4172 0 : GEN qf = powuu(p, w);
4173 0 : GEN Pl = FpX_convol(FpX_Newton(Flx_to_ZX(P), n, qf), FpX_Newton(Flx_to_ZX(Q), n, qf), qf);
4174 0 : R = ZX_to_Flx(FpX_fromNewton(Pl, qf), p);
4175 : }
4176 0 : return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
4177 :
4178 : }
4179 :
4180 : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
4181 : static GEN
4182 0 : Fl_Xp1_powu(ulong n, ulong p, long v)
4183 : {
4184 0 : ulong k, d = (n + 1) >> 1;
4185 0 : GEN C, V = identity_zv(d);
4186 :
4187 0 : Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
4188 0 : C = cgetg(n+3, t_VECSMALL);
4189 0 : C[1] = v;
4190 0 : uel(C,2) = 1UL;
4191 0 : uel(C,3) = n%p;
4192 0 : uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
4193 : /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
4194 0 : if (SMALL_ULONG(p))
4195 0 : for (k = 3; k <= d; k++)
4196 0 : uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
4197 : else
4198 : {
4199 0 : ulong pi = get_Fl_red(p);
4200 0 : for (k = 3; k <= d; k++)
4201 0 : uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
4202 : }
4203 0 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4204 0 : return C; /* normalized */
4205 : }
4206 :
4207 : /* p arbitrary */
4208 : GEN
4209 28236 : Flx_translate1_basecase(GEN P, ulong p)
4210 : {
4211 28236 : GEN R = Flx_copy(P);
4212 28236 : long i, k, n = degpol(P);
4213 654893 : for (i = 1; i <= n; i++)
4214 14846873 : for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
4215 28236 : return R;
4216 : }
4217 :
4218 : static int
4219 41401 : translate_basecase(long n, ulong p)
4220 : {
4221 : #ifdef LONG_IS_64BIT
4222 36102 : if (p <= 19) return n < 40;
4223 29910 : if (p < 1UL<<30) return n < 58;
4224 0 : if (p < 1UL<<59) return n < 100;
4225 0 : if (p < 1UL<<62) return n < 120;
4226 0 : if (p < 1UL<<63) return n < 240;
4227 0 : return n < 250;
4228 : #else
4229 5299 : if (p <= 13) return n < 18;
4230 4136 : if (p <= 17) return n < 22;
4231 4078 : if (p <= 29) return n < 39;
4232 3886 : if (p <= 67) return n < 69;
4233 3667 : if (p < 1UL<< 15) return n < 80;
4234 2047 : if (p < 1UL<< 16) return n < 100;
4235 0 : if (p < 1UL<< 28) return n < 300;
4236 0 : return n < 650;
4237 : #endif
4238 : }
4239 : /* assume p prime */
4240 : GEN
4241 16142 : Flx_translate1(GEN P, ulong p)
4242 : {
4243 16142 : long d, n = degpol(P);
4244 : GEN R, Q, S;
4245 16142 : if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
4246 : /* n > 0 */
4247 1148 : d = n >> 1;
4248 1148 : if ((ulong)n < p)
4249 : {
4250 0 : R = Flx_translate1(Flxn_red(P, d), p);
4251 0 : Q = Flx_translate1(Flx_shift(P, -d), p);
4252 0 : S = Fl_Xp1_powu(d, p, P[1]);
4253 0 : return Flx_add(Flx_mul(Q, S, p), R, p);
4254 : }
4255 : else
4256 : {
4257 : ulong q;
4258 1148 : if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
4259 1148 : R = Flx_translate1(Flxn_red(P, q), p);
4260 1148 : Q = Flx_translate1(Flx_shift(P, -q), p);
4261 1148 : S = Flx_add(Flx_shift(Q, q), Q, p);
4262 1148 : return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
4263 : }
4264 : }
4265 :
4266 : static GEN
4267 12017 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
4268 : {
4269 12017 : ulong k, d = n >> 1, c, v = 0;
4270 12017 : GEN C, V, W, U = upowers(p, e-1);
4271 12017 : init_invlaplace(d, p, &V, &W);
4272 12017 : Flv_inv_inplace(V, q);
4273 12017 : C = cgetg(n+3, t_VECSMALL);
4274 12017 : C[1] = vs;
4275 12017 : uel(C,2) = 1UL;
4276 12017 : uel(C,3) = n%q;
4277 12017 : v = u_lvalrem(n, p, &c);
4278 1355682 : for (k = 2; k <= d; k++)
4279 : {
4280 : ulong w;
4281 1343665 : v += u_lvalrem(n-k+1, p, &w) - W[k];
4282 1343665 : c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
4283 1343665 : uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
4284 : }
4285 1374521 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4286 12017 : return C; /* normalized */
4287 : }
4288 :
4289 : GEN
4290 25259 : zlx_translate1(GEN P, ulong p, long e)
4291 : {
4292 25259 : ulong d, q = upowuu(p,e), n = degpol(P);
4293 : GEN R, Q, S;
4294 25259 : if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
4295 : /* n > 0 */
4296 12017 : d = n >> 1;
4297 12017 : R = zlx_translate1(Flxn_red(P, d), p, e);
4298 12017 : Q = zlx_translate1(Flx_shift(P, -d), p, e);
4299 12017 : S = zl_Xp1_powu(d, p, q, e, P[1]);
4300 12017 : return Flx_add(Flx_mul(Q, S, q), R, q);
4301 : }
4302 :
4303 : /***********************************************************************/
4304 : /** Fl2 **/
4305 : /***********************************************************************/
4306 : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
4307 : * a nonsquare D. */
4308 :
4309 : INLINE GEN
4310 7195645 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
4311 :
4312 : /* allow pi = 0 */
4313 : GEN
4314 1917322 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4315 : {
4316 : ulong xaya, xbyb, Db2, mid, z1, z2;
4317 1917322 : ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
4318 1917322 : if (pi)
4319 : {
4320 1917336 : xaya = Fl_mul_pre(x1,y1,p,pi);
4321 1917766 : if (x2==0 && y2==0) return mkF2(xaya,0);
4322 1849039 : if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
4323 1824351 : if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
4324 1824044 : xbyb = Fl_mul_pre(x2,y2,p,pi);
4325 1823957 : mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
4326 1824193 : Db2 = Fl_mul_pre(D, xbyb, p,pi);
4327 : }
4328 0 : else if (p & HIGHMASK)
4329 : {
4330 0 : xaya = Fl_mul(x1,y1,p);
4331 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4332 0 : if (x2==0) return mkF2(xaya,Fl_mul(x1,y2,p));
4333 0 : if (y2==0) return mkF2(xaya,Fl_mul(x2,y1,p));
4334 0 : xbyb = Fl_mul(x2,y2,p);
4335 0 : mid = Fl_mul(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p);
4336 0 : Db2 = Fl_mul(D, xbyb, p);
4337 : }
4338 : else
4339 : {
4340 0 : xaya = (x1 * y1) % p;
4341 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4342 0 : if (x2==0) return mkF2(xaya, (x1 * y2) % p);
4343 0 : if (y2==0) return mkF2(xaya, (x2 * y1) % p);
4344 0 : xbyb = (x2 * y2) % p;
4345 0 : mid = (Fl_add(x1,x2,p) * Fl_add(y1,y2,p)) % p;
4346 0 : Db2 = (D * xbyb) % p;
4347 : }
4348 1824095 : z1 = Fl_add(xaya,Db2,p);
4349 1824065 : z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
4350 1824070 : return mkF2(z1,z2);
4351 : }
4352 :
4353 : /* allow pi = 0 */
4354 : GEN
4355 4825500 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
4356 : {
4357 4825500 : ulong a = x[1], b = x[2];
4358 : ulong a2, Db2, ab;
4359 4825500 : if (pi)
4360 : {
4361 4826212 : a2 = Fl_sqr_pre(a,p,pi);
4362 4828696 : if (b==0) return mkF2(a2,0);
4363 4616202 : Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
4364 4616276 : ab = Fl_mul_pre(a,b,p,pi);
4365 : }
4366 0 : else if (p & HIGHMASK)
4367 : {
4368 0 : a2 = Fl_sqr(a,p);
4369 0 : if (b==0) return mkF2(a2,0);
4370 0 : Db2= Fl_mul(D, Fl_sqr(b,p), p);
4371 0 : ab = Fl_mul(a,b,p);
4372 : }
4373 : else
4374 : {
4375 0 : a2 = (a * a) % p;
4376 0 : if (b==0) return mkF2(a2,0);
4377 0 : Db2= (D * ((b * b) % p)) % p;
4378 0 : ab = (a * b) % p;
4379 : }
4380 4615579 : return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
4381 : }
4382 :
4383 : /* allow pi = 0 */
4384 : ulong
4385 126016 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
4386 : {
4387 126016 : ulong a = x[1], b = x[2], a2;
4388 126016 : if (pi)
4389 : {
4390 72373 : a2 = Fl_sqr_pre(a,p,pi);
4391 72373 : return b? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p): a2;
4392 : }
4393 53643 : else if (p & HIGHMASK)
4394 : {
4395 0 : a2 = Fl_sqr(a,p);
4396 0 : return b? Fl_sub(a2, Fl_mul(D, Fl_sqr(b, p), p), p): a2;
4397 : }
4398 : else
4399 : {
4400 53643 : a2 = (a * a) % p;
4401 53643 : return b? Fl_sub(a2, (D * ((b * b) % p)) % p, p): a2;
4402 : }
4403 : }
4404 :
4405 : /* allow pi = 0 */
4406 : GEN
4407 192851 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
4408 : {
4409 192851 : ulong a = x[1], b = x[2], n, ni;
4410 192851 : if (b == 0) return mkF2(Fl_inv(a,p), 0);
4411 162245 : b = Fl_neg(b, p);
4412 162245 : if (pi)
4413 : {
4414 162245 : n = Fl_sub(Fl_sqr_pre(a, p,pi),
4415 : Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p);
4416 162248 : ni = Fl_inv(n,p);
4417 162250 : return mkF2(Fl_mul_pre(a, ni, p,pi), Fl_mul_pre(b, ni, p,pi));
4418 : }
4419 0 : else if (p & HIGHMASK)
4420 : {
4421 0 : n = Fl_sub(Fl_sqr(a, p), Fl_mul(D, Fl_sqr(b, p), p), p);
4422 0 : ni = Fl_inv(n,p);
4423 0 : return mkF2(Fl_mul(a, ni, p), Fl_mul(b, ni, p));
4424 : }
4425 : else
4426 : {
4427 0 : n = Fl_sub((a * a) % p, (D * ((b * b) % p)) % p, p);
4428 0 : ni = Fl_inv(n,p);
4429 0 : return mkF2((a * ni) % p, (b * ni) % p);
4430 : }
4431 : }
4432 :
4433 : int
4434 439646 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
4435 :
4436 : struct _Fl2 {
4437 : ulong p, pi, D;
4438 : };
4439 :
4440 : static GEN
4441 4826223 : _Fl2_sqr(void *data, GEN x)
4442 : {
4443 4826223 : struct _Fl2 *D = (struct _Fl2*)data;
4444 4826223 : return Fl2_sqr_pre(x, D->D, D->p, D->pi);
4445 : }
4446 : static GEN
4447 1888896 : _Fl2_mul(void *data, GEN x, GEN y)
4448 : {
4449 1888896 : struct _Fl2 *D = (struct _Fl2*)data;
4450 1888896 : return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
4451 : }
4452 :
4453 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL; allow pi = 0 */
4454 : GEN
4455 656843 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
4456 : {
4457 656843 : pari_sp av = avma;
4458 : struct _Fl2 d;
4459 : GEN y;
4460 656843 : long s = signe(n);
4461 656843 : if (!s) return mkF2(1,0);
4462 582173 : if (s < 0)
4463 192852 : x = Fl2_inv_pre(x,D,p,pi);
4464 582168 : if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
4465 428998 : d.p = p; d.pi = pi; d.D=D;
4466 428998 : y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
4467 429044 : return gerepileuptoleaf(av, y);
4468 : }
4469 :
4470 : static GEN
4471 656829 : _Fl2_pow(void *data, GEN x, GEN n)
4472 : {
4473 656829 : struct _Fl2 *D = (struct _Fl2*)data;
4474 656829 : return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
4475 : }
4476 :
4477 : static GEN
4478 111218 : _Fl2_rand(void *data)
4479 : {
4480 111218 : struct _Fl2 *D = (struct _Fl2*)data;
4481 111218 : ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
4482 111220 : return mkF2(a,b);
4483 : }
4484 :
4485 : GEN
4486 67676 : Fl2_sqrt_pre(GEN z, ulong D, ulong p, ulong pi)
4487 : {
4488 67676 : ulong a = uel(z,1), b = uel(z,2), as2, u, v, s;
4489 67676 : ulong y = Fl_2gener_pre_i(D, p, pi);
4490 67676 : if (b == 0)
4491 19383 : return krouu(a, p)==1 ? mkF2(Fl_sqrt_pre_i(a, y, p, pi), 0)
4492 19383 : : mkF2(0, Fl_sqrt_pre_i(Fl_div(a, D, p), y, p, pi));
4493 54473 : s = Fl_sqrt_pre_i(Fl2_norm_pre(z, D, p, pi), y, p, pi);
4494 54473 : if (s==~0UL) return NULL;
4495 51299 : as2 = Fl_halve(Fl_add(a, s, p), p);
4496 51299 : if (krouu(as2, p)==-1) as2 = Fl_sub(as2, s, p);
4497 51299 : u = Fl_sqrt_pre_i(as2, y, p, pi);
4498 51299 : v = Fl_div(b, Fl_double(u, p), p);
4499 51299 : return mkF2(u,v);
4500 : }
4501 :
4502 : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
4503 : hash_GEN, zv_equal, Fl2_equal1, NULL};
4504 :
4505 : /* allow pi = 0 */
4506 : GEN
4507 74669 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
4508 : {
4509 : struct _Fl2 E;
4510 : GEN o;
4511 74669 : if (a[1]==0 && a[2]==0)
4512 : {
4513 0 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
4514 0 : if (zeta) *zeta=mkF2(1,0);
4515 0 : return zv_copy(a);
4516 : }
4517 74669 : E.p=p; E.pi = pi; E.D = D;
4518 74669 : o = subiu(powuu(p,2), 1);
4519 74668 : return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
4520 : }
4521 :
4522 : /* allow pi = 0 */
4523 : GEN
4524 10528 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4525 : {
4526 : GEN p1;
4527 10528 : long i = lg(x)-1;
4528 10528 : if (i <= 2)
4529 2086 : return mkF2(i == 2? x[2]: 0, 0);
4530 8442 : p1 = mkF2(x[i], 0);
4531 36876 : for (i--; i>=2; i--)
4532 : {
4533 28434 : p1 = Fl2_mul_pre(p1, y, D, p, pi);
4534 28434 : uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
4535 : }
4536 8442 : return p1;
4537 : }
4538 :
4539 : /***********************************************************************/
4540 : /** FlxV **/
4541 : /***********************************************************************/
4542 : /* FlxV are t_VEC with Flx coefficients. */
4543 :
4544 : GEN
4545 34482 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
4546 : {
4547 34482 : pari_sp ltop=avma;
4548 : long i;
4549 34482 : GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
4550 257068 : for(i=2;i<lg(V);i++)
4551 222586 : z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
4552 34482 : return gerepileuptoleaf(ltop,z);
4553 : }
4554 :
4555 : GEN
4556 0 : ZXV_to_FlxV(GEN x, ulong p)
4557 0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
4558 :
4559 : GEN
4560 3799868 : ZXT_to_FlxT(GEN x, ulong p)
4561 : {
4562 3799868 : if (typ(x) == t_POL)
4563 3740666 : return ZX_to_Flx(x, p);
4564 : else
4565 194279 : pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
4566 : }
4567 :
4568 : GEN
4569 171947 : FlxV_to_Flm(GEN x, long n)
4570 927093 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
4571 :
4572 : GEN
4573 0 : FlxV_red(GEN x, ulong p)
4574 0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4575 :
4576 : GEN
4577 300208 : FlxT_red(GEN x, ulong p)
4578 : {
4579 300208 : if (typ(x) == t_VECSMALL)
4580 201954 : return Flx_red(x, p);
4581 : else
4582 329420 : pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4583 : }
4584 :
4585 : GEN
4586 113589 : FlxqV_dotproduct_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
4587 : {
4588 113589 : long i, lx = lg(x);
4589 : pari_sp av;
4590 : GEN c;
4591 113589 : if (lx == 1) return pol0_Flx(get_Flx_var(T));
4592 113589 : av = avma; c = Flx_mul_pre(gel(x,1),gel(y,1), p, pi);
4593 464499 : for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4594 113589 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4595 : }
4596 : GEN
4597 0 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4598 0 : { return FlxqV_dotproduct_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4599 :
4600 : GEN
4601 1918 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4602 : {
4603 1918 : long i, l = minss(lg(x), lg(y));
4604 : ulong pi;
4605 : pari_sp av;
4606 : GEN c;
4607 1918 : if (l == 2) return pol0_Flx(get_Flx_var(T));
4608 1905 : av = avma; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4609 1905 : c = Flx_mul_pre(gel(x,2),gel(y,2), p, pi);
4610 6202 : for (i=3; i<l; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4611 1905 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4612 : }
4613 :
4614 : /* allow pi = 0 */
4615 : GEN
4616 251071 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4617 : {
4618 251071 : long i, l = lg(z);
4619 251071 : GEN y = cgetg(l, t_VECSMALL);
4620 12741271 : for (i=1; i<l; i++) uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4621 251077 : return y;
4622 : }
4623 :
4624 : /***********************************************************************/
4625 : /** FlxM **/
4626 : /***********************************************************************/
4627 : /* allow pi = 0 */
4628 : GEN
4629 19452 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4630 : {
4631 19452 : long i, l = lg(z);
4632 19452 : GEN y = cgetg(l, t_MAT);
4633 270524 : for (i=1; i<l; i++) gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4634 19452 : return y;
4635 : }
4636 :
4637 : GEN
4638 0 : zero_FlxC(long n, long sv)
4639 : {
4640 0 : GEN x = cgetg(n + 1, t_COL), z = zero_Flx(sv);
4641 : long i;
4642 0 : for (i = 1; i <= n; i++) gel(x, i) = z;
4643 0 : return x;
4644 : }
4645 :
4646 : GEN
4647 0 : FlxC_neg(GEN x, ulong p)
4648 0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4649 :
4650 : GEN
4651 0 : FlxC_sub(GEN x, GEN y, ulong p)
4652 0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4653 :
4654 : GEN
4655 0 : zero_FlxM(long r, long c, long sv)
4656 : {
4657 0 : GEN x = cgetg(c + 1, t_MAT), z = zero_FlxC(r, sv);
4658 : long j;
4659 0 : for (j = 1; j <= c; j++) gel(x, j) = z;
4660 0 : return x;
4661 : }
4662 :
4663 : GEN
4664 0 : FlxM_neg(GEN x, ulong p)
4665 0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4666 :
4667 : GEN
4668 0 : FlxM_sub(GEN x, GEN y, ulong p)
4669 0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4670 :
4671 : GEN
4672 0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4673 0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4674 :
4675 : GEN
4676 0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4677 0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4678 :
4679 : static GEN
4680 47343 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4681 : long i, j, l, lc;
4682 47343 : GEN N = cgetg_copy(M, &l), x;
4683 47343 : if (l == 1)
4684 0 : return N;
4685 47343 : lc = lgcols(M);
4686 206352 : for (j = 1; j < l; j++) {
4687 159009 : gel(N, j) = cgetg(lc, t_COL);
4688 905441 : for (i = 1; i < lc; i++) {
4689 746432 : x = gcoeff(M, i, j);
4690 746432 : gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4691 : }
4692 : }
4693 47343 : return N;
4694 : }
4695 :
4696 : static GEN
4697 689211 : kron_pack_Flx_spec_half(GEN x, long l) {
4698 689211 : if (l == 0) return gen_0;
4699 458497 : return Flx_to_int_halfspec(x, l);
4700 : }
4701 :
4702 : static GEN
4703 53832 : kron_pack_Flx_spec(GEN x, long l) {
4704 : long i;
4705 : GEN w, y;
4706 53832 : if (l == 0)
4707 10072 : return gen_0;
4708 43760 : y = cgetipos(l + 2);
4709 159479 : for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4710 115719 : *w = x[i];
4711 43760 : return y;
4712 : }
4713 :
4714 : static GEN
4715 3389 : kron_pack_Flx_spec_2(GEN x, long l) { return Flx_eval2BILspec(x, 2, l); }
4716 :
4717 : static GEN
4718 0 : kron_pack_Flx_spec_3(GEN x, long l) { return Flx_eval2BILspec(x, 3, l); }
4719 :
4720 : static GEN
4721 42953 : kron_unpack_Flx(GEN z, ulong p)
4722 : {
4723 42953 : long i, l = lgefint(z);
4724 42953 : GEN x = cgetg(l, t_VECSMALL), w;
4725 201969 : for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4726 159016 : x[i] = ((ulong) *w) % p;
4727 42953 : return Flx_renormalize(x, l);
4728 : }
4729 :
4730 : static GEN
4731 2930 : kron_unpack_Flx_2(GEN x, ulong p) {
4732 2930 : long d = (lgefint(x)-1)/2 - 1;
4733 2930 : return Z_mod2BIL_Flx_2(x, d, p);
4734 : }
4735 :
4736 : static GEN
4737 0 : kron_unpack_Flx_3(GEN x, ulong p) {
4738 0 : long d = lgefint(x)/3 - 1;
4739 0 : return Z_mod2BIL_Flx_3(x, d, p);
4740 : }
4741 :
4742 : static GEN
4743 116095 : FlxM_pack_ZM_bits(GEN M, long b)
4744 : {
4745 : long i, j, l, lc;
4746 116095 : GEN N = cgetg_copy(M, &l), x;
4747 116095 : if (l == 1)
4748 0 : return N;
4749 116095 : lc = lgcols(M);
4750 478790 : for (j = 1; j < l; j++) {
4751 362695 : gel(N, j) = cgetg(lc, t_COL);
4752 5949368 : for (i = 1; i < lc; i++) {
4753 5586673 : x = gcoeff(M, i, j);
4754 5586673 : gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
4755 : }
4756 : }
4757 116095 : return N;
4758 : }
4759 :
4760 : static GEN
4761 23675 : ZM_unpack_FlxqM(GEN M, GEN T, ulong p, ulong pi, GEN (*unpack)(GEN, ulong))
4762 : {
4763 23675 : long i, j, l, lc, sv = get_Flx_var(T);
4764 23675 : GEN N = cgetg_copy(M, &l), x;
4765 23675 : if (l == 1)
4766 0 : return N;
4767 23675 : lc = lgcols(M);
4768 111660 : for (j = 1; j < l; j++) {
4769 87985 : gel(N, j) = cgetg(lc, t_COL);
4770 634641 : for (i = 1; i < lc; i++) {
4771 546656 : x = unpack(gcoeff(M, i, j), p);
4772 546656 : x[1] = sv;
4773 546656 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4774 : }
4775 : }
4776 23675 : return N;
4777 : }
4778 :
4779 : static GEN
4780 58088 : ZM_unpack_FlxqM_bits(GEN M, long b, GEN T, ulong p, ulong pi)
4781 : {
4782 58088 : long i, j, l, lc, sv = get_Flx_var(T);
4783 58088 : GEN N = cgetg_copy(M, &l), x;
4784 58088 : if (l == 1)
4785 0 : return N;
4786 58088 : lc = lgcols(M);
4787 58088 : if (b < BITS_IN_LONG) {
4788 194986 : for (j = 1; j < l; j++) {
4789 138551 : gel(N, j) = cgetg(lc, t_COL);
4790 3244565 : for (i = 1; i < lc; i++) {
4791 3106014 : x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
4792 3106014 : x[1] = sv;
4793 3106014 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4794 : }
4795 : }
4796 : } else {
4797 1653 : ulong pi = get_Fl_red(p);
4798 9784 : for (j = 1; j < l; j++) {
4799 8131 : gel(N, j) = cgetg(lc, t_COL);
4800 175175 : for (i = 1; i < lc; i++) {
4801 167044 : x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
4802 167044 : x[1] = sv;
4803 167044 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4804 : }
4805 : }
4806 : }
4807 58088 : return N;
4808 : }
4809 :
4810 : GEN
4811 81763 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
4812 : {
4813 81763 : pari_sp av = avma;
4814 81763 : long b, d = get_Flx_degree(T), n = lg(A) - 1;
4815 : GEN C, D, z;
4816 : GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
4817 81763 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4818 81763 : int is_sqr = A==B;
4819 :
4820 81763 : z = muliu(muliu(sqru(p - 1), d), n);
4821 81763 : b = expi(z) + 1;
4822 : /* only do expensive bit-packing if it saves at least 1 limb */
4823 81763 : if (b <= BITS_IN_HALFULONG)
4824 77394 : { if (nbits2nlong(d*b) == (d + 1)/2) b = BITS_IN_HALFULONG; }
4825 : else
4826 : {
4827 4369 : long l = lgefint(z) - 2;
4828 4369 : if (nbits2nlong(d*b) == d*l) b = l*BITS_IN_LONG;
4829 : }
4830 81763 : set_avma(av);
4831 :
4832 81763 : switch (b) {
4833 22608 : case BITS_IN_HALFULONG:
4834 22608 : pack = kron_pack_Flx_spec_half;
4835 22608 : unpack = int_to_Flx_half;
4836 22608 : break;
4837 1018 : case BITS_IN_LONG:
4838 1018 : pack = kron_pack_Flx_spec;
4839 1018 : unpack = kron_unpack_Flx;
4840 1018 : break;
4841 49 : case 2*BITS_IN_LONG:
4842 49 : pack = kron_pack_Flx_spec_2;
4843 49 : unpack = kron_unpack_Flx_2;
4844 49 : break;
4845 0 : case 3*BITS_IN_LONG:
4846 0 : pack = kron_pack_Flx_spec_3;
4847 0 : unpack = kron_unpack_Flx_3;
4848 0 : break;
4849 58088 : default:
4850 58088 : A = FlxM_pack_ZM_bits(A, b);
4851 58088 : B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
4852 58088 : C = ZM_mul(A, B);
4853 58088 : D = ZM_unpack_FlxqM_bits(C, b, T, p, pi);
4854 58088 : return gerepilecopy(av, D);
4855 : }
4856 23675 : A = FlxM_pack_ZM(A, pack);
4857 23675 : B = is_sqr? A: FlxM_pack_ZM(B, pack);
4858 23675 : C = ZM_mul(A, B);
4859 23675 : D = ZM_unpack_FlxqM(C, T, p, pi, unpack);
4860 23675 : return gerepilecopy(av, D);
4861 : }
|
