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 972659895 : get_Flx_red(GEN T, GEN *B)
22 : {
23 972659895 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 681538 : *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 37133584 : Flx_to_ZX(GEN z)
43 : {
44 37133584 : long i, l = lg(z);
45 37133584 : GEN x = cgetg(l,t_POL);
46 242492537 : for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
47 37120874 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
48 : }
49 :
50 : GEN
51 71396 : Flx_to_FlxX(GEN z, long sv)
52 : {
53 71396 : long i, l = lg(z);
54 71396 : GEN x = cgetg(l,t_POL);
55 278351 : for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
56 71399 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
57 : }
58 :
59 : /* same as Flx_to_ZX, in place */
60 : GEN
61 36455024 : Flx_to_ZX_inplace(GEN z)
62 : {
63 36455024 : long i, l = lg(z);
64 227412636 : for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
65 36450766 : 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 65863905 : Flx_to_Flv(GEN x, long N)
71 : {
72 65863905 : GEN z = cgetg(N+1,t_VECSMALL);
73 65858459 : long i, l = lg(x)-1;
74 65858459 : x++;
75 705364619 : for (i=1; i<l ; i++) z[i]=x[i];
76 328359788 : for ( ; i<=N; i++) z[i]=0;
77 65858459 : return z;
78 : }
79 :
80 : /*Flv_to_Flx=zv_to_zx*/
81 : GEN
82 25262798 : Flv_to_Flx(GEN x, long sv)
83 : {
84 25262798 : long i, l=lg(x)+1;
85 25262798 : GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
86 25258573 : x--;
87 278423956 : for (i=2; i<l ; i++) z[i]=x[i];
88 25258573 : 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 103996 : FlxC_to_ZXC(GEN x)
99 527283 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
100 :
101 : /*FlxC_to_ZXC=zxV_to_ZXV*/
102 : GEN
103 612675 : FlxV_to_ZXV(GEN x)
104 2478977 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
105 :
106 : void
107 2929057 : FlxV_to_ZXV_inplace(GEN v)
108 : {
109 : long i;
110 7779630 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
111 2928971 : }
112 :
113 : /*FlxM_to_ZXM=zxM_to_ZXM*/
114 : GEN
115 2429 : FlxM_to_ZXM(GEN x)
116 8183 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
117 :
118 : GEN
119 398465 : FlxV_to_FlxX(GEN x, long v)
120 : {
121 398465 : long i, l = lg(x)+1;
122 398465 : GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
123 398465 : x--;
124 5000952 : for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
125 398465 : 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 19856089 : Fl_to_Flx(ulong x, long sv) { return x? mkvecsmall2(sv, x): pol0_Flx(sv); }
156 :
157 : /* a X^d */
158 : GEN
159 913933 : monomial_Flx(ulong a, long d, long vs)
160 : {
161 : GEN P;
162 913933 : if (a==0) return pol0_Flx(vs);
163 913933 : P = const_vecsmall(d+2, 0);
164 913945 : P[1] = vs; P[d+2] = a; return P;
165 : }
166 :
167 : GEN
168 2597811 : Z_to_Flx(GEN x, ulong p, long sv)
169 : {
170 2597811 : long u = umodiu(x,p);
171 2597811 : 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 167522649 : ZX_to_Flx(GEN x, ulong p)
177 : {
178 167522649 : long i, lx = lg(x);
179 167522649 : GEN a = cgetg(lx, t_VECSMALL);
180 167484515 : a[1]=((ulong)x[1])&VARNBITS;
181 1111309001 : for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
182 167487652 : 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 6070885 : zx_to_Flx(GEN x, ulong p)
188 : {
189 6070885 : long i, lx = lg(x);
190 6070885 : GEN a = cgetg(lx, t_VECSMALL);
191 6065115 : a[1] = x[1];
192 18614324 : for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
193 6064191 : return Flx_renormalize(a,lx);
194 : }
195 :
196 : ulong
197 73773268 : Rg_to_Fl(GEN x, ulong p)
198 : {
199 73773268 : switch(typ(x))
200 : {
201 48778822 : case t_INT: return umodiu(x, p);
202 454332 : case t_FRAC: {
203 454332 : ulong z = umodiu(gel(x,1), p);
204 454334 : if (!z) return 0;
205 444572 : return Fl_div(z, umodiu(gel(x,2), p), p);
206 : }
207 205965 : case t_PADIC: return padic_to_Fl(x, p);
208 24334159 : case t_INTMOD: {
209 24334159 : GEN q = gel(x,1), a = gel(x,2);
210 24334159 : 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 1706761 : Rg_to_F2(GEN x)
221 : {
222 1706761 : switch(typ(x))
223 : {
224 273952 : case t_INT: return mpodd(x);
225 0 : case t_FRAC:
226 0 : if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
227 0 : return mpodd(gel(x,1));
228 0 : case t_PADIC:
229 0 : if (!absequaliu(padic_p(x),2)) pari_err_OP("",x, mkintmodu(1,2));
230 0 : if (valp(x) < 0) (void)Fl_inv(0,2);
231 0 : return valp(x) & 1;
232 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 2355069 : RgX_to_Flx(GEN x, ulong p)
244 : {
245 2355069 : long i, lx = lg(x);
246 2355069 : GEN a = cgetg(lx, t_VECSMALL);
247 2355069 : a[1]=((ulong)x[1])&VARNBITS;
248 20435241 : for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
249 2355069 : 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 3567813 : Rg_to_Flxq(GEN x, GEN T, ulong p)
259 : {
260 3567813 : long ta, tx = typ(x), v = get_Flx_var(T);
261 : ulong pi;
262 : GEN a, b;
263 3567810 : if (is_const_t(tx))
264 : {
265 3317038 : if (tx == t_FFELT) return FF_to_Flxq(x);
266 2586030 : return Fl_to_Flx(Rg_to_Fl(x, p), v);
267 : }
268 250772 : 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 242196 : case t_POL:
280 242196 : x = RgX_to_Flx(x,p);
281 242196 : if (x[1] != v) break;
282 242196 : 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 2108015237 : Flx_renormalize(GEN /*in place*/ x, long lx)
298 : {
299 : long i;
300 2356873867 : for (i = lx-1; i>1; i--)
301 2262361854 : if (x[i]) break;
302 2108015237 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
303 2106969664 : setlg(x, i+1); return x;
304 : }
305 :
306 : GEN
307 1877101 : Flx_red(GEN z, ulong p)
308 : {
309 1877101 : long i, l = lg(z);
310 1877101 : GEN x = cgetg(l, t_VECSMALL);
311 1876955 : x[1] = z[1];
312 33193492 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
313 1876955 : return Flx_renormalize(x,l);
314 : }
315 :
316 : int
317 27283519 : Flx_equal(GEN V, GEN W)
318 : {
319 27283519 : long l = lg(V);
320 27283519 : if (lg(W) != l) return 0;
321 28280486 : while (--l > 1) /* do not compare variables, V[1] */
322 27176462 : if (V[l] != W[l]) return 0;
323 1104024 : return 1;
324 : }
325 :
326 : GEN
327 2593073 : random_Flx(long d1, long vs, ulong p)
328 : {
329 2593073 : long i, d = d1+2;
330 2593073 : GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
331 17938808 : for (i=2; i<d; i++) y[i] = random_Fl(p);
332 2593302 : return Flx_renormalize(y,d);
333 : }
334 :
335 : static GEN
336 7163592 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
337 : {
338 : long i,lz;
339 : GEN z;
340 :
341 7163592 : if (ly>lx) swapspec(x,y, lx,ly);
342 7163592 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
343 106303861 : for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
344 90019447 : for ( ; i<lx; i++) z[i+2] = x[i];
345 7163592 : z[1] = 0; return Flx_renormalize(z, lz);
346 : }
347 :
348 : GEN
349 62606003 : Flx_add(GEN x, GEN y, ulong p)
350 : {
351 : long i,lz;
352 : GEN z;
353 62606003 : long lx=lg(x);
354 62606003 : long ly=lg(y);
355 62606003 : if (ly>lx) swapspec(x,y, lx,ly);
356 62606003 : lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
357 572562836 : for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
358 127968690 : for ( ; i<lx; i++) z[i] = x[i];
359 62582000 : return Flx_renormalize(z, lz);
360 : }
361 :
362 : GEN
363 9943291 : Flx_Fl_add(GEN y, ulong x, ulong p)
364 : {
365 : GEN z;
366 : long lz, i;
367 9943291 : if (!lgpol(y))
368 229689 : return Fl_to_Flx(x,y[1]);
369 9714925 : lz=lg(y);
370 9714925 : z=cgetg(lz,t_VECSMALL);
371 9714802 : z[1]=y[1];
372 9714802 : z[2] = Fl_add(y[2],x,p);
373 46962204 : for(i=3;i<lz;i++)
374 37247846 : z[i] = y[i];
375 9714358 : if (lz==3) z = Flx_renormalize(z,lz);
376 9714045 : return z;
377 : }
378 :
379 : static GEN
380 896706 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
381 : {
382 : long i,lz;
383 : GEN z;
384 :
385 896706 : if (ly <= lx)
386 : {
387 896716 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
388 53723084 : for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
389 1446970 : 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 896371 : z[1] = 0; return Flx_renormalize(z, lz);
398 : }
399 :
400 : GEN
401 137976139 : Flx_sub(GEN x, GEN y, ulong p)
402 : {
403 137976139 : long i,lz,lx = lg(x), ly = lg(y);
404 : GEN z;
405 :
406 137976139 : if (ly <= lx)
407 : {
408 87946646 : lz = lx; z = cgetg(lz, t_VECSMALL);
409 456294570 : for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
410 175826410 : for ( ; i<lx; i++) z[i] = x[i];
411 : }
412 : else
413 : {
414 50029493 : lz = ly; z = cgetg(lz, t_VECSMALL);
415 259590707 : for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
416 232106757 : for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
417 : }
418 137966525 : z[1]=x[1]; return Flx_renormalize(z, lz);
419 : }
420 :
421 : GEN
422 151422 : Flx_Fl_sub(GEN y, ulong x, ulong p)
423 : {
424 : GEN z;
425 151422 : long lz = lg(y), i;
426 151422 : if (lz==2)
427 513 : return Fl_to_Flx(Fl_neg(x, p),y[1]);
428 150909 : z = cgetg(lz, t_VECSMALL);
429 150909 : z[1] = y[1];
430 150909 : z[2] = Fl_sub(uel(y,2), x, p);
431 751873 : for(i=3; i<lz; i++)
432 600964 : z[i] = y[i];
433 150909 : if (lz==3) z = Flx_renormalize(z,lz);
434 150909 : return z;
435 : }
436 :
437 : static GEN
438 3264301 : Flx_negspec(GEN x, ulong p, long l)
439 : {
440 : long i;
441 3264301 : GEN z = cgetg(l+2, t_VECSMALL) + 2;
442 20973225 : for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
443 3264227 : return z-2;
444 : }
445 :
446 : GEN
447 3264312 : Flx_neg(GEN x, ulong p)
448 : {
449 3264312 : GEN z = Flx_negspec(x+2, p, lgpol(x));
450 3264452 : z[1] = x[1];
451 3264452 : return z;
452 : }
453 :
454 : GEN
455 1755948 : Flx_neg_inplace(GEN x, ulong p)
456 : {
457 1755948 : long i, l = lg(x);
458 52206744 : for (i=2; i<l; i++)
459 50450796 : if (x[i]) x[i] = p - x[i];
460 1755948 : return x;
461 : }
462 :
463 : GEN
464 2444279 : Flx_double(GEN y, ulong p)
465 : {
466 : long i, l;
467 2444279 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
468 20388866 : for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
469 2444279 : return Flx_renormalize(z, l);
470 : }
471 : GEN
472 1049257 : Flx_triple(GEN y, ulong p)
473 : {
474 : long i, l;
475 1049257 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
476 8305948 : for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
477 1049257 : return Flx_renormalize(z, l);
478 : }
479 :
480 : GEN
481 18257131 : Flx_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
482 : {
483 : GEN z;
484 : long i, l;
485 18257131 : if (!x) return pol0_Flx(y[1]);
486 17474498 : z = cgetg_copy(y, &l); z[1] = y[1];
487 17474337 : if (pi==0)
488 : {
489 15292965 : if (HIGHWORD(x | p))
490 0 : for(i=2; i<l; i++) z[i] = Fl_mul(uel(y,i), x, p);
491 : else
492 92227090 : for(i=2; i<l; i++) z[i] = (uel(y,i) * x) % p;
493 : } else
494 17963631 : for(i=2; i<l; i++) z[i] = Fl_mul_pre(uel(y,i), x, p, pi);
495 17475642 : return Flx_renormalize(z, l);
496 : }
497 :
498 : GEN
499 7322154 : Flx_Fl_mul(GEN x, ulong y, ulong p)
500 7322154 : { 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 11961681 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
515 : {
516 : GEN z;
517 : long i, l;
518 11961681 : z = cgetg_copy(y, &l); z[1] = y[1];
519 11958638 : if (HIGHWORD(x | p))
520 5408854 : for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
521 : else
522 26826723 : for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
523 11958633 : 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 26810107 : Flx_shift(GEN a, long n)
529 : {
530 26810107 : long i, l = lg(a);
531 : GEN b;
532 26810107 : if (l==2 || !n) return Flx_copy(a);
533 26467183 : if (l+n<=2) return pol0_Flx(a[1]);
534 26252612 : b = cgetg(l+n, t_VECSMALL);
535 26250595 : b[1] = a[1];
536 26250595 : if (n < 0)
537 71656106 : for (i=2-n; i<l; i++) b[i+n] = a[i];
538 : else
539 : {
540 50897245 : for (i=0; i<n; i++) b[2+i] = 0;
541 148278553 : for (i=2; i<l; i++) b[i+n] = a[i];
542 : }
543 26250595 : return b;
544 : }
545 :
546 : GEN
547 62158662 : Flx_normalize(GEN z, ulong p)
548 : {
549 62158662 : long l = lg(z)-1;
550 62158662 : ulong p1 = z[l]; /* leading term */
551 62158662 : if (p1 == 1) return z;
552 11941020 : 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 3680970 : Flx_addshift(GEN x, GEN y, ulong p, long d)
558 : {
559 3680970 : GEN xd,yd,zd = (GEN)avma;
560 3680970 : long a,lz,ny = lgpol(y), nx = lgpol(x);
561 3680970 : long vs = x[1];
562 3680970 : if (nx == 0) return y;
563 3679118 : x += 2; y += 2; a = ny-d;
564 3679118 : if (a <= 0)
565 : {
566 84992 : lz = (a>nx)? ny+2: nx+d+2;
567 84992 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
568 1726350 : while (xd > x) *--zd = *--xd;
569 84992 : x = zd + a;
570 163542 : while (zd > x) *--zd = 0;
571 : }
572 : else
573 : {
574 3594126 : xd = new_chunk(d); yd = y+d;
575 3594126 : x = Flx_addspec(x,yd,p, nx,a);
576 3594126 : lz = (a>nx)? ny+2: lg(x)+d;
577 132376043 : x += 2; while (xd > x) *--zd = *--xd;
578 : }
579 60201202 : while (yd > y) *--zd = *--yd;
580 3679118 : *--zd = vs;
581 3679118 : *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
582 : }
583 :
584 : /* shift polynomial + gerepile */
585 : /* Do not set evalvarn*/
586 : static GEN
587 625442639 : Flx_shiftip(pari_sp av, GEN x, long v)
588 : {
589 625442639 : long i, lx = lg(x), ly;
590 : GEN y;
591 625442639 : if (!v || lx==2) return gerepileuptoleaf(av, x);
592 174037219 : ly = lx + v; /* result length */
593 174037219 : (void)new_chunk(ly); /* check that result fits */
594 173941009 : x += lx; y = (GEN)av;
595 1231815075 : for (i = 2; i<lx; i++) *--y = *--x;
596 701119206 : for (i = 0; i< v; i++) *--y = 0;
597 173941009 : y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
598 174084392 : return gc_const((pari_sp)y, y);
599 : }
600 :
601 : static long
602 2289348907 : get_Fl_threshold(ulong p, long mul, long mul2)
603 : {
604 2289348907 : 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 8336640 : maxbitcoeffpol(ulong p, long n)
615 : {
616 8336640 : GEN z = muliu(sqru(p - 1), n);
617 8334087 : long b = expi(z) + 1;
618 : /* only do expensive bit-packing if it saves at least 1 limb */
619 8334775 : if (b <= BITS_IN_QUARTULONG)
620 : {
621 873145 : if (nbits2nlong(n*b) == (n + 3)>>2)
622 107501 : b = BITS_IN_QUARTULONG;
623 : }
624 7461630 : else if (b <= BITS_IN_HALFULONG)
625 : {
626 1545662 : if (nbits2nlong(n*b) == (n + 1)>>1)
627 5590 : b = BITS_IN_HALFULONG;
628 : }
629 : else
630 : {
631 5915968 : long l = lgefint(z) - 2;
632 5915968 : if (nbits2nlong(n*b) == n*l)
633 307341 : b = l*BITS_IN_LONG;
634 : }
635 8334544 : return b;
636 : }
637 :
638 : INLINE ulong
639 3363145485 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
640 : { /* Assume OK_ULONG*/
641 3363145485 : ulong p1 = 0;
642 : long i;
643 15938774852 : for (i=a; i<b; i++)
644 12575629367 : if (y[i])
645 : {
646 10566221532 : p1 += y[i] * x[-i];
647 10566221532 : if (p1 & HIGHBIT) p1 %= p;
648 : }
649 3363145485 : return p1 % p;
650 : }
651 :
652 : INLINE ulong
653 1164926141 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
654 : {
655 1164926141 : ulong p1 = 0;
656 : long i;
657 3693282796 : for (i=a; i<b; i++)
658 2529104963 : if (y[i])
659 2489687661 : p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
660 1164177833 : return p1;
661 : }
662 :
663 : /* assume nx >= ny > 0 */
664 : static GEN
665 338423131 : 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 338423131 : lz = nx+ny+1; nz = lz-2;
671 338423131 : z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
672 338219813 : if (!pi)
673 : {
674 1133211147 : for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
675 725858242 : for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
676 882539907 : for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
677 : }
678 : else
679 : {
680 308077387 : for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
681 211403335 : for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
682 222039304 : for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
683 : }
684 338134645 : 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 1060205 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
693 12298 : 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 10039 : return int_to_Flx(z,p);
701 : }
702 :
703 : static GEN
704 470235 : Flx_to_int_halfspec(GEN a, long na)
705 : {
706 : long j;
707 470235 : long n = (na+1)>>1UL;
708 470235 : GEN V = cgetipos(2+n);
709 : GEN w;
710 1380211 : for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
711 909976 : *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
712 470235 : if (j<na)
713 319584 : *w = a[j];
714 470235 : return V;
715 : }
716 :
717 : static GEN
718 507447 : int_to_Flx_half(GEN z, ulong p)
719 : {
720 : long i;
721 507447 : long lx = (lgefint(z)-2)*2+2;
722 507447 : GEN w, x = cgetg(lx, t_VECSMALL);
723 1912035 : for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
724 : {
725 1404588 : x[i] = LOWWORD((ulong)*w)%p;
726 1404588 : x[i+1] = HIGHWORD((ulong)*w)%p;
727 : }
728 507447 : 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 204783 : Flx_to_int_quartspec(GEN a, long na)
742 : {
743 : long j;
744 204783 : long n = (na+3)>>2UL;
745 204783 : GEN V = cgetipos(2+n);
746 : GEN w;
747 4380273 : for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
748 4175488 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
749 204785 : switch (na-j)
750 : {
751 116384 : case 3:
752 116384 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
753 116384 : break;
754 34458 : case 2:
755 34458 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
756 34458 : break;
757 27371 : case 1:
758 27371 : *w = a[j];
759 27371 : break;
760 26575 : case 0:
761 26575 : break;
762 : }
763 204785 : return V;
764 : }
765 :
766 : static GEN
767 107504 : int_to_Flx_quart(GEN z, ulong p)
768 : {
769 : long i;
770 107504 : long lx = (lgefint(z)-2)*4+2;
771 107504 : GEN w, x = cgetg(lx, t_VECSMALL);
772 4876315 : for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
773 : {
774 4768811 : x[i] = LLQUARTWORD((ulong)*w)%p;
775 4768811 : x[i+1] = HLQUARTWORD((ulong)*w)%p;
776 4768811 : x[i+2] = LHQUARTWORD((ulong)*w)%p;
777 4768811 : x[i+3] = HHQUARTWORD((ulong)*w)%p;
778 : }
779 107504 : return Flx_renormalize(x, lx);
780 : }
781 :
782 : static GEN
783 97283 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
784 : {
785 97283 : GEN A = Flx_to_int_quartspec(a,na);
786 97283 : GEN B = Flx_to_int_quartspec(b,nb);
787 97286 : GEN z = mulii(A,B);
788 97286 : 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 581977 : Flx_eval2BILspec(GEN x, long k, long l)
794 : {
795 581977 : long i, lz = k*l, ki;
796 581977 : GEN pz = cgetipos(2+lz);
797 16363781 : for (i=0; i < lz; i++)
798 15781804 : *int_W(pz,i) = 0UL;
799 8472879 : for (i=0, ki=0; i<l; i++, ki+=k)
800 7890902 : *int_W(pz,ki) = x[i];
801 581977 : return int_normalize(pz,0);
802 : }
803 :
804 : static GEN
805 297973 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
806 : {
807 297973 : long i, offset, lm = lgefint(x)-2, l = d+3;
808 297973 : ulong pi = get_Fl_red(p);
809 297973 : GEN pol = cgetg(l, t_VECSMALL);
810 297973 : pol[1] = 0;
811 8007356 : for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
812 7709383 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
813 297973 : if (offset < lm)
814 225080 : pol[i+2] = (*int_W(x,offset)) % p;
815 297973 : 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 295043 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
837 : {
838 295043 : return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
839 : }
840 :
841 : static GEN
842 283545 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
843 : {
844 283545 : pari_sp av = avma;
845 283545 : GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
846 283545 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
847 : }
848 :
849 : static GEN
850 20736700 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
851 : GEN y;
852 : long i;
853 20736700 : if (l == 0)
854 3429546 : return gen_0;
855 17307154 : y = cgetg(l + 1, t_VECSMALL);
856 812302504 : for(i = 1; i <= l; i++)
857 795000202 : y[i] = x[l - i];
858 17302302 : return nv_fromdigits_2k(y, b);
859 : }
860 :
861 : /* assume b < BITS_IN_LONG */
862 : static GEN
863 5646458 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
864 5646458 : GEN v = binary_2k_nv(z, b), x;
865 5646498 : long i, l = lg(v) + 1;
866 5646498 : x = cgetg(l, t_VECSMALL);
867 620454224 : for (i = 2; i < l; i++)
868 614807668 : x[i] = v[l - i] % p;
869 5646556 : return Flx_renormalize(x, l);
870 : }
871 :
872 : static GEN
873 5548339 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
874 5548339 : GEN v = binary_2k(z, b), x, y;
875 5548731 : long i, l = lg(v) + 1, ly;
876 5548731 : x = cgetg(l, t_VECSMALL);
877 232941189 : for (i = 2; i < l; i++) {
878 227395341 : y = gel(v, l - i);
879 227395341 : ly = lgefint(y);
880 227395341 : switch (ly) {
881 6285308 : case 2: x[i] = 0; break;
882 29319708 : case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
883 175895169 : case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
884 31789857 : case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
885 15895156 : *int_W_lg(y, 0, ly), p, pi); break;
886 0 : default: x[i] = umodiu(gel(v, l - i), p);
887 : }
888 : }
889 5545848 : return Flx_renormalize(x, l);
890 : }
891 :
892 : static GEN
893 7229308 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
894 : {
895 : GEN C, D;
896 7229308 : pari_sp av = avma;
897 7229308 : A = kron_pack_Flx_spec_bits(A, b, lA);
898 7235546 : B = kron_pack_Flx_spec_bits(B, b, lB);
899 7235677 : C = gerepileuptoint(av, mulii(A, B));
900 7233908 : if (b < BITS_IN_LONG)
901 2059082 : D = kron_unpack_Flx_bits_narrow(C, b, p);
902 : else
903 : {
904 5174826 : ulong pi = get_Fl_red(p);
905 5173503 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
906 : }
907 7231752 : return D;
908 : }
909 :
910 : static GEN
911 683763 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
912 : {
913 : GEN C, D;
914 683763 : A = kron_pack_Flx_spec_bits(A, b, lA);
915 683848 : C = sqri(A);
916 683862 : if (b < BITS_IN_LONG)
917 475718 : D = kron_unpack_Flx_bits_narrow(C, b, p);
918 : else
919 : {
920 208144 : ulong pi = get_Fl_red(p);
921 208139 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
922 : }
923 683839 : 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 375715596 : Flx_mulspec(GEN a, GEN b, ulong p, ulong pi, long na, long nb)
932 : {
933 : GEN a0,c,c0;
934 375715596 : long n0, n0a, i, v = 0;
935 : pari_sp av;
936 :
937 479954840 : while (na && !a[0]) { a++; na--; v++; }
938 560326244 : while (nb && !b[0]) { b++; nb--; v++; }
939 375715596 : if (na < nb) swapspec(a,b, na,nb);
940 375715596 : if (!nb) return pol0_Flx(0);
941 :
942 347578143 : av = avma;
943 347578143 : if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
944 : {
945 7629144 : long m = maxbitcoeffpol(p,nb);
946 7625385 : switch (m)
947 : {
948 97283 : case BITS_IN_QUARTULONG:
949 97283 : 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 283545 : case 2*BITS_IN_LONG:
955 283545 : 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 7229069 : default:
959 7229069 : return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
960 : }
961 : }
962 340135722 : if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
963 338351770 : return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,pi,na,nb), v);
964 1809863 : i=(na>>1); n0=na-i; na=i;
965 1809863 : a0=a+n0; n0a=n0;
966 2575665 : while (n0a && !a[n0a-1]) n0a--;
967 :
968 1809863 : if (nb > n0)
969 : {
970 : GEN b0,c1,c2;
971 : long n0b;
972 :
973 1755948 : nb -= n0; b0 = b+n0; n0b = n0;
974 2835403 : while (n0b && !b[n0b-1]) n0b--;
975 1755948 : c = Flx_mulspec(a,b,p,pi,n0a,n0b);
976 1755948 : c0 = Flx_mulspec(a0,b0,p,pi,na,nb);
977 :
978 1755948 : c2 = Flx_addspec(a0,a,p,na,n0a);
979 1755948 : c1 = Flx_addspec(b0,b,p,nb,n0b);
980 :
981 1755948 : c1 = Flx_mul_pre(c1,c2,p,pi);
982 1755948 : c2 = Flx_add(c0,c,p);
983 :
984 1755948 : c2 = Flx_neg_inplace(c2,p);
985 1755948 : c2 = Flx_add(c1,c2,p);
986 1755948 : c0 = Flx_addshift(c0,c2 ,p, n0);
987 : }
988 : else
989 : {
990 53915 : c = Flx_mulspec(a,b,p,pi,n0a,nb);
991 53915 : c0 = Flx_mulspec(a0,b,p,pi,na,nb);
992 : }
993 1809863 : c0 = Flx_addshift(c0,c,p,n0);
994 1809863 : return Flx_shiftip(av,c0, v);
995 : }
996 :
997 : GEN
998 370084521 : Flx_mul_pre(GEN x, GEN y, ulong p, ulong pi)
999 : {
1000 370084521 : GEN z = Flx_mulspec(x+2,y+2,p, pi, lgpol(x),lgpol(y));
1001 370273308 : z[1] = x[1]; return z;
1002 : }
1003 : GEN
1004 27680786 : Flx_mul(GEN x, GEN y, ulong p)
1005 27680786 : { return Flx_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1006 :
1007 : static GEN
1008 278097900 : 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 278097900 : if (!nx) return pol0_Flx(0);
1015 278097900 : lz = (nx << 1) + 1, nz = lz-2;
1016 278097900 : z = cgetg(lz, t_VECSMALL) + 2;
1017 277639005 : if (!pi)
1018 : {
1019 213531522 : z[0] = x[0]*x[0]%p;
1020 913803033 : for (i=1; i<nx; i++)
1021 : {
1022 700389728 : p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1023 700271511 : p1 <<= 1;
1024 700271511 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1025 700271511 : z[i] = p1 % p;
1026 : }
1027 917857033 : for ( ; i<nz; i++)
1028 : {
1029 703822151 : p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1030 704443728 : p1 <<= 1;
1031 704443728 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1032 704443728 : z[i] = p1 % p;
1033 : }
1034 : }
1035 : else
1036 : {
1037 64107483 : z[0] = Fl_sqr_pre(x[0], p, pi);
1038 411072568 : for (i=1; i<nx; i++)
1039 : {
1040 347059900 : p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1041 347183252 : p1 = Fl_add(p1, p1, p);
1042 346715728 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1043 346908171 : z[i] = p1;
1044 : }
1045 411183186 : for ( ; i<nz; i++)
1046 : {
1047 347104434 : p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1048 347698707 : p1 = Fl_add(p1, p1, p);
1049 347308811 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1050 347170518 : z[i] = p1;
1051 : }
1052 : }
1053 278113634 : z -= 2; return Flx_renormalize(z, lz);
1054 : }
1055 :
1056 : static GEN
1057 2264 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
1058 : {
1059 2264 : 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 10218 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1072 : {
1073 10218 : GEN z = sqri(Flx_to_int_quartspec(a,na));
1074 10218 : return int_to_Flx_quart(z,p);
1075 : }
1076 :
1077 : static GEN
1078 11498 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1079 : {
1080 11498 : pari_sp av = avma;
1081 11498 : GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1082 11498 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1083 : }
1084 :
1085 : static GEN
1086 278290630 : Flx_sqrspec(GEN a, ulong p, ulong pi, long na)
1087 : {
1088 : GEN a0, c, c0;
1089 278290630 : long n0, n0a, i, v = 0, m;
1090 : pari_sp av;
1091 :
1092 399328330 : while (na && !a[0]) { a++; na--; v += 2; }
1093 278290630 : if (!na) return pol0_Flx(0);
1094 :
1095 278045307 : av = avma;
1096 278045307 : if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1097 : {
1098 707847 : m = maxbitcoeffpol(p,na);
1099 707873 : switch(m)
1100 : {
1101 10218 : case BITS_IN_QUARTULONG:
1102 10218 : 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 2264 : case BITS_IN_LONG:
1106 2264 : return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1107 11498 : case 2*BITS_IN_LONG:
1108 11498 : 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 683757 : default:
1112 683757 : return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1113 : }
1114 : }
1115 277773898 : if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1116 277732988 : return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,pi,na), v);
1117 57589 : i=(na>>1); n0=na-i; na=i;
1118 57589 : a0=a+n0; n0a=n0;
1119 72329 : while (n0a && !a[n0a-1]) n0a--;
1120 :
1121 57589 : c = Flx_sqrspec(a,p,pi,n0a);
1122 57589 : c0= Flx_sqrspec(a0,p,pi,na);
1123 57589 : if (p == 2) n0 *= 2;
1124 : else
1125 : {
1126 57570 : GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1127 57570 : t = Flx_sqr_pre(t,p,pi);
1128 57570 : c1= Flx_add(c0,c, p);
1129 57570 : c1= Flx_sub(t, c1, p);
1130 57570 : c0 = Flx_addshift(c0,c1,p,n0);
1131 : }
1132 57589 : c0 = Flx_addshift(c0,c,p,n0);
1133 57589 : return Flx_shiftip(av,c0,v);
1134 : }
1135 :
1136 : GEN
1137 278070888 : Flx_sqr_pre(GEN x, ulong p, ulong pi)
1138 : {
1139 278070888 : GEN z = Flx_sqrspec(x+2,p, pi, lgpol(x));
1140 279328986 : z[1] = x[1]; return z;
1141 : }
1142 : GEN
1143 356222 : Flx_sqr(GEN x, ulong p)
1144 356222 : { return Flx_sqr_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1145 :
1146 : GEN
1147 7778 : Flx_powu_pre(GEN x, ulong n, ulong p, ulong pi)
1148 : {
1149 7778 : GEN y = pol1_Flx(x[1]), z;
1150 : ulong m;
1151 7777 : if (n == 0) return y;
1152 7777 : m = n; z = x;
1153 : for (;;)
1154 : {
1155 30000 : if (m&1UL) y = Flx_mul_pre(y,z, p, pi);
1156 30004 : m >>= 1; if (!m) return y;
1157 22227 : 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 14223 : Flx_halve(GEN y, ulong p)
1169 : {
1170 : GEN z;
1171 : long i, l;
1172 14223 : z = cgetg_copy(y, &l); z[1] = y[1];
1173 59438 : 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 7126814 : Flx_recipspec(GEN x, long l, long n)
1179 : {
1180 : long i;
1181 7126814 : GEN z=cgetg(n+2,t_VECSMALL)+2;
1182 115434656 : for(i=0; i<l; i++)
1183 108309503 : z[n-i-1] = x[i];
1184 15594518 : for( ; i<n; i++)
1185 8469365 : z[n-i-1] = 0;
1186 7125153 : 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 P(x * h) */
1198 : GEN
1199 0 : Flx_unscale(GEN P, ulong h, ulong p)
1200 : {
1201 : long i, l;
1202 0 : ulong hi = 1UL;
1203 0 : GEN Q = cgetg_copy(P, &l);
1204 0 : Q[1] = P[1];
1205 0 : if (l == 2) return Q;
1206 0 : uel(Q,2) = uel(P,2);
1207 0 : for (i=3; i<l; i++)
1208 : {
1209 0 : hi = Fl_mul(hi, h ,p);
1210 0 : uel(Q,i) = Fl_mul(uel(P,i), hi, p);
1211 : }
1212 0 : return Q;
1213 : }
1214 : /* Return h^degpol(P) P(x / h) */
1215 : GEN
1216 1117 : Flx_rescale(GEN P, ulong h, ulong p)
1217 : {
1218 1117 : long i, l = lg(P);
1219 1117 : GEN Q = cgetg(l,t_VECSMALL);
1220 1117 : ulong hi = h;
1221 1117 : Q[l-1] = P[l-1];
1222 12538 : for (i=l-2; i>=2; i--)
1223 : {
1224 12538 : Q[i] = Fl_mul(P[i], hi, p);
1225 12538 : if (i == 2) break;
1226 11421 : hi = Fl_mul(hi,h, p);
1227 : }
1228 1117 : Q[1] = P[1]; return Q;
1229 : }
1230 :
1231 : /* x/polrecip(P)+O(x^n); allow pi = 0 */
1232 : static GEN
1233 134238 : Flx_invBarrett_basecase(GEN T, ulong p, ulong pi)
1234 : {
1235 134238 : long i, l=lg(T)-1, lr=l-1, k;
1236 134238 : GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1237 134238 : r[2] = 1;
1238 134238 : if (!pi)
1239 764563 : for (i=3;i<lr;i++)
1240 : {
1241 757568 : ulong u = uel(T, l-i+2);
1242 45415241 : for (k=3; k<i; k++)
1243 44657673 : { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1244 757568 : r[i] = Fl_neg(u % p, p);
1245 : }
1246 : else
1247 2109749 : for (i=3;i<lr;i++)
1248 : {
1249 1982506 : ulong u = Fl_neg(uel(T,l-i+2), p);
1250 59522845 : for (k=3; k<i; k++)
1251 : {
1252 57540339 : ulong t = Fl_neg(uel(T,l-i+k), p);
1253 57540339 : u = Fl_addmul_pre(u, t, uel(r,k), p, pi);
1254 : }
1255 1982506 : r[i] = u;
1256 : }
1257 134238 : return Flx_renormalize(r,lr);
1258 : }
1259 :
1260 : /* Return new lgpol */
1261 : static long
1262 2129853 : Flx_lgrenormalizespec(GEN x, long lx)
1263 : {
1264 : long i;
1265 7434932 : for (i = lx-1; i>=0; i--)
1266 7434095 : if (x[i]) break;
1267 2129853 : return i+1;
1268 : }
1269 : /* allow pi = 0 */
1270 : static GEN
1271 23115 : Flx_invBarrett_Newton(GEN T, ulong p, ulong pi)
1272 : {
1273 23115 : long nold, lx, lz, lq, l = degpol(T), lQ;
1274 23115 : GEN q, y, z, x = zero_zv(l+1) + 2;
1275 23116 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1276 : pari_sp av;
1277 :
1278 23116 : y = T+2;
1279 23116 : q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1280 23115 : av = avma;
1281 : /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1282 :
1283 : /* initialize */
1284 23115 : x[0] = Fl_inv(q[0], p);
1285 23115 : if (lQ>1 && q[1])
1286 5105 : {
1287 5105 : ulong u = q[1];
1288 5105 : if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1289 5105 : x[1] = p - u; lx = 2;
1290 : }
1291 : else
1292 18010 : lx = 1;
1293 23115 : nold = 1;
1294 158676 : for (; mask > 1; set_avma(av))
1295 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1296 135598 : long i, lnew, nnew = nold << 1;
1297 :
1298 135598 : if (mask & 1) nnew--;
1299 135598 : mask >>= 1;
1300 :
1301 135598 : lnew = nnew + 1;
1302 135598 : lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1303 135601 : z = Flx_mulspec(x, q, p, pi, lx, lq); /* FIXME: high product */
1304 135561 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1305 135561 : z += 2;
1306 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1307 290632 : for (i = nold; i < lz; i++) if (z[i]) break;
1308 135561 : nold = nnew;
1309 135561 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1310 :
1311 : /* z + i represents (x*q - 1) / t^i */
1312 100752 : lz = Flx_lgrenormalizespec (z+i, lz-i);
1313 100751 : z = Flx_mulspec(x, z+i, p, pi, lx, lz); /* FIXME: low product */
1314 100754 : lz = lgpol(z); z += 2;
1315 100754 : if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1316 :
1317 100752 : lx = lz+ i;
1318 100752 : y = x + i; /* x -= z * t^i, in place */
1319 915379 : for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1320 : }
1321 23116 : x -= 2; setlg(x, lx + 2); x[1] = T[1];
1322 23116 : return x;
1323 : }
1324 :
1325 : /* allow pi = 0 */
1326 : static GEN
1327 158653 : Flx_invBarrett_pre(GEN T, ulong p, ulong pi)
1328 : {
1329 158653 : pari_sp ltop = avma;
1330 158653 : long l = lgpol(T);
1331 : GEN r;
1332 158653 : if (l < 3) return pol0_Flx(T[1]);
1333 157353 : if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1334 : {
1335 134238 : ulong c = T[l+1];
1336 134238 : if (c != 1)
1337 : {
1338 98118 : ulong ci = Fl_inv(c,p);
1339 98118 : T = Flx_Fl_mul_pre(T, ci, p, pi);
1340 98118 : r = Flx_invBarrett_basecase(T, p, pi);
1341 98118 : r = Flx_Fl_mul_pre(r, ci, p, pi);
1342 : }
1343 : else
1344 36120 : r = Flx_invBarrett_basecase(T, p, pi);
1345 : }
1346 : else
1347 23115 : r = Flx_invBarrett_Newton(T, p, pi);
1348 157354 : return gerepileuptoleaf(ltop, r);
1349 : }
1350 : GEN
1351 0 : Flx_invBarrett(GEN T, ulong p)
1352 0 : { return Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1353 :
1354 : /* allow pi = 0 */
1355 : GEN
1356 96871288 : Flx_get_red_pre(GEN T, ulong p, ulong pi)
1357 : {
1358 96871288 : if (typ(T)!=t_VECSMALL
1359 96834866 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1360 : Flx_BARRETT2_LIMIT))
1361 96854501 : return T;
1362 7613 : retmkvec2(Flx_invBarrett_pre(T, p, pi),T);
1363 : }
1364 : GEN
1365 14299788 : Flx_get_red(GEN T, ulong p)
1366 : {
1367 14299788 : if (typ(T)!=t_VECSMALL
1368 14299686 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1369 : Flx_BARRETT2_LIMIT))
1370 14294047 : return T;
1371 5194 : retmkvec2(Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)),T);
1372 : }
1373 :
1374 : /* separate from Flx_divrem for maximal speed. */
1375 : static GEN
1376 785631026 : Flx_rem_basecase(GEN x, GEN y, ulong p, ulong pi)
1377 : {
1378 : pari_sp av;
1379 : GEN z, c;
1380 : long dx,dy,dy1,dz,i,j;
1381 : ulong p1,inv;
1382 785631026 : long vs=x[1];
1383 :
1384 785631026 : dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1385 750193821 : dx = degpol(x);
1386 750089384 : dz = dx-dy; if (dz < 0) return Flx_copy(x);
1387 750089384 : x += 2; y += 2;
1388 750089384 : inv = y[dy];
1389 750089384 : if (inv != 1UL) inv = Fl_inv(inv,p);
1390 898987376 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1391 :
1392 751350283 : c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1393 749986951 : z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1394 :
1395 748314916 : if (!pi)
1396 : {
1397 479011259 : z[dz] = (inv*x[dx]) % p;
1398 1796081056 : for (i=dx-1; i>=dy; --i)
1399 : {
1400 1317069797 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1401 10439946025 : for (j=i-dy1; j<=i && j<=dz; j++)
1402 : {
1403 9122876228 : p1 += z[j]*y[i-j];
1404 9122876228 : if (p1 & HIGHBIT) p1 %= p;
1405 : }
1406 1317069797 : p1 %= p;
1407 1317069797 : z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1408 : }
1409 3268104068 : for (i=0; i<dy; i++)
1410 : {
1411 2789608984 : p1 = z[0]*y[i];
1412 14426735862 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1413 : {
1414 11637126878 : p1 += z[j]*y[i-j];
1415 11637126878 : if (p1 & HIGHBIT) p1 %= p;
1416 : }
1417 2789520521 : c[i] = Fl_sub(x[i], p1%p, p);
1418 : }
1419 : }
1420 : else
1421 : {
1422 269303657 : z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1423 833268278 : for (i=dx-1; i>=dy; --i)
1424 : {
1425 563632546 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1426 2362074341 : for (j=i-dy1; j<=i && j<=dz; j++)
1427 1797950224 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1428 564124117 : z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1429 : }
1430 2002153516 : for (i=0; i<dy; i++)
1431 : {
1432 1733340893 : p1 = Fl_mul_pre(z[0],y[i],p,pi);
1433 4685167604 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1434 2941962128 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1435 1719658150 : c[i] = Fl_sub(x[i], p1, p);
1436 : }
1437 : }
1438 914379158 : i = dy-1; while (i>=0 && !c[i]) i--;
1439 747307707 : set_avma(av); return Flx_renormalize(c-2, i+3);
1440 : }
1441 :
1442 : /* as FpX_divrem but working only on ulong types.
1443 : * if relevant, *pr is the last object on stack */
1444 : static GEN
1445 61795468 : Flx_divrem_basecase(GEN x, GEN y, ulong p, ulong pi, GEN *pr)
1446 : {
1447 : GEN z,q,c;
1448 : long dx,dy,dy1,dz,i,j;
1449 : ulong p1,inv;
1450 61795468 : long sv=x[1];
1451 :
1452 61795468 : dy = degpol(y);
1453 61793754 : if (dy<0) pari_err_INV("Flx_divrem",y);
1454 61793891 : if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p, pi);
1455 61793493 : if (!dy)
1456 : {
1457 7140466 : if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1458 7140427 : if (y[2] == 1UL) return Flx_copy(x);
1459 5135197 : return Flx_Fl_mul_pre(x, Fl_inv(y[2], p), p, pi);
1460 : }
1461 54653027 : dx = degpol(x);
1462 54655607 : dz = dx-dy;
1463 54655607 : if (dz < 0)
1464 : {
1465 1029382 : q = pol0_Flx(sv);
1466 1029373 : if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1467 1029372 : return q;
1468 : }
1469 53626225 : x += 2;
1470 53626225 : y += 2;
1471 53626225 : z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1472 53625037 : inv = uel(y, dy);
1473 53625037 : if (inv != 1UL) inv = Fl_inv(inv,p);
1474 78963694 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1475 :
1476 53627480 : if (SMALL_ULONG(p))
1477 : {
1478 51748096 : z[dz] = (inv*x[dx]) % p;
1479 131392101 : for (i=dx-1; i>=dy; --i)
1480 : {
1481 79644005 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1482 257565741 : for (j=i-dy1; j<=i && j<=dz; j++)
1483 : {
1484 177921736 : p1 += z[j]*y[i-j];
1485 177921736 : if (p1 & HIGHBIT) p1 %= p;
1486 : }
1487 79644005 : p1 %= p;
1488 79644005 : z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1489 : }
1490 : }
1491 : else
1492 : {
1493 1879384 : z[dz] = Fl_mul(inv, x[dx], p);
1494 9249040 : for (i=dx-1; i>=dy; --i)
1495 : { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1496 7370001 : p1 = p - uel(x,i);
1497 26363298 : for (j=i-dy1; j<=i && j<=dz; j++)
1498 18993297 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1499 7370001 : z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1500 : }
1501 : }
1502 53627135 : q = Flx_renormalize(z-2, dz+3);
1503 53626329 : if (!pr) return q;
1504 :
1505 26474735 : c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1506 26476577 : if (SMALL_ULONG(p))
1507 : {
1508 225445531 : for (i=0; i<dy; i++)
1509 : {
1510 200608113 : p1 = (ulong)z[0]*y[i];
1511 470333240 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1512 : {
1513 269725127 : p1 += (ulong)z[j]*y[i-j];
1514 269725127 : if (p1 & HIGHBIT) p1 %= p;
1515 : }
1516 200607767 : c[i] = Fl_sub(x[i], p1%p, p);
1517 : }
1518 : }
1519 : else
1520 : {
1521 16044506 : for (i=0; i<dy; i++)
1522 : {
1523 14405986 : p1 = Fl_mul(z[0],y[i],p);
1524 50246521 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1525 35840537 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1526 14405983 : c[i] = Fl_sub(x[i], p1, p);
1527 : }
1528 : }
1529 35603583 : i=dy-1; while (i>=0 && !c[i]) i--;
1530 26475938 : c = Flx_renormalize(c-2, i+3);
1531 26476721 : if (pr == ONLY_DIVIDES)
1532 427 : { if (lg(c) != 2) return NULL; }
1533 : else
1534 26476294 : *pr = c;
1535 26476581 : return q;
1536 : }
1537 :
1538 : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1539 : * and mg is the Barrett inverse of T. */
1540 : static GEN
1541 904146 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1542 : {
1543 : GEN q, r;
1544 904146 : long lt = degpol(T); /*We discard the leading term*/
1545 : long ld, lm, lT, lmg;
1546 904080 : ld = l-lt;
1547 904080 : lm = minss(ld, lgpol(mg));
1548 904406 : lT = Flx_lgrenormalizespec(T+2,lt);
1549 904537 : lmg = Flx_lgrenormalizespec(mg+2,lm);
1550 904460 : q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1551 903646 : q = Flx_mulspec(q+2,mg+2,p,pi,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1552 904579 : q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1553 903598 : if (!pr) return q;
1554 895892 : r = Flx_mulspec(q+2,T+2,p,pi,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1555 896863 : r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1556 896535 : if (pr == ONLY_REM) return r;
1557 427928 : *pr = r; return q;
1558 : }
1559 :
1560 : static GEN
1561 603718 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1562 : {
1563 603718 : GEN q = NULL, r = Flx_copy(x);
1564 603729 : long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1565 : long i;
1566 603730 : if (l <= lt)
1567 : {
1568 0 : if (pr == ONLY_REM) return Flx_copy(x);
1569 0 : if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1570 0 : if (pr) *pr = Flx_copy(x);
1571 0 : return pol0_Flx(v);
1572 : }
1573 603730 : if (lt <= 1)
1574 1300 : return Flx_divrem_basecase(x,T,p,pi,pr);
1575 602430 : if (pr != ONLY_REM && l>lm)
1576 28920 : { q = zero_zv(l-lt+1); q[1] = T[1]; }
1577 905784 : while (l>lm)
1578 : {
1579 303372 : GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,pi,&zr);
1580 303416 : long lz = lgpol(zr);
1581 303354 : if (pr != ONLY_REM)
1582 : {
1583 58086 : long lq = lgpol(zq);
1584 872756 : for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1585 : }
1586 4391895 : for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1587 303354 : l = l-lm+lz;
1588 : }
1589 602412 : if (pr == ONLY_REM)
1590 : {
1591 468665 : if (l > lt)
1592 468623 : r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi,ONLY_REM);
1593 : else
1594 42 : r = Flx_renormalize(r, l+2);
1595 468649 : r[1] = v; return r;
1596 : }
1597 133747 : if (l > lt)
1598 : {
1599 132190 : GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi, pr ? &r: NULL);
1600 132190 : if (!q) q = zq;
1601 : else
1602 : {
1603 27346 : long lq = lgpol(zq);
1604 158692 : for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1605 : }
1606 : }
1607 1557 : else if (pr)
1608 1535 : r = Flx_renormalize(r, l+2);
1609 133747 : q[1] = v; q = Flx_renormalize(q, lg(q));
1610 133764 : if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1611 133764 : if (pr) { r[1] = v; *pr = r; }
1612 133764 : return q;
1613 : }
1614 :
1615 : /* allow pi = 0 (SMALL_ULONG) */
1616 : GEN
1617 79243306 : Flx_divrem_pre(GEN x, GEN T, ulong p, ulong pi, GEN *pr)
1618 : {
1619 : GEN B, y;
1620 : long dy, dx, d;
1621 79243306 : if (pr==ONLY_REM) return Flx_rem_pre(x, T, p, pi);
1622 61920238 : y = get_Flx_red(T, &B);
1623 61930101 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1624 61926958 : if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1625 61793107 : return Flx_divrem_basecase(x,y,p,pi,pr);
1626 : else
1627 : {
1628 134666 : pari_sp av = avma;
1629 134666 : GEN mg = B? B: Flx_invBarrett_pre(y, p, pi);
1630 134666 : GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pi,pr);
1631 134666 : if (!q1) return gc_NULL(av);
1632 134666 : if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
1633 126366 : return gc_all(av, 2, &q1, pr);
1634 : }
1635 : }
1636 : GEN
1637 30302312 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1638 30302312 : { return Flx_divrem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p), pr); }
1639 :
1640 : GEN
1641 909210950 : Flx_rem_pre(GEN x, GEN T, ulong p, ulong pi)
1642 : {
1643 909210950 : GEN B, y = get_Flx_red(T, &B);
1644 909038418 : long d = degpol(x) - degpol(y);
1645 908655651 : if (d < 0) return Flx_copy(x);
1646 786128680 : if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1647 785540224 : return Flx_rem_basecase(x,y,p, pi);
1648 : else
1649 : {
1650 469053 : pari_sp av=avma;
1651 469053 : GEN mg = B ? B: Flx_invBarrett_pre(y, p, pi);
1652 469053 : GEN r = Flx_divrem_Barrett(x, mg, y, p, pi, ONLY_REM);
1653 469047 : return gerepileuptoleaf(av, r);
1654 : }
1655 : }
1656 : GEN
1657 41841936 : Flx_rem(GEN x, GEN T, ulong p)
1658 41841936 : { return Flx_rem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1659 :
1660 : /* reduce T mod (X^n - 1, p). Shallow function */
1661 : GEN
1662 5106845 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1663 : {
1664 5106845 : long i, j, L = lg(T), l = n+2;
1665 : GEN S;
1666 5106845 : if (L <= l || n & ~LGBITS) return T;
1667 3450 : S = cgetg(l, t_VECSMALL);
1668 3450 : S[1] = T[1];
1669 14013 : for (i = 2; i < l; i++) S[i] = T[i];
1670 9420 : for (j = 2; i < L; i++) {
1671 5970 : S[j] = Fl_add(S[j], T[i], p);
1672 5970 : if (++j == l) j = 2;
1673 : }
1674 3450 : return Flx_renormalize(S, l);
1675 : }
1676 : /* reduce T mod (X^n + 1, p). Shallow function */
1677 : GEN
1678 30321 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
1679 : {
1680 30321 : long i, j, L = lg(T), l = n+2;
1681 : GEN S;
1682 30321 : if (L <= l || n & ~LGBITS) return T;
1683 2682 : S = cgetg(l, t_VECSMALL);
1684 2682 : S[1] = T[1];
1685 11347 : for (i = 2; i < l; i++) S[i] = T[i];
1686 6974 : for (j = 2; i < L; i++) {
1687 4292 : S[j] = Fl_sub(S[j], T[i], p);
1688 4292 : if (++j == l) j = 2;
1689 : }
1690 2682 : return Flx_renormalize(S, l);
1691 : }
1692 :
1693 : struct _Flxq {
1694 : GEN aut, T;
1695 : ulong p, pi;
1696 : };
1697 : /* allow pi = 0 */
1698 : static void
1699 69461215 : set_Flxq_pre(struct _Flxq *D, GEN T, ulong p, ulong pi)
1700 : {
1701 69461215 : D->p = p;
1702 69461215 : D->pi = pi;
1703 69461215 : D->T = Flx_get_red_pre(T, p, pi);
1704 69457233 : }
1705 : static void
1706 69005 : set_Flxq(struct _Flxq *D, GEN T, ulong p)
1707 69005 : { set_Flxq_pre(D, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1708 :
1709 : static GEN
1710 0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1711 : {
1712 0 : struct _Flxq *D = (struct _Flxq*) E;
1713 0 : return Flx_divrem_pre(x, y, D->p, D->pi, r);
1714 : }
1715 : static GEN
1716 389834 : _Flx_add(void * E, GEN x, GEN y) {
1717 389834 : struct _Flxq *D = (struct _Flxq*) E;
1718 389834 : return Flx_add(x, y, D->p);
1719 : }
1720 : static GEN
1721 10482861 : _Flx_mul(void *E, GEN x, GEN y) {
1722 10482861 : struct _Flxq *D = (struct _Flxq*) E;
1723 10482861 : return Flx_mul_pre(x, y, D->p, D->pi);
1724 : }
1725 : static GEN
1726 0 : _Flx_sqr(void *E, GEN x) {
1727 0 : struct _Flxq *D = (struct _Flxq*) E;
1728 0 : return Flx_sqr_pre(x, D->p, D->pi);
1729 : }
1730 :
1731 : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1732 :
1733 : GEN
1734 0 : Flx_digits(GEN x, GEN T, ulong p)
1735 : {
1736 : struct _Flxq D;
1737 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
1738 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1739 0 : return gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1740 : }
1741 :
1742 : GEN
1743 0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1744 : {
1745 : struct _Flxq D;
1746 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1747 0 : return gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1748 : }
1749 :
1750 : long
1751 4162412 : Flx_val(GEN x)
1752 : {
1753 4162412 : long i, l=lg(x);
1754 4162412 : if (l==2) return LONG_MAX;
1755 4171342 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1756 4162412 : return i-2;
1757 : }
1758 : long
1759 26292982 : Flx_valrem(GEN x, GEN *Z)
1760 : {
1761 26292982 : long v, i, l=lg(x);
1762 : GEN y;
1763 26292982 : if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1764 28475758 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1765 26292982 : v = i-2;
1766 26292982 : if (v == 0) { *Z = x; return 0; }
1767 1031784 : l -= v;
1768 1031784 : y = cgetg(l, t_VECSMALL); y[1] = x[1];
1769 2637923 : for (i=2; i<l; i++) y[i] = x[i+v];
1770 1029230 : *Z = y; return v;
1771 : }
1772 :
1773 : GEN
1774 21171973 : Flx_deriv(GEN z, ulong p)
1775 : {
1776 21171973 : long i,l = lg(z)-1;
1777 : GEN x;
1778 21171973 : if (l < 2) l = 2;
1779 21171973 : x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1780 21170440 : if (HIGHWORD(l | p))
1781 57463544 : for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1782 : else
1783 85481995 : for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1784 21170562 : return Flx_renormalize(x,l);
1785 : }
1786 :
1787 : static GEN
1788 422849 : Flx_integXn(GEN x, long n, ulong p)
1789 : {
1790 422849 : long i, lx = lg(x);
1791 : GEN y;
1792 422849 : if (lx == 2) return Flx_copy(x);
1793 413036 : y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1794 2097086 : for (i=2; i<lx; i++)
1795 : {
1796 1683469 : ulong xi = uel(x,i);
1797 1683469 : if (xi == 0)
1798 13345 : uel(y,i) = 0;
1799 : else
1800 : {
1801 1670124 : ulong j = n+i-1;
1802 1670124 : ulong d = ugcd(j, xi);
1803 1670050 : if (d==1)
1804 1018469 : uel(y,i) = Fl_div(xi, j, p);
1805 : else
1806 651581 : uel(y,i) = Fl_div(xi/d, j/d, p);
1807 : }
1808 : }
1809 413617 : return Flx_renormalize(y, lx);;
1810 : }
1811 :
1812 : GEN
1813 0 : Flx_integ(GEN x, ulong p)
1814 : {
1815 0 : long i, lx = lg(x);
1816 : GEN y;
1817 0 : if (lx == 2) return Flx_copy(x);
1818 0 : y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1819 0 : uel(y,2) = 0;
1820 0 : for (i=3; i<=lx; i++)
1821 0 : uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1822 0 : return Flx_renormalize(y, lx+1);;
1823 : }
1824 :
1825 : /* assume p prime */
1826 : GEN
1827 13482 : Flx_diff1(GEN P, ulong p)
1828 : {
1829 13482 : return Flx_sub(Flx_translate1(P, p), P, p);
1830 : }
1831 :
1832 : GEN
1833 420541 : Flx_deflate(GEN x0, long d)
1834 : {
1835 : GEN z, y, x;
1836 420541 : long i,id, dy, dx = degpol(x0);
1837 420541 : if (d == 1 || dx <= 0) return Flx_copy(x0);
1838 357052 : dy = dx/d;
1839 357052 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1840 357051 : z = y + 2;
1841 357051 : x = x0+ 2;
1842 1160930 : for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1843 357051 : return y;
1844 : }
1845 :
1846 : GEN
1847 160035 : Flx_inflate(GEN x0, long d)
1848 : {
1849 160035 : long i, id, dy, dx = degpol(x0);
1850 160033 : GEN x = x0 + 2, z, y;
1851 160033 : if (dx <= 0) return Flx_copy(x0);
1852 158972 : dy = dx*d;
1853 158972 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1854 158962 : z = y + 2;
1855 8749755 : for (i=0; i<=dy; i++) z[i] = 0;
1856 4255151 : for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1857 158962 : return y;
1858 : }
1859 :
1860 : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1861 : GEN
1862 147099 : Flx_splitting(GEN p, long k)
1863 : {
1864 147099 : long n = degpol(p), v = p[1], m, i, j, l;
1865 : GEN r;
1866 :
1867 147099 : m = n/k;
1868 147099 : r = cgetg(k+1,t_VEC);
1869 678611 : for(i=1; i<=k; i++)
1870 : {
1871 531516 : gel(r,i) = cgetg(m+3, t_VECSMALL);
1872 531512 : mael(r,i,1) = v;
1873 : }
1874 4429832 : for (j=1, i=0, l=2; i<=n; i++)
1875 : {
1876 4282737 : mael(r,j,l) = p[2+i];
1877 4282737 : if (j==k) { j=1; l++; } else j++;
1878 : }
1879 678614 : for(i=1; i<=k; i++)
1880 531522 : gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1881 147092 : return r;
1882 : }
1883 :
1884 : /* ux + vy */
1885 : static GEN
1886 416536 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p, ulong pi)
1887 416536 : { return Flx_add(Flx_mul_pre(u,x, p,pi), Flx_mul_pre(v,y, p,pi), p); }
1888 :
1889 : static GEN
1890 24679 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p, ulong pi)
1891 : {
1892 24679 : GEN res = cgetg(3, t_COL);
1893 24679 : gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p, pi);
1894 24680 : gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p, pi);
1895 24679 : return res;
1896 : }
1897 :
1898 : #if 0
1899 : static GEN
1900 : FlxM_mul2_old(GEN M, GEN N, ulong p)
1901 : {
1902 : GEN res = cgetg(3, t_MAT);
1903 : gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1904 : gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1905 : return res;
1906 : }
1907 : #endif
1908 : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1909 : static GEN
1910 6481 : FlxM_mul2(GEN A, GEN B, ulong p, ulong pi)
1911 : {
1912 6481 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1913 6481 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1914 6481 : GEN M1 = Flx_mul_pre(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p, pi);
1915 6481 : GEN M2 = Flx_mul_pre(Flx_add(A21,A22, p), B11, p, pi);
1916 6481 : GEN M3 = Flx_mul_pre(A11, Flx_sub(B12,B22, p), p, pi);
1917 6481 : GEN M4 = Flx_mul_pre(A22, Flx_sub(B21,B11, p), p, pi);
1918 6481 : GEN M5 = Flx_mul_pre(Flx_add(A11,A12, p), B22, p, pi);
1919 6481 : GEN M6 = Flx_mul_pre(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p, pi);
1920 6481 : GEN M7 = Flx_mul_pre(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p, pi);
1921 6481 : GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1922 6481 : GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1923 6481 : retmkmat22(Flx_add(T1,T2, p), Flx_add(M3,M5, p),
1924 : Flx_add(M2,M4, p), Flx_add(T3,T4, p));
1925 : }
1926 :
1927 : /* Return [0,1;1,-q]*M */
1928 : static GEN
1929 6309 : Flx_FlxM_qmul(GEN q, GEN M, ulong p, ulong pi)
1930 : {
1931 6309 : GEN u = Flx_mul_pre(gcoeff(M,2,1), q, p, pi);
1932 6309 : GEN v = Flx_mul_pre(gcoeff(M,2,2), q, p, pi);
1933 6309 : retmkmat22(gcoeff(M,2,1), gcoeff(M,2,2),
1934 : Flx_sub(gcoeff(M,1,1), u, p), Flx_sub(gcoeff(M,1,2), v, p));
1935 : }
1936 :
1937 : static GEN
1938 895 : matid2_FlxM(long v)
1939 895 : { retmkmat22(pol1_Flx(v),pol0_Flx(v),pol0_Flx(v),pol1_Flx(v)); }
1940 :
1941 : static GEN
1942 13 : matJ2_FlxM(long v)
1943 13 : { retmkmat22(pol0_Flx(v),pol1_Flx(v),pol1_Flx(v),pol0_Flx(v)); }
1944 :
1945 : struct Flx_res
1946 : {
1947 : ulong res, lc;
1948 : long deg0, deg1, off;
1949 : };
1950 :
1951 : INLINE void
1952 9405 : Flx_halfres_update_pre(long da, long db, long dr, ulong p, ulong pi, struct Flx_res *res)
1953 : {
1954 9405 : if (dr >= 0)
1955 : {
1956 9405 : if (res->lc != 1)
1957 : {
1958 7596 : if (pi)
1959 : {
1960 3127 : res->lc = Fl_powu_pre(res->lc, da - dr, p, pi);
1961 3127 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1962 : } else
1963 : {
1964 4469 : res->lc = Fl_powu(res->lc, da - dr, p);
1965 4469 : res->res = Fl_mul(res->res, res->lc, p);
1966 : }
1967 : }
1968 9405 : if (both_odd(da + res->off, db + res->off))
1969 63 : res->res = Fl_neg(res->res, p);
1970 : } else
1971 : {
1972 0 : if (db == 0)
1973 : {
1974 0 : if (res->lc != 1)
1975 : {
1976 0 : if (pi)
1977 : {
1978 0 : res->lc = Fl_powu_pre(res->lc, da, p, pi);
1979 0 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1980 : } else
1981 : {
1982 0 : res->lc = Fl_powu(res->lc, da, p);
1983 0 : res->res = Fl_mul(res->res, res->lc, p);
1984 : }
1985 : }
1986 : } else
1987 0 : res->res = 0;
1988 : }
1989 9405 : }
1990 :
1991 : static GEN
1992 1108953 : Flx_halfres_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *pa, GEN *pb, struct Flx_res *res)
1993 : {
1994 1108953 : pari_sp av = avma;
1995 : GEN u, u1, v, v1, M;
1996 1108953 : long vx = a[1], n = lgpol(a)>>1;
1997 1108949 : u1 = v = pol0_Flx(vx);
1998 1108944 : u = v1 = pol1_Flx(vx);
1999 6844671 : while (lgpol(b)>n)
2000 : {
2001 : GEN r, q;
2002 5735717 : q = Flx_divrem_pre(a,b,p,pi, &r);
2003 5735780 : if (res)
2004 : {
2005 8362 : long da = degpol(a), db=degpol(b), dr = degpol(r);
2006 8362 : res->lc = b[db+2];
2007 8362 : if (dr >= n)
2008 7133 : Flx_halfres_update_pre(da, db, dr, p, pi, res);
2009 : else
2010 : {
2011 1229 : res->deg0 = da;
2012 1229 : res->deg1 = db;
2013 : }
2014 : }
2015 5735780 : a = b; b = r; swap(u,u1); swap(v,v1);
2016 5735780 : u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
2017 5735680 : v1 = Flx_sub(v1, Flx_mul(v, q, p), p);
2018 5735730 : if (gc_needed(av,2))
2019 : {
2020 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
2021 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
2022 : }
2023 : }
2024 1108804 : M = mkmat22(u,v,u1,v1); *pa = a; *pb = b;
2025 1108935 : return gc_all(av,3, &M, pa, pb);
2026 : }
2027 :
2028 : static GEN Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res);
2029 :
2030 : static GEN
2031 19245 : Flx_halfres_split(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2032 : {
2033 19245 : pari_sp av = avma;
2034 : GEN R, S, T, V1, V2;
2035 : GEN x1, y1, r, q;
2036 19245 : long l = lgpol(x), n = l>>1, k;
2037 19245 : if (lgpol(y) <= n)
2038 855 : { *a = Flx_copy(x); *b = Flx_copy(y); return matid2_FlxM(x[1]); }
2039 18390 : if (res)
2040 : {
2041 3262 : res->lc = Flx_lead(y);
2042 3262 : res->deg0 -= n;
2043 3262 : res->deg1 -= n;
2044 3262 : res->off += n;
2045 : }
2046 18390 : R = Flx_halfres_i(Flx_shift(x,-n),Flx_shift(y,-n),p,pi,a,b,res);
2047 18390 : if (res)
2048 : {
2049 3262 : res->off -= n;
2050 3262 : res->deg0 += n;
2051 3262 : res->deg1 += n;
2052 : }
2053 18390 : V1 = FlxM_Flx_mul2(R, Flxn_red(x,n), Flxn_red(y,n), p, pi);
2054 18391 : x1 = Flx_add(Flx_shift(*a,n), gel(V1,1), p);
2055 18390 : y1 = Flx_add(Flx_shift(*b,n), gel(V1,2), p);
2056 18390 : if (lgpol(y1) <= n)
2057 12101 : { *a = x1; *b = y1; return gc_all(av, 3, &R, a, b); }
2058 6289 : k = 2*n-degpol(y1);
2059 6289 : q = Flx_divrem_pre(x1, y1, p, pi, &r);
2060 6289 : if (res)
2061 : {
2062 1043 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
2063 1043 : if (dy1 < degpol(y))
2064 185 : Flx_halfres_update_pre(res->deg0, res->deg1, dy1, p, pi, res);
2065 1043 : res->lc = uel(y1, dy1+2);
2066 1043 : res->deg0 = dx1;
2067 1043 : res->deg1 = dy1;
2068 1043 : if (dr >= n)
2069 : {
2070 1043 : Flx_halfres_update_pre(dx1, dy1, dr, p, pi, res);
2071 1043 : res->deg0 = dy1;
2072 1043 : res->deg1 = dr;
2073 : }
2074 1043 : res->deg0 -= k;
2075 1043 : res->deg1 -= k;
2076 1043 : res->off += k;
2077 : }
2078 6289 : S = Flx_halfres_i(Flx_shift(y1,-k), Flx_shift(r,-k), p, pi, a, b, res);
2079 6289 : if (res)
2080 : {
2081 1043 : res->deg0 += k;
2082 1043 : res->deg1 += k;
2083 1043 : res->off -= k;
2084 : }
2085 6289 : T = FlxM_mul2(S, Flx_FlxM_qmul(q, R, p,pi), p, pi);
2086 6289 : V2 = FlxM_Flx_mul2(S, Flxn_red(y1,k), Flxn_red(r,k), p, pi);
2087 6288 : *a = Flx_add(Flx_shift(*a,k), gel(V2,1), p);
2088 6289 : *b = Flx_add(Flx_shift(*b,k), gel(V2,2), p);
2089 6289 : return gc_all(av, 3, &T, a, b);
2090 : }
2091 :
2092 : static GEN
2093 1128199 : Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, struct Flx_res *res)
2094 : {
2095 1128199 : if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
2096 1108953 : return Flx_halfres_basecase(x, y, p, pi, a, b, res);
2097 19245 : return Flx_halfres_split(x, y, p, pi, a, b, res);
2098 : }
2099 :
2100 : static GEN
2101 1102475 : Flx_halfgcd_all_i(GEN x, GEN y, ulong p, ulong pi, GEN *pa, GEN *pb)
2102 : {
2103 : GEN a, b, R;
2104 1102475 : R = Flx_halfres_i(x, y, p, pi, &a, &b, NULL);
2105 1102490 : if (pa) *pa = a;
2106 1102490 : if (pb) *pb = b;
2107 1102490 : return R;
2108 : }
2109 :
2110 : /* Return M in GL_2(Fl[X]) such that:
2111 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
2112 : */
2113 :
2114 : GEN
2115 1102484 : Flx_halfgcd_all_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b)
2116 : {
2117 : pari_sp av;
2118 : GEN R, q ,r;
2119 1102484 : long lx = lgpol(x), ly = lgpol(y);
2120 1102476 : if (!lx)
2121 : {
2122 0 : if (a) *a = Flx_copy(y);
2123 0 : if (b) *b = Flx_copy(x);
2124 0 : return matJ2_FlxM(x[1]);
2125 : }
2126 1102476 : if (ly < lx) return Flx_halfgcd_all_i(x, y, p, pi, a, b);
2127 8209 : av = avma;
2128 8209 : q = Flx_divrem(y,x,p,&r);
2129 8209 : R = Flx_halfgcd_all_i(x, r, p, pi, a, b);
2130 8209 : gcoeff(R,1,1) = Flx_sub(gcoeff(R,1,1), Flx_mul_pre(q,gcoeff(R,1,2), p,pi), p);
2131 8209 : gcoeff(R,2,1) = Flx_sub(gcoeff(R,2,1), Flx_mul_pre(q,gcoeff(R,2,2), p,pi), p);
2132 8209 : return !a && b ? gc_all(av, 2, &R, b): gc_all(av, 1+!!a+!!b, &R, a, b);
2133 : }
2134 :
2135 : GEN
2136 154 : Flx_halfgcd_all(GEN x, GEN y, ulong p, GEN *a, GEN *b)
2137 154 : { return Flx_halfgcd_all_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), a, b); }
2138 :
2139 : GEN
2140 846890 : Flx_halfgcd_pre(GEN x, GEN y, ulong p, ulong pi)
2141 846890 : { return Flx_halfgcd_all_pre(x, y, p, pi, NULL, NULL); }
2142 :
2143 : GEN
2144 0 : Flx_halfgcd(GEN x, GEN y, ulong p)
2145 0 : { return Flx_halfgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2146 :
2147 : /*Do not garbage collect*/
2148 : static GEN
2149 82990050 : Flx_gcd_basecase(GEN a, GEN b, ulong p, ulong pi)
2150 : {
2151 82990050 : pari_sp av = avma;
2152 82990050 : ulong iter = 0;
2153 82990050 : if (lg(b) > lg(a)) swap(a, b);
2154 286492249 : while (lgpol(b))
2155 : {
2156 203091885 : GEN c = Flx_rem_pre(a,b,p,pi);
2157 203502199 : iter++; a = b; b = c;
2158 203502199 : if (gc_needed(av,2))
2159 : {
2160 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
2161 0 : gerepileall(av,2, &a,&b);
2162 : }
2163 : }
2164 82939094 : return iter < 2 ? Flx_copy(a) : a;
2165 : }
2166 :
2167 : GEN
2168 84641771 : Flx_gcd_pre(GEN x, GEN y, ulong p, ulong pi)
2169 : {
2170 84641771 : pari_sp av = avma;
2171 : long lim;
2172 84641771 : if (!lgpol(x)) return Flx_copy(y);
2173 83000019 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2174 83004146 : while (lgpol(y) >= lim)
2175 : {
2176 144 : if (lgpol(y)<=(lgpol(x)>>1))
2177 : {
2178 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2179 0 : x = y; y = r;
2180 : }
2181 144 : (void) Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2182 144 : if (gc_needed(av,2))
2183 : {
2184 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
2185 0 : gerepileall(av,2,&x,&y);
2186 : }
2187 : }
2188 82987981 : return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p,pi));
2189 : }
2190 : GEN
2191 32455094 : Flx_gcd(GEN x, GEN y, ulong p)
2192 32455094 : { return Flx_gcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2193 :
2194 : int
2195 8540411 : Flx_is_squarefree(GEN z, ulong p)
2196 : {
2197 8540411 : pari_sp av = avma;
2198 8540411 : GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
2199 8540309 : return gc_bool(av, degpol(d) == 0);
2200 : }
2201 :
2202 : static long
2203 126944 : Flx_is_smooth_squarefree(GEN f, long r, ulong p, ulong pi)
2204 : {
2205 126944 : pari_sp av = avma;
2206 : long i;
2207 126944 : GEN sx = polx_Flx(f[1]), a = sx;
2208 535115 : for(i=1;;i++)
2209 : {
2210 535115 : if (degpol(f)<=r) return gc_long(av,1);
2211 513229 : a = Flxq_powu_pre(Flx_rem_pre(a,f,p,pi), p, f, p, pi);
2212 513052 : if (Flx_equal(a, sx)) return gc_long(av,1);
2213 509488 : if (i==r) return gc_long(av,0);
2214 407740 : f = Flx_div_pre(f, Flx_gcd_pre(Flx_sub(a,sx,p),f,p,pi),p,pi);
2215 : }
2216 : }
2217 :
2218 : static long
2219 8199 : Flx_is_l_pow(GEN x, ulong p)
2220 : {
2221 8199 : ulong i, lx = lgpol(x);
2222 16380 : for (i=1; i<lx; i++)
2223 14697 : if (x[i+2] && i%p) return 0;
2224 1683 : return 1;
2225 : }
2226 :
2227 : int
2228 126931 : Flx_is_smooth_pre(GEN g, long r, ulong p, ulong pi)
2229 : {
2230 : while (1)
2231 8197 : {
2232 126931 : GEN f = Flx_gcd_pre(g, Flx_deriv(g, p), p, pi);
2233 126758 : if (!Flx_is_smooth_squarefree(Flx_div_pre(g, f, p, pi), r, p, pi))
2234 101751 : return 0;
2235 25217 : if (degpol(f)==0) return 1;
2236 8185 : g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2237 : }
2238 : }
2239 : int
2240 74256 : Flx_is_smooth(GEN g, long r, ulong p)
2241 74256 : { return Flx_is_smooth_pre(g, r, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2242 :
2243 : static GEN
2244 6266697 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2245 : {
2246 6266697 : pari_sp av=avma;
2247 : GEN u,v,u1,v1;
2248 6266697 : long vx = a[1];
2249 6266697 : v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2250 6266536 : if (ptu) { u = pol1_Flx(vx); u1 = pol0_Flx(vx); }
2251 28014991 : while (lgpol(b))
2252 : {
2253 21747411 : GEN r, q = Flx_divrem_pre(a,b,p,pi, &r);
2254 21748703 : a = b; b = r;
2255 21748703 : if (ptu)
2256 : {
2257 2435900 : swap(u,u1);
2258 2435900 : u1 = Flx_sub(u1, Flx_mul_pre(u, q, p, pi), p);
2259 : }
2260 21748719 : swap(v,v1);
2261 21748719 : v1 = Flx_sub(v1, Flx_mul_pre(v, q, p, pi), p);
2262 21748457 : if (gc_needed(av,2))
2263 : {
2264 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(a));
2265 0 : gerepileall(av,ptu ? 6: 4, &a,&b,&v,&v1,&u,&u1);
2266 : }
2267 : }
2268 6266577 : if (ptu) *ptu = u;
2269 6266577 : *ptv = v;
2270 6266577 : return a;
2271 : }
2272 :
2273 : static GEN
2274 147433 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2275 : {
2276 : GEN u, v;
2277 147433 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2278 147433 : GEN V = cgetg(expu(lgpol(y))+2,t_VEC);
2279 147433 : long i, n = 0, vs = x[1];
2280 401414 : while (lgpol(y) >= lim)
2281 : {
2282 253981 : if (lgpol(y)<=(lgpol(x)>>1))
2283 : {
2284 26 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2285 26 : x = y; y = r;
2286 26 : gel(V,++n) = mkmat22(pol0_Flx(vs),pol1_Flx(vs),pol1_Flx(vs),Flx_neg(q,p));
2287 : } else
2288 253955 : gel(V,++n) = Flx_halfgcd_all_pre(x, y, p, pi, &x, &y);
2289 : }
2290 147432 : y = Flx_extgcd_basecase(x,y,p,pi,&u,&v);
2291 253981 : for (i = n; i>1; i--)
2292 : {
2293 106548 : GEN R = gel(V,i);
2294 106548 : GEN u1 = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2295 106548 : GEN v1 = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2296 106548 : u = u1; v = v1;
2297 : }
2298 : {
2299 147433 : GEN R = gel(V,1);
2300 147433 : if (ptu)
2301 6543 : *ptu = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2302 147433 : *ptv = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2303 : }
2304 147433 : return y;
2305 : }
2306 :
2307 : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2308 : * ux + vy = gcd (mod p) */
2309 : GEN
2310 6266685 : Flx_extgcd_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2311 : {
2312 6266685 : pari_sp av = avma;
2313 : GEN d;
2314 6266685 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2315 6266694 : if (lgpol(y) >= lim)
2316 147433 : d = Flx_extgcd_halfgcd(x, y, p, pi, ptu, ptv);
2317 : else
2318 6119248 : d = Flx_extgcd_basecase(x, y, p, pi, ptu, ptv);
2319 6266580 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
2320 : }
2321 : GEN
2322 855131 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2323 855131 : { return Flx_extgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptu, ptv); }
2324 :
2325 : static GEN
2326 1043 : Flx_halfres_pre(GEN x, GEN y, ulong p, ulong pi, GEN *a, GEN *b, ulong *r)
2327 : {
2328 : struct Flx_res res;
2329 : GEN R;
2330 : long dB;
2331 :
2332 1043 : res.res = *r;
2333 1043 : res.lc = Flx_lead(y);
2334 1043 : res.deg0 = degpol(x);
2335 1043 : res.deg1 = degpol(y);
2336 1043 : res.off = 0;
2337 1043 : R = Flx_halfres_i(x, y, p, pi, a, b, &res);
2338 1044 : dB = degpol(*b);
2339 1044 : if (dB < degpol(y))
2340 1044 : Flx_halfres_update_pre(res.deg0, res.deg1, dB, p, pi, &res);
2341 1044 : *r = res.res;
2342 1044 : return R;
2343 : }
2344 :
2345 : static ulong
2346 10269330 : Flx_resultant_basecase_pre(GEN a, GEN b, ulong p, ulong pi)
2347 : {
2348 : pari_sp av;
2349 : long da,db,dc;
2350 10269330 : ulong lb, res = 1UL;
2351 : GEN c;
2352 :
2353 10269330 : da = degpol(a);
2354 10269195 : db = degpol(b);
2355 10269299 : if (db > da)
2356 : {
2357 0 : swapspec(a,b, da,db);
2358 0 : if (both_odd(da,db)) res = p-res;
2359 : }
2360 10269299 : else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2361 10269299 : av = avma;
2362 107321046 : while (db)
2363 : {
2364 97071842 : lb = b[db+2];
2365 97071842 : c = Flx_rem_pre(a,b, p,pi);
2366 96712612 : a = b; b = c; dc = degpol(c);
2367 96684230 : if (dc < 0) return gc_long(av,0);
2368 :
2369 96678752 : if (both_odd(da,db)) res = p - res;
2370 96681214 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, da - dc, p, pi), p);
2371 97064698 : if (gc_needed(av,2))
2372 : {
2373 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant (da = %ld)",da);
2374 0 : gerepileall(av,2, &a,&b);
2375 : }
2376 97051747 : da = db; /* = degpol(a) */
2377 97051747 : db = dc; /* = degpol(b) */
2378 : }
2379 10249204 : return gc_ulong(av, Fl_mul(res, Fl_powu_pre(b[2], da, p, pi), p));
2380 : }
2381 :
2382 : ulong
2383 10271290 : Flx_resultant_pre(GEN x, GEN y, ulong p, ulong pi)
2384 : {
2385 10271290 : pari_sp av = avma;
2386 : long lim;
2387 10271290 : ulong res = 1;
2388 10271290 : long dx = degpol(x), dy = degpol(y);
2389 10270860 : if (dx < 0 || dy < 0) return 0;
2390 10269418 : if (dx < dy)
2391 : {
2392 1065302 : swap(x,y);
2393 1065302 : if (both_odd(dx, dy))
2394 1906 : res = Fl_neg(res, p);
2395 : }
2396 10269414 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2397 10270284 : while (lgpol(y) >= lim)
2398 : {
2399 851 : if (lgpol(y)<=(lgpol(x)>>1))
2400 : {
2401 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2402 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2403 0 : ulong ly = y[dy+2];
2404 0 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2405 0 : if (both_odd(dx, dy))
2406 0 : res = Fl_neg(res, p);
2407 0 : x = y; y = r;
2408 : }
2409 851 : (void) Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2410 852 : if (gc_needed(av,2))
2411 : {
2412 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_res (y = %ld)",degpol(y));
2413 0 : gerepileall(av,2,&x,&y);
2414 : }
2415 : }
2416 10269349 : return gc_ulong(av, Fl_mul(res, Flx_resultant_basecase_pre(x, y, p, pi), p));
2417 : }
2418 :
2419 : ulong
2420 4732378 : Flx_resultant(GEN a, GEN b, ulong p)
2421 4732378 : { return Flx_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2422 :
2423 : /* If resultant is 0, *ptU and *ptV are not set */
2424 : static ulong
2425 53 : Flx_extresultant_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2426 : {
2427 53 : GEN z,q,u,v, x = a, y = b;
2428 53 : ulong lb, res = 1UL;
2429 53 : pari_sp av = avma;
2430 : long dx, dy, dz;
2431 53 : long vs = a[1];
2432 :
2433 53 : u = pol0_Flx(vs);
2434 53 : v = pol1_Flx(vs); /* v = 1 */
2435 53 : dx = degpol(x);
2436 53 : dy = degpol(y);
2437 764 : while (dy)
2438 : { /* b u = x (a), b v = y (a) */
2439 711 : lb = y[dy+2];
2440 711 : q = Flx_divrem_pre(x,y, p, pi, &z);
2441 711 : x = y; y = z; /* (x,y) = (y, x - q y) */
2442 711 : dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2443 711 : z = Flx_sub(u, Flx_mul_pre(q,v, p, pi), p);
2444 711 : u = v; v = z; /* (u,v) = (v, u - q v) */
2445 :
2446 711 : if (both_odd(dx,dy)) res = p - res;
2447 711 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, dx-dz, p, pi), p);
2448 711 : dx = dy; /* = degpol(x) */
2449 711 : dy = dz; /* = degpol(y) */
2450 : }
2451 53 : res = Fl_mul(res, Fl_powu_pre(y[2], dx, p, pi), p);
2452 53 : lb = Fl_mul(res, Fl_inv(y[2],p), p);
2453 53 : v = gerepileuptoleaf(av, Flx_Fl_mul_pre(v, lb, p, pi));
2454 53 : av = avma;
2455 53 : u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul_pre(b,v,p,pi), p);
2456 53 : u = gerepileuptoleaf(av, Flx_div_pre(u,a,p,pi)); /* = (res - b v) / a */
2457 53 : *ptU = u;
2458 53 : *ptV = v; return res;
2459 : }
2460 :
2461 : ulong
2462 53 : Flx_extresultant_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2463 : {
2464 53 : pari_sp av=avma;
2465 : GEN u, v, R;
2466 53 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2467 53 : ulong res = 1, res1;
2468 53 : long dx = degpol(x), dy = degpol(y);
2469 53 : if (dy > dx)
2470 : {
2471 13 : swap(x,y); lswap(dx,dy);
2472 13 : if (both_odd(dx,dy)) res = p-res;
2473 13 : R = matJ2_FlxM(x[1]);
2474 40 : } else R = matid2_FlxM(x[1]);
2475 53 : if (dy < 0) return 0;
2476 245 : while (lgpol(y) >= lim)
2477 : {
2478 : GEN M;
2479 192 : if (lgpol(y)<=(lgpol(x)>>1))
2480 : {
2481 20 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2482 20 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2483 20 : ulong ly = y[dy+2];
2484 20 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2485 20 : if (both_odd(dx, dy))
2486 0 : res = Fl_neg(res, p);
2487 20 : x = y; y = r;
2488 20 : R = Flx_FlxM_qmul(q, R, p,pi);
2489 : }
2490 192 : M = Flx_halfres_pre(x, y, p, pi, &x, &y, &res);
2491 192 : if (!res) return gc_ulong(av, 0);
2492 192 : R = FlxM_mul2(M, R, p, pi);
2493 192 : gerepileall(av,3,&x,&y,&R);
2494 : }
2495 53 : res1 = Flx_extresultant_basecase(x,y,p,pi,&u,&v);
2496 53 : if (!res1) return gc_ulong(av, 0);
2497 53 : *ptU = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi), res, p, pi);
2498 53 : *ptV = Flx_Fl_mul_pre(Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi), res, p, pi);
2499 53 : gerepileall(av, 2, ptU, ptV);
2500 53 : return Fl_mul(res1,res,p);
2501 : }
2502 :
2503 : ulong
2504 53 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2505 53 : { return Flx_extresultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptU, ptV); }
2506 :
2507 : /* allow pi = 0 (SMALL_ULONG) */
2508 : ulong
2509 43789989 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2510 : {
2511 43789989 : ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2512 :
2513 43789989 : if (lx == 1) return 0;
2514 41035826 : x++;
2515 41035826 : if (pi)
2516 : {
2517 : LOCAL_OVERFLOW;
2518 : LOCAL_HIREMAINDER;
2519 40973611 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2520 97815171 : while (++i < lx)
2521 : {
2522 56841560 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2523 56841560 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2524 : }
2525 81118 : return v1? remlll_pre(v1, h1, l1, p, pi)
2526 41054729 : : remll_pre(h1, l1, p, pi);
2527 : }
2528 : else
2529 : {
2530 62215 : l1 = x[i] * y[i];
2531 30928708 : while (++i < lx) { l1 += x[i] * y[i]; if (l1 & HIGHBIT) l1 %= p; }
2532 62215 : return l1 % p;
2533 : }
2534 : }
2535 :
2536 : /* allow pi = 0 (SMALL_ULONG) */
2537 : ulong
2538 100665774 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2539 : {
2540 100665774 : long i, n = degpol(x);
2541 : ulong t;
2542 100668293 : if (n <= 0) return n? 0: x[2];
2543 32935350 : if (n > 15)
2544 : {
2545 180435 : pari_sp av = avma;
2546 180435 : GEN v = Fl_powers_pre(y, n, p, pi);
2547 180428 : return gc_ulong(av, Flx_eval_powers_pre(x, v, p, pi));
2548 : }
2549 32754915 : i = n+2; t = x[i];
2550 32754915 : if (pi)
2551 : {
2552 123084343 : for (i--; i>=2; i--) t = Fl_addmul_pre(uel(x, i), t, y, p, pi);
2553 31645429 : return t;
2554 : }
2555 2675827 : for (i--; i>=2; i--) t = (t * y + x[i]) % p;
2556 1118399 : return t %= p;
2557 : }
2558 : ulong
2559 20383452 : Flx_eval(GEN x, ulong y, ulong p)
2560 20383452 : { return Flx_eval_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2561 :
2562 : ulong
2563 3255 : Flv_prod_pre(GEN x, ulong p, ulong pi)
2564 : {
2565 3255 : pari_sp ltop = avma;
2566 : GEN v;
2567 3255 : long i,k,lx = lg(x);
2568 3255 : if (lx == 1) return 1UL;
2569 3255 : if (lx == 2) return uel(x,1);
2570 3003 : v = cgetg(1+(lx << 1), t_VECSMALL);
2571 3003 : k = 1;
2572 28593 : for (i=1; i<lx-1; i+=2)
2573 25590 : uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2574 3003 : if (i < lx) uel(v,k++) = uel(x,i);
2575 13529 : while (k > 2)
2576 : {
2577 10526 : lx = k; k = 1;
2578 36116 : for (i=1; i<lx-1; i+=2)
2579 25590 : uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2580 10526 : if (i < lx) uel(v,k++) = uel(v,i);
2581 : }
2582 3003 : return gc_ulong(ltop, uel(v,1));
2583 : }
2584 :
2585 : ulong
2586 0 : Flv_prod(GEN v, ulong p)
2587 : {
2588 0 : return Flv_prod_pre(v, p, get_Fl_red(p));
2589 : }
2590 :
2591 : GEN
2592 0 : FlxV_prod(GEN V, ulong p)
2593 : {
2594 : struct _Flxq D;
2595 0 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2596 0 : return gen_product(V, (void *)&D, &_Flx_mul);
2597 : }
2598 :
2599 : /* compute prod (x - a[i]) */
2600 : GEN
2601 737357 : Flv_roots_to_pol(GEN a, ulong p, long vs)
2602 : {
2603 : struct _Flxq D;
2604 737357 : long i,k,lx = lg(a);
2605 : GEN p1;
2606 737357 : if (lx == 1) return pol1_Flx(vs);
2607 737357 : p1 = cgetg(lx, t_VEC);
2608 11905092 : for (k=1,i=1; i<lx-1; i+=2)
2609 11166513 : gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2610 11167985 : Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2611 737107 : if (i < lx)
2612 58112 : gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2613 737092 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2614 737091 : setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2615 : }
2616 :
2617 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2618 : INLINE void
2619 18869285 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2620 : {
2621 18869285 : pari_sp av = avma;
2622 18869285 : long n = lg(w), i;
2623 : ulong u;
2624 : GEN c;
2625 :
2626 18869285 : if (n == 1) return;
2627 18869285 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2628 80475387 : for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2629 19029248 : i = n-1; u = Fl_inv(c[i], p);
2630 80769011 : for ( ; i > 1; --i)
2631 : {
2632 61694185 : ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2633 61672173 : u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2634 : }
2635 19074826 : v[1] = u; set_avma(av);
2636 : }
2637 :
2638 : void
2639 18260094 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2640 :
2641 : GEN
2642 10848 : Flv_inv_pre(GEN w, ulong p, ulong pi)
2643 10848 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2644 :
2645 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2646 : INLINE void
2647 49724 : Flv_inv_indir(GEN w, GEN v, ulong p)
2648 : {
2649 49724 : pari_sp av = avma;
2650 49724 : long n = lg(w), i;
2651 : ulong u;
2652 : GEN c;
2653 :
2654 49724 : if (n == 1) return;
2655 49724 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2656 1718434 : for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2657 49721 : i = n-1; u = Fl_inv(c[i], p);
2658 1718439 : for ( ; i > 1; --i)
2659 : {
2660 1668712 : ulong t = Fl_mul(u, c[i-1], p);
2661 1668712 : u = Fl_mul(u, w[i], p); v[i] = t;
2662 : }
2663 49727 : v[1] = u; set_avma(av);
2664 : }
2665 : static void
2666 635845 : Flv_inv_i(GEN v, GEN w, ulong p)
2667 : {
2668 635845 : if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2669 586121 : else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2670 635848 : }
2671 : void
2672 12017 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2673 : GEN
2674 623829 : Flv_inv(GEN w, ulong p)
2675 623829 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2676 :
2677 : GEN
2678 33022721 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2679 : {
2680 33022721 : long l = lg(a), i;
2681 : GEN a0, z0, z;
2682 33022721 : if (l <= 3)
2683 : {
2684 0 : if (rem) *rem = l == 2? 0: a[2];
2685 0 : return zero_Flx(a[1]);
2686 : }
2687 33022721 : z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2688 32873491 : a0 = a + l-1;
2689 32873491 : z0 = z + l-2; *z0 = *a0--;
2690 32873491 : if (SMALL_ULONG(p))
2691 : {
2692 79707482 : for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2693 : {
2694 59061244 : ulong t = (*a0-- + x * *z0--) % p;
2695 59061244 : *z0 = (long)t;
2696 : }
2697 20646238 : if (rem) *rem = (*a0 + x * *z0) % p;
2698 : }
2699 : else
2700 : {
2701 48245258 : for (i=l-3; i>1; i--)
2702 : {
2703 35992309 : ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2704 36018005 : *z0 = (long)t;
2705 : }
2706 12252949 : if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2707 : }
2708 32900631 : return z;
2709 : }
2710 :
2711 : /* xa, ya = t_VECSMALL */
2712 : static GEN
2713 625031 : Flv_producttree(GEN xa, GEN s, ulong p, ulong pi, long vs)
2714 : {
2715 625031 : long n = lg(xa)-1;
2716 625031 : long m = n==1 ? 1: expu(n-1)+1;
2717 625031 : long i, j, k, ls = lg(s);
2718 625031 : GEN T = cgetg(m+1, t_VEC);
2719 625027 : GEN t = cgetg(ls, t_VEC);
2720 7835045 : for (j=1, k=1; j<ls; k+=s[j++])
2721 7209893 : gel(t, j) = s[j] == 1 ?
2722 7210021 : mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2723 1516152 : mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2724 1516154 : Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2725 625024 : gel(T,1) = t;
2726 2356793 : for (i=2; i<=m; i++)
2727 : {
2728 1731827 : GEN u = gel(T, i-1);
2729 1731827 : long n = lg(u)-1;
2730 1731827 : GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2731 8316041 : for (j=1, k=1; k<n; j++, k+=2)
2732 6584272 : gel(t, j) = Flx_mul_pre(gel(u, k), gel(u, k+1), p, pi);
2733 1731769 : gel(T, i) = t;
2734 : }
2735 624966 : return T;
2736 : }
2737 :
2738 : static GEN
2739 665335 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p, ulong pi)
2740 : {
2741 : long i,j,k;
2742 665335 : long m = lg(T)-1;
2743 665335 : GEN R = cgetg(lg(xa), t_VECSMALL);
2744 665329 : GEN Tp = cgetg(m+1, t_VEC), t;
2745 665326 : gel(Tp, m) = mkvec(P);
2746 2582494 : for (i=m-1; i>=1; i--)
2747 : {
2748 1917165 : GEN u = gel(T, i), v = gel(Tp, i+1);
2749 1917165 : long n = lg(u)-1;
2750 1917165 : t = cgetg(n+1, t_VEC);
2751 9532109 : for (j=1, k=1; k<n; j++, k+=2)
2752 : {
2753 7614932 : gel(t, k) = Flx_rem_pre(gel(v, j), gel(u, k), p, pi);
2754 7614951 : gel(t, k+1) = Flx_rem_pre(gel(v, j), gel(u, k+1), p, pi);
2755 : }
2756 1917177 : gel(Tp, i) = t;
2757 : }
2758 : {
2759 665329 : GEN u = gel(T, i+1), v = gel(Tp, i+1);
2760 665329 : long n = lg(u)-1;
2761 8947760 : for (j=1, k=1; j<=n; j++)
2762 : {
2763 8282427 : long c, d = degpol(gel(u,j));
2764 18330278 : for (c=1; c<=d; c++, k++) R[k] = Flx_eval_pre(gel(v, j), xa[k], p, pi);
2765 : }
2766 665333 : return gc_const((pari_sp)R, R);
2767 : }
2768 : }
2769 :
2770 : static GEN
2771 1386411 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, ulong pi, long vs)
2772 : {
2773 1386411 : pari_sp av = avma;
2774 1386411 : long m = lg(T)-1;
2775 1386411 : long i, j, k, ls = lg(s);
2776 1386411 : GEN Tp = cgetg(m+1, t_VEC);
2777 1386057 : GEN t = cgetg(ls, t_VEC);
2778 24941764 : for (j=1, k=1; j<ls; k+=s[j++])
2779 23555952 : if (s[j]==2)
2780 : {
2781 6936532 : ulong a = Fl_mul(ya[k], R[k], p);
2782 6936140 : ulong b = Fl_mul(ya[k+1], R[k+1], p);
2783 6942342 : gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2784 6936507 : Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2785 6939905 : gel(t, j) = Flx_renormalize(gel(t, j), 4);
2786 : }
2787 : else
2788 16619420 : gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2789 1385812 : gel(Tp, 1) = t;
2790 6387303 : for (i=2; i<=m; i++)
2791 : {
2792 5001333 : GEN u = gel(T, i-1);
2793 5001333 : GEN t = cgetg(lg(gel(T,i)), t_VEC);
2794 4998709 : GEN v = gel(Tp, i-1);
2795 4998709 : long n = lg(v)-1;
2796 27118401 : for (j=1, k=1; k<n; j++, k+=2)
2797 22109880 : gel(t, j) = Flx_add(Flx_mul_pre(gel(u, k), gel(v, k+1), p, pi),
2798 22116910 : Flx_mul_pre(gel(u, k+1), gel(v, k), p, pi), p);
2799 5001491 : gel(Tp, i) = t;
2800 : }
2801 1385970 : return gerepileuptoleaf(av, gmael(Tp,m,1));
2802 : }
2803 :
2804 : GEN
2805 0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2806 : {
2807 0 : pari_sp av = avma;
2808 0 : GEN s = producttree_scheme(lg(xa)-1);
2809 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2810 0 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2811 0 : return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p, pi));
2812 : }
2813 :
2814 : static GEN
2815 2471 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p, ulong pi)
2816 45248 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p, pi)) }
2817 :
2818 : GEN
2819 2471 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2820 : {
2821 2471 : pari_sp av = avma;
2822 2471 : GEN s = producttree_scheme(lg(xa)-1);
2823 2471 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2824 2471 : GEN T = Flv_producttree(xa, s, p, pi, P[1]);
2825 2471 : return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p, pi));
2826 : }
2827 :
2828 : GEN
2829 368455 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2830 : {
2831 368455 : pari_sp av = avma;
2832 368455 : GEN s = producttree_scheme(lg(xa)-1);
2833 368461 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2834 368461 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2835 368461 : long m = lg(T)-1;
2836 368461 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2837 368460 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2838 368462 : return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, pi, vs));
2839 : }
2840 :
2841 : GEN
2842 101105 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2843 : {
2844 101105 : pari_sp av = avma;
2845 101105 : GEN s = producttree_scheme(lg(xa)-1);
2846 101105 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2847 101105 : GEN T = Flv_producttree(xa, s, p, pi, vs);
2848 101102 : long i, m = lg(T)-1, l = lg(ya)-1;
2849 101102 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2850 101103 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p, pi), p);
2851 101103 : GEN M = cgetg(l+1, t_VEC);
2852 1118859 : for (i=1; i<=l; i++)
2853 1017759 : gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, pi, vs);
2854 101100 : return gerepileupto(av, M);
2855 : }
2856 :
2857 : GEN
2858 152995 : Flv_invVandermonde(GEN L, ulong den, ulong p)
2859 : {
2860 152995 : pari_sp av = avma;
2861 152995 : long i, n = lg(L);
2862 : GEN M, R;
2863 152995 : GEN s = producttree_scheme(n-1);
2864 152995 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2865 152995 : GEN tree = Flv_producttree(L, s, p, pi, 0);
2866 152995 : long m = lg(tree)-1;
2867 152995 : GEN T = gmael(tree, m, 1);
2868 152995 : R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p, pi), p);
2869 152995 : if (den!=1) R = Flv_Fl_mul(R, den, p);
2870 152995 : M = cgetg(n, t_MAT);
2871 600537 : for (i = 1; i < n; i++)
2872 : {
2873 447542 : GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2874 447542 : gel(M,i) = Flx_to_Flv(P, n-1);
2875 : }
2876 152995 : return gerepilecopy(av, M);
2877 : }
2878 :
2879 : /***********************************************************************/
2880 : /** Flxq **/
2881 : /***********************************************************************/
2882 : /* Flxq objects are Flx modulo another Flx called q. */
2883 :
2884 : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2885 : GEN
2886 189090451 : Flxq_mul_pre(GEN x,GEN y,GEN T,ulong p,ulong pi)
2887 189090451 : { return Flx_rem_pre(Flx_mul_pre(x,y,p,pi),T,p,pi); }
2888 : GEN
2889 13192904 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2890 13192904 : { return Flxq_mul_pre(x,y,T,p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2891 :
2892 : GEN
2893 277144544 : Flxq_sqr_pre(GEN x,GEN T,ulong p,ulong pi)
2894 277144544 : { return Flx_rem_pre(Flx_sqr_pre(x, p,pi), T, p,pi); }
2895 : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2896 : GEN
2897 2759482 : Flxq_sqr(GEN x,GEN T,ulong p)
2898 2759482 : { return Flxq_sqr_pre(x,T,p,SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2899 :
2900 : static GEN
2901 1551566 : _Flxq_red(void *E, GEN x)
2902 1551566 : { struct _Flxq *s = (struct _Flxq *)E;
2903 1551566 : return Flx_rem_pre(x, s->T, s->p, s->pi); }
2904 : #if 0
2905 : static GEN
2906 : _Flx_sub(void *E, GEN x, GEN y)
2907 : { struct _Flxq *s = (struct _Flxq *)E;
2908 : return Flx_sub(x,y,s->p); }
2909 : #endif
2910 : static GEN
2911 269294050 : _Flxq_sqr(void *data, GEN x)
2912 : {
2913 269294050 : struct _Flxq *D = (struct _Flxq*)data;
2914 269294050 : return Flxq_sqr_pre(x, D->T, D->p, D->pi);
2915 : }
2916 : static GEN
2917 148031527 : _Flxq_mul(void *data, GEN x, GEN y)
2918 : {
2919 148031527 : struct _Flxq *D = (struct _Flxq*)data;
2920 148031527 : return Flxq_mul_pre(x,y, D->T, D->p, D->pi);
2921 : }
2922 : static GEN
2923 22230712 : _Flxq_one(void *data)
2924 : {
2925 22230712 : struct _Flxq *D = (struct _Flxq*)data;
2926 22230712 : return pol1_Flx(get_Flx_var(D->T));
2927 : }
2928 :
2929 : static GEN
2930 22925514 : _Flxq_powu_i(struct _Flxq *D, GEN x, ulong n)
2931 22925514 : { return gen_powu_i(x, n, (void*)D, &_Flxq_sqr, &_Flxq_mul); }
2932 : static GEN
2933 68 : _Flxq_powu(struct _Flxq *D, GEN x, ulong n)
2934 68 : { pari_sp av = avma; return gerepileuptoleaf(av, _Flxq_powu_i(D, x, n)); }
2935 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2936 : GEN
2937 24176270 : Flxq_powu_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
2938 : {
2939 : pari_sp av;
2940 : struct _Flxq D;
2941 24176270 : switch(n)
2942 : {
2943 0 : case 0: return pol1_Flx(get_Flx_var(T));
2944 278017 : case 1: return Flx_copy(x);
2945 972714 : case 2: return Flxq_sqr_pre(x, T, p, pi);
2946 : }
2947 22925539 : av = avma; set_Flxq_pre(&D, T, p, pi);
2948 22925730 : return gerepileuptoleaf(av, _Flxq_powu_i(&D, x, n));
2949 : }
2950 : GEN
2951 488349 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2952 488349 : { return Flxq_powu_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2953 :
2954 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2955 : GEN
2956 23772538 : Flxq_pow_pre(GEN x, GEN n, GEN T, ulong p, ulong pi)
2957 : {
2958 23772538 : pari_sp av = avma;
2959 : struct _Flxq D;
2960 : GEN y;
2961 23772538 : long s = signe(n);
2962 23772538 : if (!s) return pol1_Flx(get_Flx_var(T));
2963 23695946 : if (s < 0) x = Flxq_inv_pre(x,T,p,pi);
2964 23695946 : if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2965 23176680 : set_Flxq_pre(&D, T, p, pi);
2966 23176690 : y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2967 23176624 : return gerepileuptoleaf(av, y);
2968 : }
2969 : GEN
2970 931478 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
2971 931478 : { return Flxq_pow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2972 :
2973 : GEN
2974 28 : Flxq_pow_init_pre(GEN x, GEN n, long k, GEN T, ulong p, ulong pi)
2975 : {
2976 28 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2977 28 : return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2978 : }
2979 : GEN
2980 0 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
2981 0 : { return Flxq_pow_init_pre(x, n, k, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2982 :
2983 : GEN
2984 4393 : Flxq_pow_table_pre(GEN R, GEN n, GEN T, ulong p, ulong pi)
2985 : {
2986 4393 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2987 4393 : return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
2988 : }
2989 : GEN
2990 0 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
2991 0 : { return Flxq_pow_table_pre(R, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2992 :
2993 : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
2994 : * not stack clean. */
2995 : GEN
2996 5411562 : Flxq_invsafe_pre(GEN x, GEN T, ulong p, ulong pi)
2997 : {
2998 5411562 : GEN V, z = Flx_extgcd_pre(get_Flx_mod(T), x, p, pi, NULL, &V);
2999 : ulong iz;
3000 5411673 : if (degpol(z)) return NULL;
3001 5411017 : iz = Fl_inv(uel(z,2), p);
3002 5411009 : return Flx_Fl_mul_pre(V, iz, p, pi);
3003 : }
3004 : GEN
3005 669320 : Flxq_invsafe(GEN x, GEN T, ulong p)
3006 669320 : { return Flxq_invsafe_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3007 :
3008 : GEN
3009 4283232 : Flxq_inv_pre(GEN x, GEN T, ulong p, ulong pi)
3010 : {
3011 4283232 : pari_sp av=avma;
3012 4283232 : GEN U = Flxq_invsafe_pre(x, T, p, pi);
3013 4283229 : if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
3014 4283222 : return gerepileuptoleaf(av, U);
3015 : }
3016 : GEN
3017 335812 : Flxq_inv(GEN x, GEN T, ulong p)
3018 335812 : { return Flxq_inv_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3019 :
3020 : GEN
3021 2416584 : Flxq_div_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
3022 : {
3023 2416584 : pari_sp av = avma;
3024 2416584 : return gerepileuptoleaf(av, Flxq_mul_pre(x,Flxq_inv_pre(y,T,p,pi),T,p,pi));
3025 : }
3026 : GEN
3027 237703 : Flxq_div(GEN x, GEN y, GEN T, ulong p)
3028 237703 : { return Flxq_div_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3029 :
3030 : GEN
3031 22231002 : Flxq_powers_pre(GEN x, long l, GEN T, ulong p, ulong pi)
3032 : {
3033 22231002 : int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
3034 22228215 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3035 22226581 : return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
3036 : }
3037 : GEN
3038 232142 : Flxq_powers(GEN x, long l, GEN T, ulong p)
3039 232142 : { return Flxq_powers_pre(x, l, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3040 :
3041 : GEN
3042 170603 : Flxq_matrix_pow_pre(GEN y, long n, long m, GEN P, ulong l, ulong li)
3043 170603 : { return FlxV_to_Flm(Flxq_powers_pre(y,m-1,P,l,li),n); }
3044 : GEN
3045 399 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
3046 399 : { return Flxq_matrix_pow_pre(y, n, m, P, l, SMALL_ULONG(l)? 0: get_Fl_red(l)); }
3047 :
3048 : GEN
3049 13712252 : Flx_Frobenius_pre(GEN T, ulong p, ulong pi)
3050 13712252 : { return Flxq_powu_pre(polx_Flx(get_Flx_var(T)), p, T, p, pi); }
3051 : GEN
3052 86486 : Flx_Frobenius(GEN T, ulong p)
3053 86486 : { return Flx_Frobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3054 :
3055 : GEN
3056 86529 : Flx_matFrobenius_pre(GEN T, ulong p, ulong pi)
3057 : {
3058 86529 : long n = get_Flx_degree(T);
3059 86529 : return Flxq_matrix_pow_pre(Flx_Frobenius_pre(T, p, pi), n, n, T, p, pi);
3060 : }
3061 : GEN
3062 0 : Flx_matFrobenius(GEN T, ulong p)
3063 0 : { return Flx_matFrobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3064 :
3065 : static GEN
3066 12811908 : Flx_blocks_Flm(GEN P, long n, long m)
3067 : {
3068 12811908 : GEN z = cgetg(m+1,t_MAT);
3069 12811269 : long i,j, k=2, l = lg(P);
3070 36709691 : for(i=1; i<=m; i++)
3071 : {
3072 23901811 : GEN zi = cgetg(n+1,t_VECSMALL);
3073 23898422 : gel(z,i) = zi;
3074 110933067 : for(j=1; j<=n; j++)
3075 87034645 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3076 : }
3077 12807880 : return z;
3078 : }
3079 :
3080 : GEN
3081 516905 : Flx_blocks(GEN P, long n, long m)
3082 : {
3083 516905 : GEN z = cgetg(m+1,t_VEC);
3084 516611 : long i,j, k=2, l = lg(P);
3085 1547906 : for(i=1; i<=m; i++)
3086 : {
3087 1031620 : GEN zi = cgetg(n+2,t_VECSMALL);
3088 1030795 : zi[1] = P[1];
3089 1030795 : gel(z,i) = zi;
3090 6462312 : for(j=2; j<n+2; j++)
3091 5431517 : uel(zi, j) = k==l ? 0 : uel(P,k++);
3092 1030795 : zi = Flx_renormalize(zi, n+2);
3093 : }
3094 516286 : return z;
3095 : }
3096 :
3097 : static GEN
3098 12812679 : FlxV_to_Flm_lg(GEN x, long m, long n)
3099 : {
3100 : long i;
3101 12812679 : GEN y = cgetg(n+1, t_MAT);
3102 60880194 : for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
3103 12810194 : return y;
3104 : }
3105 :
3106 : /* allow pi = 0 (SMALL_ULONG) */
3107 : GEN
3108 13011455 : Flx_FlxqV_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3109 : {
3110 13011455 : pari_sp btop, av = avma;
3111 13011455 : long sv = get_Flx_var(T), m = get_Flx_degree(T);
3112 13011606 : long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
3113 : GEN A, B, C, S, g;
3114 13012376 : if (lQ == 0) return pol0_Flx(sv);
3115 12813468 : if (lQ <= l)
3116 : {
3117 6348896 : n = l;
3118 6348896 : d = 1;
3119 : }
3120 : else
3121 : {
3122 6464572 : n = l-1;
3123 6464572 : d = (lQ+n-1)/n;
3124 : }
3125 12813468 : A = FlxV_to_Flm_lg(x, m, n);
3126 12811798 : B = Flx_blocks_Flm(Q, n, d);
3127 12811002 : C = gerepileupto(av, Flm_mul(A, B, p));
3128 12813923 : g = gel(x, l);
3129 12813923 : if (pi && SMALL_ULONG(p)) pi = 0;
3130 12813923 : T = Flx_get_red_pre(T, p, pi);
3131 12813552 : btop = avma;
3132 12813552 : S = Flv_to_Flx(gel(C, d), sv);
3133 23906113 : for (i = d-1; i>0; i--)
3134 : {
3135 11093654 : S = Flx_add(Flxq_mul_pre(S, g, T, p, pi), Flv_to_Flx(gel(C,i), sv), p);
3136 11093430 : if (gc_needed(btop,1))
3137 0 : S = gerepileuptoleaf(btop, S);
3138 : }
3139 12812459 : return gerepileuptoleaf(av, S);
3140 : }
3141 : GEN
3142 5082 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
3143 5082 : { return Flx_FlxqV_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3144 :
3145 : /* allow pi = 0 (SMALL_ULONG) */
3146 : GEN
3147 2405483 : Flx_Flxq_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3148 : {
3149 2405483 : pari_sp av = avma;
3150 : GEN z, V;
3151 2405483 : long d = degpol(Q), rtd;
3152 2405486 : if (d < 0) return pol0_Flx(get_Flx_var(T));
3153 2405395 : rtd = (long) sqrt((double)d);
3154 2405395 : T = Flx_get_red_pre(T, p, pi);
3155 2405410 : V = Flxq_powers_pre(x, rtd, T, p, pi);
3156 2405453 : z = Flx_FlxqV_eval_pre(Q, V, T, p, pi);
3157 2405380 : return gerepileupto(av, z);
3158 : }
3159 : GEN
3160 790011 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
3161 790011 : { return Flx_Flxq_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3162 :
3163 : /* allow pi = 0 (SMALL_ULONG) */
3164 : GEN
3165 0 : FlxC_FlxqV_eval_pre(GEN x, GEN v, GEN T, ulong p, ulong pi)
3166 0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval_pre(gel(x,i), v, T, p, pi)) }
3167 : GEN
3168 0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
3169 0 : { return FlxC_FlxqV_eval_pre(x, v, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3170 :
3171 : /* allow pi = 0 (SMALL_ULONG) */
3172 : GEN
3173 0 : FlxC_Flxq_eval_pre(GEN x, GEN F, GEN T, ulong p, ulong pi)
3174 : {
3175 0 : long d = brent_kung_optpow(get_Flx_degree(T)-1,lg(x)-1,1);
3176 0 : GEN Fp = Flxq_powers_pre(F, d, T, p, pi);
3177 0 : return FlxC_FlxqV_eval_pre(x, Fp, T, p, pi);
3178 : }
3179 : GEN
3180 0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
3181 0 : { return FlxC_Flxq_eval_pre(x, F, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3182 :
3183 : #if 0
3184 : static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
3185 : _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
3186 : #endif
3187 :
3188 : static GEN
3189 46249 : Flxq_autpow_sqr(void *E, GEN x)
3190 : {
3191 46249 : struct _Flxq *D = (struct _Flxq*)E;
3192 46249 : return Flx_Flxq_eval_pre(x, x, D->T, D->p, D->pi);
3193 : }
3194 : static GEN
3195 20696 : Flxq_autpow_msqr(void *E, GEN x)
3196 : {
3197 20696 : struct _Flxq *D = (struct _Flxq*)E;
3198 20696 : return Flx_FlxqV_eval_pre(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p, D->pi);
3199 : }
3200 :
3201 : GEN
3202 67489 : Flxq_autpow_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3203 : {
3204 67489 : pari_sp av = avma;
3205 : struct _Flxq D;
3206 : long d;
3207 67489 : if (n==0) return Flx_rem_pre(polx_Flx(x[1]), T, p, pi);
3208 67482 : if (n==1) return Flx_rem_pre(x, T, p, pi);
3209 31377 : set_Flxq_pre(&D, T, p, pi);
3210 31377 : d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
3211 31377 : D.aut = Flxq_powers_pre(x, d, T, p, D.pi);
3212 31377 : x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
3213 31377 : return gerepilecopy(av, x);
3214 : }
3215 : GEN
3216 7 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
3217 7 : { return Flxq_autpow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3218 :
3219 : GEN
3220 1667 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
3221 : {
3222 1667 : long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
3223 : ulong i, pi;
3224 1667 : pari_sp av = avma;
3225 1667 : GEN xp, V = cgetg(l+2,t_VEC);
3226 1667 : gel(V,1) = polx_Flx(vT); if (l==0) return V;
3227 1667 : gel(V,2) = gcopy(x); if (l==1) return V;
3228 1667 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3229 1667 : T = Flx_get_red_pre(T, p, pi);
3230 1667 : d = brent_kung_optpow(dT-1, l-1, 1);
3231 1667 : xp = Flxq_powers_pre(x, d, T, p, pi);
3232 6998 : for(i = 3; i < l+2; i++)
3233 5331 : gel(V,i) = Flx_FlxqV_eval_pre(gel(V,i-1), xp, T, p, pi);
3234 1667 : return gerepilecopy(av, V);
3235 : }
3236 :
3237 : static GEN
3238 112451 : Flxq_autsum_mul(void *E, GEN x, GEN y)
3239 : {
3240 112451 : struct _Flxq *D = (struct _Flxq*)E;
3241 112451 : GEN T = D->T;
3242 112451 : ulong p = D->p, pi = D->pi;
3243 112451 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3244 112451 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3245 112451 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3246 112451 : GEN V2 = Flxq_powers_pre(phi2, d, T, p, pi);
3247 112451 : GEN phi3 = Flx_FlxqV_eval_pre(phi1, V2, T, p, pi);
3248 112451 : GEN aphi = Flx_FlxqV_eval_pre(a1, V2, T, p, pi);
3249 112451 : GEN a3 = Flxq_mul_pre(aphi, a2, T, p, pi);
3250 112451 : return mkvec2(phi3, a3);
3251 : }
3252 : static GEN
3253 105089 : Flxq_autsum_sqr(void *E, GEN x)
3254 105089 : { return Flxq_autsum_mul(E, x, x); }
3255 :
3256 : static GEN
3257 98754 : Flxq_autsum_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3258 : {
3259 98754 : pari_sp av = avma;
3260 98754 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3261 98754 : x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
3262 98754 : return gerepilecopy(av, x);
3263 : }
3264 : GEN
3265 0 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
3266 0 : { return Flxq_autsum_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3267 :
3268 : static GEN
3269 763246 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
3270 : {
3271 763246 : struct _Flxq *D = (struct _Flxq*)E;
3272 763246 : GEN T = D->T;
3273 763246 : ulong p = D->p, pi = D->pi;
3274 763246 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3275 763246 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3276 763246 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3277 763263 : GEN V1 = Flxq_powers_pre(phi1, d, T, p, pi);
3278 763221 : GEN phi3 = Flx_FlxqV_eval_pre(phi2, V1, T, p, pi);
3279 763212 : GEN aphi = Flx_FlxqV_eval_pre(a2, V1, T, p, pi);
3280 763235 : GEN a3 = Flx_add(a1, aphi, p);
3281 763245 : return mkvec2(phi3, a3);
3282 : }
3283 :
3284 : static GEN
3285 636212 : Flxq_auttrace_sqr(void *E, GEN x)
3286 636212 : { return Flxq_auttrace_mul(E, x, x); }
3287 :
3288 : GEN
3289 936330 : Flxq_auttrace_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3290 : {
3291 936330 : pari_sp av = avma;
3292 : struct _Flxq D;
3293 936330 : set_Flxq_pre(&D, T, p, pi);
3294 936332 : x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
3295 936310 : return gerepilecopy(av, x);
3296 : }
3297 : GEN
3298 0 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
3299 0 : { return Flxq_auttrace_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3300 :
3301 : static long
3302 393864 : bounded_order(ulong p, GEN b, long k)
3303 : {
3304 393864 : GEN a = modii(utoipos(p), b);
3305 : long i;
3306 809293 : for(i = 1; i < k; i++)
3307 : {
3308 514822 : if (equali1(a)) return i;
3309 415429 : a = modii(muliu(a,p),b);
3310 : }
3311 294471 : return 0;
3312 : }
3313 :
3314 : /* n = (p^d-a)\b
3315 : * b = bb*p^vb
3316 : * p^k = 1 [bb]
3317 : * d = m*k+r+vb
3318 : * u = (p^k-1)/bb;
3319 : * v = (p^(r+vb)-a)/b;
3320 : * w = (p^(m*k)-1)/(p^k-1)
3321 : * n = p^r*w*u+v
3322 : * w*u = p^vb*(p^(m*k)-1)/b
3323 : * n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b */
3324 : static GEN
3325 22732434 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p, ulong pi)
3326 : {
3327 22732434 : pari_sp av=avma;
3328 22732434 : long d = get_Flx_degree(T);
3329 22732435 : GEN an = absi_shallow(n), z, q;
3330 22732436 : if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow_pre(x, n, T, p, pi);
3331 394226 : q = powuu(p, d);
3332 394226 : if (dvdii(q, n))
3333 : {
3334 314 : long vn = logint(an, utoipos(p));
3335 314 : GEN autvn = vn==1 ? aut: Flxq_autpow_pre(aut,vn,T,p,pi);
3336 314 : z = Flx_Flxq_eval_pre(x,autvn,T,p,pi);
3337 : } else
3338 : {
3339 393912 : GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
3340 : GEN bb, u, v, autk;
3341 393912 : long vb = Z_lvalrem(b,p,&bb);
3342 393912 : long m, r, k = is_pm1(bb)? 1: bounded_order(p,bb,d);
3343 393912 : if (!k || d-vb < k) return Flxq_pow_pre(x,n, T,p,pi);
3344 99434 : m = (d-vb)/k; r = (d-vb)%k;
3345 99434 : u = diviiexact(subiu(powuu(p,k),1),bb);
3346 99434 : v = diviiexact(subii(powuu(p,r+vb),a),b);
3347 99434 : autk = k==1 ? aut: Flxq_autpow_pre(aut,k,T,p,pi);
3348 99434 : if (r)
3349 : {
3350 487 : GEN autr = r==1 ? aut: Flxq_autpow_pre(aut,r,T,p,pi);
3351 487 : z = Flx_Flxq_eval_pre(x,autr,T,p,pi);
3352 98947 : } else z = x;
3353 99434 : if (m > 1) z = gel(Flxq_autsum_pre(mkvec2(autk, z), m, T, p, pi), 2);
3354 99434 : if (!is_pm1(u)) z = Flxq_pow_pre(z, u, T, p, pi);
3355 99434 : if (signe(v)) z = Flxq_mul_pre(z, Flxq_pow_pre(x, v, T, p, pi), T, p, pi);
3356 : }
3357 99748 : return gerepileupto(av,signe(n)>0 ? z : Flxq_inv_pre(z,T,p,pi));
3358 : }
3359 :
3360 : static GEN
3361 22725056 : _Flxq_pow(void *data, GEN x, GEN n)
3362 : {
3363 22725056 : struct _Flxq *D = (struct _Flxq*)data;
3364 22725056 : return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p, D->pi);
3365 : }
3366 :
3367 : static GEN
3368 5465 : _Flxq_rand(void *data)
3369 : {
3370 5465 : pari_sp av=avma;
3371 5465 : struct _Flxq *D = (struct _Flxq*)data;
3372 : GEN z;
3373 : do
3374 : {
3375 5478 : set_avma(av);
3376 5478 : z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3377 5478 : } while (lgpol(z)==0);
3378 5465 : return z;
3379 : }
3380 :
3381 : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3382 : static GEN
3383 35443 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3384 : {
3385 35443 : pari_sp av = avma;
3386 : GEN q,n_q,ord,ordp, op;
3387 :
3388 35443 : if (a == 1UL) return gen_0;
3389 : /* p > 2 */
3390 :
3391 35443 : ordp = utoi(p - 1);
3392 35443 : ord = get_arith_Z(o);
3393 35443 : if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3394 35443 : if (a == p - 1) /* -1 */
3395 7739 : return gerepileuptoint(av, shifti(ord,-1));
3396 27704 : ordp = gcdii(ordp, ord);
3397 27704 : op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3398 :
3399 27704 : q = NULL;
3400 27704 : if (T)
3401 : { /* we want < g > = Fp^* */
3402 27704 : if (!equalii(ord,ordp)) {
3403 11906 : q = diviiexact(ord,ordp);
3404 11906 : g = Flxq_pow(g,q,T,p);
3405 : }
3406 : }
3407 27704 : n_q = Fp_log(utoi(a), utoipos(uel(g,2)), op, utoipos(p));
3408 27704 : if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
3409 27704 : if (q) n_q = mulii(q, n_q);
3410 27704 : return gerepileuptoint(av, n_q);
3411 : }
3412 :
3413 : static GEN
3414 519071 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3415 : {
3416 519071 : struct _Flxq *f = (struct _Flxq *)E;
3417 519071 : GEN T = f->T;
3418 519071 : ulong p = f->p;
3419 519071 : long d = get_Flx_degree(T);
3420 519071 : if (Flx_equal1(a)) return gen_0;
3421 359307 : if (Flx_equal(a,g)) return gen_1;
3422 174231 : if (!degpol(a))
3423 35443 : return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3424 138788 : if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3425 138760 : return NULL;
3426 28 : return Flxq_log_index(a, g, ord, T, p);
3427 : }
3428 :
3429 : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3430 :
3431 : const struct bb_group *
3432 281012 : get_Flxq_star(void **E, GEN T, ulong p)
3433 : {
3434 281012 : struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3435 281012 : e->T = T; e->p = p; e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3436 281012 : e->aut = Flx_Frobenius_pre(T, p, e->pi);
3437 281012 : *E = (void*)e; return &Flxq_star;
3438 : }
3439 :
3440 : GEN
3441 97253 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3442 : {
3443 : void *E;
3444 97253 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3445 97253 : return gen_order(a,ord,E,S);
3446 : }
3447 :
3448 : GEN
3449 164497 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3450 : {
3451 : void *E;
3452 164497 : pari_sp av = avma;
3453 164497 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3454 164497 : GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3455 164497 : if (lg(F) > 1 && Flxq_log_use_index(veclast(F), T, p))
3456 24381 : v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3457 164497 : return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
3458 : }
3459 :
3460 : static GEN
3461 292885 : Flxq_sumautsum_sqr(void *E, GEN xzd)
3462 : {
3463 292885 : struct _Flxq *D = (struct _Flxq*)E;
3464 292885 : pari_sp av = avma;
3465 : GEN xi, zeta, delta, xi2, zeta2, delta2, temp, xipow;
3466 292885 : GEN T = D->T;
3467 292885 : ulong d, p = D-> p, pi = D->pi;
3468 292885 : xi = gel(xzd, 1); zeta = gel(xzd, 2); delta = gel(xzd, 3);
3469 :
3470 292885 : d = brent_kung_optpow(get_Flx_degree(T)-1,3,1);
3471 292885 : xipow = Flxq_powers_pre(xi, d, T, p, pi);
3472 :
3473 292885 : xi2 = Flx_FlxqV_eval_pre(xi, xipow, T, p, pi);
3474 292885 : zeta2 = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(zeta, xipow, T, p, pi), T, p, pi);
3475 292885 : temp = Flxq_mul_pre(zeta, Flx_FlxqV_eval_pre(delta, xipow, T, p, pi), T, p, pi);
3476 292885 : delta2 = Flx_add(delta, temp, p);
3477 292885 : return gerepilecopy(av, mkvec3(xi2, zeta2, delta2));
3478 : }
3479 :
3480 : static GEN
3481 40683 : Flxq_sumautsum_msqr(void *E, GEN xzd)
3482 : {
3483 40683 : struct _Flxq *D = (struct _Flxq*)E;
3484 40683 : pari_sp av = avma;
3485 : GEN xii, zetai, deltai, xzd2;
3486 40683 : GEN T = D->T, xi0pow = gel(D->aut, 1), zeta0 = gel(D->aut, 2);
3487 40683 : ulong p = D-> p, pi = D->pi;
3488 40683 : xzd2 = Flxq_sumautsum_sqr(E, xzd);
3489 40683 : xii = Flx_FlxqV_eval_pre(gel(xzd2, 1), xi0pow, T, p, pi);
3490 40683 : zetai = Flxq_mul_pre(zeta0, Flx_FlxqV_eval_pre(gel(xzd2, 2), xi0pow, T, p, pi), T, p, pi);
3491 40683 : deltai = Flx_add(gel(xzd2, 3), zetai, p);
3492 :
3493 40683 : return gerepilecopy(av, mkvec3(xii, zetai, deltai));
3494 : }
3495 :
3496 : /*returns a + a^(1+s) + a^(1+s+2s) + ... + a^(1+s+...+is)
3497 : where ax = [a,s] with s an automorphism */
3498 : static GEN
3499 208737 : Flxq_sumautsum_pre(GEN ax, long i, GEN T, ulong p, ulong pi) {
3500 208737 : pari_sp av = avma;
3501 : GEN a, xi, zeta, vec, res;
3502 : struct _Flxq D;
3503 : ulong d;
3504 208737 : D.T = Flx_get_red(T, p); D.p = p; D.pi = pi;
3505 208737 : a = gel(ax, 1); xi = gel(ax,2);
3506 208737 : d = brent_kung_optpow(get_Flx_degree(T)-1,2*(hammingl(i)-1),1);
3507 208737 : zeta = Flx_Flxq_eval_pre(a, xi, T, p, pi);
3508 208737 : D.aut = mkvec2(Flxq_powers_pre(xi, d, T, p, pi), zeta);
3509 :
3510 208737 : vec = gen_powu_fold(mkvec3(xi, zeta, zeta), i, (void *)&D, Flxq_sumautsum_sqr, Flxq_sumautsum_msqr);
3511 208737 : res = Flxq_mul_pre(a, Flx_add(pol1_Flx(get_Flx_var(T)), gel(vec, 3), p), T, p, pi);
3512 :
3513 208737 : return gerepilecopy(av, res);
3514 : }
3515 :
3516 : /*algorithm from
3517 : Doliskani, J., & Schost, E. (2014).
3518 : Taking roots over high extensions of finite fields*/
3519 : static GEN
3520 35707 : Flxq_sqrtl_spec_pre(GEN z, GEN n, GEN T, ulong p, ulong pi, GEN *zetan)
3521 : {
3522 35707 : pari_sp av = avma;
3523 : GEN psn, c, b, new_z, beta, x, y, w, ax, g, zeta;
3524 35707 : long s, l, v = get_Flx_var(T), d = get_Flx_degree(T);
3525 : ulong zeta2, beta2;
3526 35707 : s = itos(Fp_order(utoi(p), stoi(d), n));
3527 35706 : if(s >= d || d % s != 0)
3528 0 : pari_err(e_MISC, "expected p's order mod n to divide the degree of T");
3529 35706 : l = d/s;
3530 35706 : if (!lgpol(z)) return pol0_Flx(get_Flx_var(T));
3531 35706 : T = Flx_get_red(T, p);
3532 35706 : ax = mkvec2(NULL, Flxq_autpow_pre(Flx_Frobenius_pre(T,p,pi), s, T, p,pi));
3533 35707 : psn = diviiexact(subiu(powuu(p, s), 1), n);
3534 : do {
3535 39505 : do c = random_Flx(d, v, p); while (!lgpol(c));
3536 39038 : new_z = Flxq_mul_pre(z, Flxq_pow_pre(c, n, T, p,pi), T, p,pi);
3537 39039 : gel(ax,1) = Flxq_pow_pre(new_z, psn, T, p,pi);
3538 :
3539 : /*If l == 2, b has to be 1 + a^((p^s-1)/n)*/
3540 39038 : if(l == 2) y = gel(ax, 1);
3541 1200 : else y = Flxq_sumautsum_pre(ax, l-2, T, p, pi);
3542 39038 : b = Flx_Fl_add(y, 1, p);
3543 39038 : } while (!lgpol(b));
3544 :
3545 35706 : x = Flxq_mul_pre(new_z, Flxq_pow_pre(b, n, T, p,pi), T, p,pi);
3546 35706 : if(s == 1) {
3547 35622 : if (degpol(x) > 0) return gc_NULL(av);
3548 35580 : beta2 = Fl_sqrtn(Flx_constant(x), umodiu(n, p), p, &zeta2);
3549 35580 : if (beta2==~0UL) return gc_NULL(av);
3550 35580 : if(zetan) *zetan = monomial_Flx(zeta2, 0, get_Flx_var(T));
3551 35580 : w = Flx_Fl_mul(Flxq_inv_pre(Flxq_mul_pre(b, c, T, p,pi), T, p,pi), beta2, p);
3552 35581 : gerepileall(av, zetan? 2: 1, &w, zetan);
3553 35581 : return w;
3554 : }
3555 84 : g = Flxq_minpoly(x, T, p);
3556 84 : if (degpol(g) != s) return gc_NULL(av);
3557 77 : beta = Flxq_sqrtn(polx_Flx(get_Flx_var(T)), n, g, p, &zeta);
3558 77 : if (!beta) return gc_NULL(av);
3559 :
3560 77 : if(zetan) *zetan = Flx_Flxq_eval(zeta, x, T, p);
3561 77 : beta = Flx_Flxq_eval(beta, x, T, p);
3562 77 : w = Flxq_mul_pre(Flxq_inv_pre(Flxq_mul_pre(b, c, T, p,pi), T, p,pi), beta, T, p,pi);
3563 77 : gerepileall(av, zetan? 2: 1, &w, zetan);
3564 77 : return w;
3565 : }
3566 :
3567 : static GEN
3568 19262 : Flxq_sqrtn_spec_pre(GEN a, GEN n, GEN T, ulong p, ulong pi, GEN q, GEN *zetan)
3569 : {
3570 19262 : pari_sp ltop = avma;
3571 : GEN z, m, u1, u2;
3572 : int is_1;
3573 19262 : if (is_pm1(n))
3574 : {
3575 847 : if (zetan) *zetan = pol1_Flx(get_Flx_var(T));
3576 847 : return signe(n) < 0? Flxq_inv_pre(a, T, p,pi): gcopy(a);
3577 : }
3578 18415 : is_1 = gequal1(a);
3579 18415 : if (is_1 && !zetan) return gcopy(a);
3580 18415 : z = pol1_Flx(get_Flx_var(T));
3581 18415 : m = bezout(n,q,&u1,&u2);
3582 18414 : if (!is_pm1(m))
3583 : {
3584 18414 : GEN F = Z_factor(m);
3585 18415 : long i, j, j2 = 0; /* -Wall */
3586 : GEN y, l;
3587 18415 : pari_sp av1 = avma;
3588 36907 : for (i = nbrows(F); i; i--)
3589 : {
3590 18541 : l = gcoeff(F,i,1);
3591 18541 : j = itos(gcoeff(F,i,2));
3592 18541 : if(zetan) {
3593 188 : a = Flxq_sqrtl_spec_pre(a,l,T,p,pi,&y);
3594 237 : if (!a) return gc_NULL(ltop);
3595 188 : j--;
3596 188 : j2 = j;
3597 : }
3598 18541 : if (!is_1 && j > 0) {
3599 : do
3600 : {
3601 35309 : a = Flxq_sqrtl_spec_pre(a,l,T,p,pi,NULL);
3602 35309 : if (!a) return gc_NULL(ltop);
3603 35260 : } while (--j);
3604 : }
3605 : /*This is below finding a's root,
3606 : so we don't spend time doing this, if a is not n-th root*/
3607 18492 : if(zetan) {
3608 391 : for(; j2>0; j2--) y = Flxq_sqrtl_spec_pre(y, l, T, p,pi,NULL);
3609 181 : z = Flxq_mul_pre(z, y, T, p,pi);
3610 : }
3611 18492 : if (gc_needed(ltop,1))
3612 : { /* n can have lots of prime factors*/
3613 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxq_sqrtn_spec");
3614 0 : gerepileall(av1, zetan? 2: 1, &a, &z);
3615 : }
3616 : }
3617 : }
3618 :
3619 18366 : if (!equalii(m, n))
3620 119 : a = Flxq_pow_pre(a,modii(u1,q), T, p,pi);
3621 18366 : if (zetan)
3622 : {
3623 181 : *zetan = z;
3624 181 : gerepileall(ltop,2,&a,zetan);
3625 : }
3626 : else /* is_1 is 0: a was modified above -> gerepileupto valid */
3627 18185 : a = gerepileupto(ltop, a);
3628 18366 : return a;
3629 : }
3630 :
3631 : GEN
3632 20461 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3633 : {
3634 20461 : if (!lgpol(a))
3635 : {
3636 7 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3637 0 : if (zeta)
3638 0 : *zeta=pol1_Flx(get_Flx_var(T));
3639 0 : return pol0_Flx(get_Flx_var(T));
3640 : }
3641 20454 : else if(p == 2) {
3642 1192 : pari_sp av = avma;
3643 : GEN z;
3644 1192 : z = F2xq_sqrtn(Flx_to_F2x(a), n, Flx_to_F2x(get_FpX_mod(T)), zeta);
3645 1192 : if (!z) return NULL;
3646 1192 : z = F2x_to_Flx(z);
3647 1192 : if (!zeta) return gerepileuptoleaf(av, z);
3648 0 : *zeta=F2x_to_Flx(*zeta);
3649 0 : return gc_all(av, 2, &z,zeta);
3650 : }
3651 : else
3652 : {
3653 : void *E;
3654 19262 : pari_sp av = avma;
3655 19262 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3656 19262 : GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3657 : GEN m, u1, u2, l, zeta2, F, n2, z;
3658 19261 : long i, s, pi, d = get_Flx_degree(T);
3659 19261 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3660 19261 : m = bezout(n,o,&u1,&u2);
3661 19261 : F = Z_factor(m);
3662 41490 : for (i = nbrows(F); i; i--)
3663 : {
3664 22229 : l = gcoeff(F,i,1);
3665 22229 : s = itos(Fp_order(utoi(p), subiu(l, 1), l));
3666 : /*Flxq_sqrtn_spec only works if d > s and s | d
3667 : for those factors of m we use Flxq_sqrtn_spec
3668 : for the other factor we stay with gen_Shanks_sqrtn*/
3669 22229 : if(d <= s || d % s != 0) {
3670 3689 : gcoeff(F,i,2) = gen_0;
3671 : }
3672 18540 : else gcoeff(F,i,2) = stoi(Z_pval(n,l));
3673 : }
3674 19261 : F = factorback(F);
3675 19262 : z = Flxq_sqrtn_spec_pre(a,F,T, p,pi,o,zeta);
3676 19262 : if(!z) return gc_NULL(av);
3677 19213 : n2 = diviiexact(n, F);
3678 19213 : if(!gequal1(n2)) {
3679 3934 : if(zeta) zeta2 = gcopy(*zeta);
3680 3934 : z = gen_Shanks_sqrtn(z, n2, o, zeta, E, S);
3681 3934 : if (!z) return gc_NULL(av);
3682 3934 : if(zeta) *zeta = Flxq_mul_pre(*zeta, zeta2, T, p,pi);
3683 : }
3684 19213 : return gc_all(av, zeta?2:1, &z, zeta);
3685 : }
3686 : }
3687 :
3688 : GEN
3689 232377 : Flxq_sqrt_pre(GEN z, GEN T, ulong p, ulong pi)
3690 : {
3691 232377 : pari_sp av = avma;
3692 : long d;
3693 232377 : if (p==2)
3694 : {
3695 0 : GEN r = F2xq_sqrt(Flx_to_F2x(z), Flx_to_F2x(get_Flx_mod(T)));
3696 0 : return gerepileupto(av, F2x_to_Flx(r));
3697 : }
3698 232377 : d = get_Flx_degree(T);
3699 232377 : if (d==2)
3700 : {
3701 67676 : GEN P = get_Flx_mod(T), s;
3702 67676 : ulong c = uel(P,2), b = uel(P,3), a = uel(P,4);
3703 67676 : ulong y = degpol(z)<1 ? 0: uel(z,3);
3704 67676 : if (a==1 && b==0)
3705 14890 : {
3706 15670 : ulong x = degpol(z)<1 ? Flx_constant(z): uel(z,2);
3707 15670 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), Fl_neg(c, p), p, pi);
3708 15670 : if (!r) return gc_NULL(av);
3709 14890 : s = mkvecsmall3(P[1], uel(r,1), uel(r,2));
3710 : }
3711 : else
3712 : {
3713 52006 : ulong b2 = Fl_halve(b, p), t = Fl_div(b2, a, p);
3714 52006 : ulong D = Fl_sub(Fl_sqr(b2, p), Fl_mul(a, c, p), p);
3715 52006 : ulong x = degpol(z)<1 ? Flx_constant(z): Fl_sub(uel(z,2), Fl_mul(uel(z,3), t, p), p);
3716 52006 : GEN r = Fl2_sqrt_pre(mkvecsmall2(x, y), D, p, pi);
3717 52006 : if (!r) return gc_NULL(av);
3718 49612 : s = mkvecsmall3(P[1], Fl_add(uel(r,1), Fl_mul(uel(r,2),t,p), p), uel(r,2));
3719 : }
3720 64502 : return gerepileuptoleaf(av, Flx_renormalize(s, 4));
3721 : }
3722 164701 : if (lgpol(z)<=1 && odd(d))
3723 : {
3724 11626 : pari_sp av = avma;
3725 11626 : ulong s = Fl_sqrt(Flx_constant(z), p);
3726 11626 : if (s==~0UL) return gc_NULL(av);
3727 11612 : return gerepilecopy(av, Fl_to_Flx(s, get_Flx_var(T)));
3728 : } else
3729 : {
3730 : GEN c, b, new_z, x, y, w, ax;
3731 : ulong p2, beta;
3732 153075 : long v = get_Flx_var(T);
3733 153075 : if (!lgpol(z)) return pol0_Flx(v);
3734 152410 : T = Flx_get_red_pre(T, p, pi);
3735 152410 : ax = mkvec2(NULL, Flx_Frobenius_pre(T, p, pi));
3736 152410 : p2 = p >> 1; /* (p-1) / 2 */
3737 : do {
3738 208209 : do c = random_Flx(d, v, p); while (!lgpol(c));
3739 :
3740 207537 : new_z = Flxq_mul_pre(z, Flxq_sqr_pre(c, T, p, pi), T, p, pi);
3741 207537 : gel(ax, 1) = Flxq_powu_pre(new_z, p2, T, p, pi);
3742 207537 : y = Flxq_sumautsum_pre(ax, d-2, T, p, pi); /* d > 2 */
3743 207537 : b = Flx_Fl_add(y, 1UL, p);
3744 207537 : } while (!lgpol(b));
3745 :
3746 152410 : x = Flxq_mul_pre(new_z, Flxq_sqr_pre(b, T, p, pi), T, p, pi);
3747 152410 : if (degpol(x) > 0) return gc_NULL(av);
3748 145368 : beta = Fl_sqrt_pre(Flx_constant(x), p, pi);
3749 145368 : if (beta==~0UL) return gc_NULL(av);
3750 145368 : w = Flx_Fl_mul(Flxq_inv_pre(Flxq_mul_pre(b, c, T,p,pi), T,p,pi), beta, p);
3751 145368 : return gerepilecopy(av, w);
3752 : }
3753 : }
3754 :
3755 : GEN
3756 232376 : Flxq_sqrt(GEN a, GEN T, ulong p)
3757 232376 : { return Flxq_sqrt_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3758 :
3759 : /* assume T irreducible mod p */
3760 : int
3761 404152 : Flxq_issquare(GEN x, GEN T, ulong p)
3762 : {
3763 404152 : if (lgpol(x) == 0 || p == 2) return 1;
3764 397845 : return krouu(Flxq_norm(x,T,p), p) == 1;
3765 : }
3766 :
3767 : /* assume T irreducible mod p */
3768 : int
3769 0 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3770 : {
3771 : pari_sp av;
3772 : GEN m;
3773 0 : if (n==1) return Flxq_issquare(x, T, p);
3774 0 : if (lgpol(x) == 0 || p == 2) return 1;
3775 0 : av = avma;
3776 0 : m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3777 0 : return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3778 : }
3779 :
3780 : GEN
3781 113589 : Flxq_lroot_fast_pre(GEN a, GEN sqx, GEN T, long p, ulong pi)
3782 : {
3783 113589 : pari_sp av=avma;
3784 113589 : GEN A = Flx_splitting(a,p);
3785 113589 : return gerepileuptoleaf(av, FlxqV_dotproduct_pre(A,sqx,T,p,pi));
3786 : }
3787 : GEN
3788 0 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
3789 0 : { return Flxq_lroot_fast_pre(a, sqx, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3790 :
3791 : GEN
3792 25053 : Flxq_lroot_pre(GEN a, GEN T, long p, ulong pi)
3793 : {
3794 25053 : pari_sp av=avma;
3795 25053 : long n = get_Flx_degree(T), d = degpol(a);
3796 : GEN sqx, V;
3797 25053 : if (n==1) return leafcopy(a);
3798 25053 : if (n==2) return Flxq_powu_pre(a, p, T, p, pi);
3799 25053 : sqx = Flxq_autpow_pre(Flx_Frobenius_pre(T, p, pi), n-1, T, p, pi);
3800 25053 : if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
3801 0 : if (d>=p)
3802 : {
3803 0 : V = Flxq_powers_pre(sqx,p-1,T,p,pi);
3804 0 : return gerepileuptoleaf(av, Flxq_lroot_fast_pre(a,V,T,p,pi));
3805 : } else
3806 0 : return gerepileuptoleaf(av, Flx_Flxq_eval_pre(a,sqx,T,p,pi));
3807 : }
3808 : GEN
3809 0 : Flxq_lroot(GEN a, GEN T, long p)
3810 0 : { return Flxq_lroot_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3811 :
3812 : ulong
3813 443348 : Flxq_norm(GEN x, GEN TB, ulong p)
3814 : {
3815 443348 : GEN T = get_Flx_mod(TB);
3816 443348 : ulong y = Flx_resultant(T, x, p), L = Flx_lead(T);
3817 443348 : if (L==1 || lgpol(x)==0) return y;
3818 0 : return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3819 : }
3820 :
3821 : ulong
3822 4696 : Flxq_trace(GEN x, GEN TB, ulong p)
3823 : {
3824 4696 : pari_sp av = avma;
3825 : ulong t;
3826 4696 : GEN T = get_Flx_mod(TB);
3827 4696 : long n = degpol(T)-1;
3828 4696 : GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3829 4696 : t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3830 4696 : return gc_ulong(av, t);
3831 : }
3832 :
3833 : /*x must be reduced*/
3834 : GEN
3835 3624 : Flxq_charpoly(GEN x, GEN TB, ulong p)
3836 : {
3837 3624 : pari_sp ltop=avma;
3838 3624 : GEN T = get_Flx_mod(TB);
3839 3624 : long vs = evalvarn(fetch_var());
3840 3624 : GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3841 3624 : GEN r = Flx_FlxY_resultant(T, xm1, p);
3842 3624 : r[1] = x[1];
3843 3624 : (void)delete_var(); return gerepileupto(ltop, r);
3844 : }
3845 :
3846 : /* Computing minimal polynomial : */
3847 : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3848 : /* in Algebraic Extensions of Finite Fields' */
3849 :
3850 : /* Let v a linear form, return the linear form z->v(tau*z)
3851 : that is, v*(M_tau) */
3852 :
3853 : static GEN
3854 1693742 : Flxq_transmul_init(GEN tau, GEN T, ulong p, ulong pi)
3855 : {
3856 : GEN bht;
3857 1693742 : GEN h, Tp = get_Flx_red(T, &h);
3858 1693735 : long n = degpol(Tp), vT = Tp[1];
3859 1693722 : GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3860 1693723 : GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3861 1693728 : ft[1] = vT; bt[1] = vT;
3862 1693728 : if (h)
3863 2688 : bht = Flxn_mul_pre(bt, h, n-1, p, pi);
3864 : else
3865 : {
3866 1691040 : GEN bh = Flx_div_pre(Flx_shift(tau, n-1), T, p, pi);
3867 1691025 : bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3868 1691025 : bht[1] = vT;
3869 : }
3870 1693713 : return mkvec3(bt, bht, ft);
3871 : }
3872 :
3873 : static GEN
3874 4086934 : Flxq_transmul(GEN tau, GEN a, long n, ulong p, ulong pi)
3875 : {
3876 4086934 : pari_sp ltop = avma;
3877 : GEN t1, t2, t3, vec;
3878 4086934 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3879 4086934 : if (lgpol(a)==0) return pol0_Flx(a[1]);
3880 4056122 : t2 = Flx_shift(Flx_mul_pre(bt, a, p, pi),1-n);
3881 4055777 : if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
3882 3059866 : t1 = Flx_shift(Flx_mul_pre(ft, a, p, pi),-n);
3883 3059872 : t3 = Flxn_mul_pre(t1, bht, n-1, p, pi);
3884 3059917 : vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3885 3059965 : return gerepileuptoleaf(ltop, vec);
3886 : }
3887 :
3888 : GEN
3889 785102 : Flxq_minpoly_pre(GEN x, GEN T, ulong p, ulong pi)
3890 : {
3891 785102 : pari_sp ltop = avma;
3892 785102 : long vT = get_Flx_var(T), n = get_Flx_degree(T);
3893 : GEN v_x;
3894 785096 : GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3895 785065 : T = Flx_get_red_pre(T, p, pi);
3896 785072 : v_x = Flxq_powers_pre(x, usqrt(2*n), T, p, pi);
3897 1631929 : while (lgpol(tau) != 0)
3898 : {
3899 : long i, j, m, k1;
3900 : GEN M, v, tr, g_prime, c;
3901 846849 : if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3902 846849 : v = random_Flx(n, vT, p);
3903 846882 : tr = Flxq_transmul_init(tau, T, p, pi);
3904 846858 : v = Flxq_transmul(tr, v, n, p, pi);
3905 846871 : m = 2*(n-degpol(g));
3906 846875 : k1 = usqrt(m);
3907 846875 : tr = Flxq_transmul_init(gel(v_x,k1+1), T, p, pi);
3908 846866 : c = cgetg(m+2,t_VECSMALL);
3909 846822 : c[1] = vT;
3910 4086741 : for (i=0; i<m; i+=k1)
3911 : {
3912 3239861 : long mj = minss(m-i, k1);
3913 12657443 : for (j=0; j<mj; j++)
3914 9417139 : uel(c,m+1-(i+j)) = Flx_dotproduct_pre(v, gel(v_x,j+1), p, pi);
3915 3240304 : v = Flxq_transmul(tr, v, n, p, pi);
3916 : }
3917 846880 : c = Flx_renormalize(c, m+2);
3918 : /* now c contains <v,x^i> , i = 0..m-1 */
3919 846880 : M = Flx_halfgcd_pre(monomial_Flx(1, m, vT), c, p, pi);
3920 846895 : g_prime = gmael(M, 2, 2);
3921 846895 : if (degpol(g_prime) < 1) continue;
3922 834880 : g = Flx_mul_pre(g, g_prime, p, pi);
3923 834865 : tau = Flxq_mul_pre(tau, Flx_FlxqV_eval_pre(g_prime, v_x, T,p,pi), T,p,pi);
3924 : }
3925 785041 : g = Flx_normalize(g,p);
3926 785084 : return gerepileuptoleaf(ltop,g);
3927 : }
3928 : GEN
3929 44554 : Flxq_minpoly(GEN x, GEN T, ulong p)
3930 44554 : { return Flxq_minpoly_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3931 :
3932 : GEN
3933 20 : Flxq_conjvec(GEN x, GEN T, ulong p)
3934 : {
3935 20 : long i, l = 1+get_Flx_degree(T);
3936 20 : GEN z = cgetg(l,t_COL);
3937 20 : struct _Flxq D; set_Flxq(&D, T, p);
3938 20 : gel(z,1) = Flx_copy(x);
3939 88 : for (i=2; i<l; i++) gel(z,i) = _Flxq_powu(&D, gel(z,i-1), p);
3940 20 : return z;
3941 : }
3942 :
3943 : GEN
3944 7201 : gener_Flxq(GEN T, ulong p, GEN *po)
3945 : {
3946 7201 : long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
3947 : ulong p_1, pi;
3948 : GEN g, L, L2, o, q, F;
3949 : pari_sp av0, av;
3950 :
3951 7201 : if (f == 1) {
3952 : GEN fa;
3953 28 : o = utoipos(p-1);
3954 28 : fa = Z_factor(o);
3955 28 : L = gel(fa,1);
3956 28 : L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
3957 28 : g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
3958 28 : if (po) *po = mkvec2(o, fa);
3959 28 : return g;
3960 : }
3961 :
3962 7173 : av0 = avma; p_1 = p - 1;
3963 7173 : q = diviuexact(subiu(powuu(p,f), 1), p_1);
3964 :
3965 7173 : L = cgetg(1, t_VECSMALL);
3966 7173 : if (p > 3)
3967 : {
3968 2371 : ulong t = p_1 >> vals(p_1);
3969 2371 : GEN P = gel(factoru(t), 1);
3970 2371 : L = cgetg_copy(P, &i);
3971 3787 : while (--i) L[i] = p_1 / P[i];
3972 : }
3973 7173 : o = factor_pn_1(utoipos(p),f);
3974 7173 : L2 = leafcopy( gel(o, 1) );
3975 19212 : for (i = j = 1; i < lg(L2); i++)
3976 : {
3977 12039 : if (umodui(p_1, gel(L2,i)) == 0) continue;
3978 6488 : gel(L2,j++) = diviiexact(q, gel(L2,i));
3979 : }
3980 7173 : setlg(L2, j); pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3981 7173 : F = Flx_Frobenius_pre(T, p, pi);
3982 17669 : for (av = avma;; set_avma(av))
3983 10496 : {
3984 : GEN tt;
3985 17669 : g = random_Flx(f, vT, p);
3986 17669 : if (degpol(g) < 1) continue;
3987 12101 : if (p == 2) tt = g;
3988 : else
3989 : {
3990 8902 : ulong t = Flxq_norm(g, T, p);
3991 8902 : if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
3992 4742 : tt = Flxq_powu_pre(g, p_1>>1, T, p, pi);
3993 : }
3994 14551 : for (i = 1; i < j; i++)
3995 : {
3996 7378 : GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p, pi);
3997 7378 : if (!degpol(a) && uel(a,2) == p_1) break;
3998 : }
3999 7941 : if (i == j) break;
4000 : }
4001 7173 : if (!po)
4002 : {
4003 187 : set_avma((pari_sp)g);
4004 187 : g = gerepileuptoleaf(av0, g);
4005 : }
4006 : else {
4007 6986 : *po = mkvec2(subiu(powuu(p,f), 1), o);
4008 6986 : gerepileall(av0, 2, &g, po);
4009 : }
4010 7173 : return g;
4011 : }
4012 :
4013 : static GEN
4014 366572 : _Flxq_neg(void *E, GEN x)
4015 366572 : { struct _Flxq *s = (struct _Flxq *)E;
4016 366572 : return Flx_neg(x,s->p); }
4017 :
4018 : static GEN
4019 1462476 : _Flxq_rmul(void *E, GEN x, GEN y)
4020 1462476 : { struct _Flxq *s = (struct _Flxq *)E;
4021 1462476 : return Flx_mul_pre(x,y,s->p,s->pi); }
4022 :
4023 : static GEN
4024 9460 : _Flxq_inv(void *E, GEN x)
4025 9460 : { struct _Flxq *s = (struct _Flxq *)E;
4026 9460 : return Flxq_inv(x,s->T,s->p); }
4027 :
4028 : static int
4029 69139 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
4030 :
4031 : static GEN
4032 6567 : _Flxq_s(void *E, long x)
4033 6567 : { struct _Flxq *s = (struct _Flxq *)E;
4034 6567 : ulong u = x<0 ? s->p+x: (ulong)x;
4035 6567 : return Fl_to_Flx(u, get_Flx_var(s->T));
4036 : }
4037 :
4038 : static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
4039 : _Flxq_inv,_Flxq_equal0,_Flxq_s};
4040 :
4041 68985 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
4042 : {
4043 68985 : GEN z = new_chunk(sizeof(struct _Flxq));
4044 68985 : set_Flxq((struct _Flxq *)z, T, p); *E = (void*)z; return &Flxq_field;
4045 : }
4046 :
4047 : /***********************************************************************/
4048 : /** Flxn **/
4049 : /***********************************************************************/
4050 :
4051 : GEN
4052 54411 : Flx_invLaplace(GEN x, ulong p)
4053 : {
4054 54411 : long i, d = degpol(x);
4055 : ulong t;
4056 : GEN y;
4057 54409 : if (d <= 1) return Flx_copy(x);
4058 54409 : t = Fl_inv(factorial_Fl(d, p), p);
4059 54463 : y = cgetg(d+3, t_VECSMALL);
4060 54412 : y[1] = x[1];
4061 1327255 : for (i=d; i>=2; i--)
4062 : {
4063 1272810 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
4064 1272797 : t = Fl_mul(t, i, p);
4065 : }
4066 54445 : uel(y,3) = uel(x,3);
4067 54445 : uel(y,2) = uel(x,2);
4068 54445 : return y;
4069 : }
4070 :
4071 : GEN
4072 27373 : Flx_Laplace(GEN x, ulong p)
4073 : {
4074 27373 : long i, d = degpol(x);
4075 27370 : ulong t = 1;
4076 : GEN y;
4077 27370 : if (d <= 1) return Flx_copy(x);
4078 27370 : y = cgetg(d+3, t_VECSMALL);
4079 27351 : y[1] = x[1];
4080 27351 : uel(y,2) = uel(x,2);
4081 27351 : uel(y,3) = uel(x,3);
4082 755871 : for (i=2; i<=d; i++)
4083 : {
4084 728489 : t = Fl_mul(t, i%p, p);
4085 728504 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
4086 : }
4087 27382 : return y;
4088 : }
4089 :
4090 : GEN
4091 6232185 : Flxn_red(GEN a, long n)
4092 : {
4093 6232185 : long i, L, l = lg(a);
4094 : GEN b;
4095 6232185 : if (l == 2 || !n) return zero_Flx(a[1]);
4096 5842264 : L = n+2; if (L > l) L = l;
4097 5842264 : b = cgetg(L, t_VECSMALL); b[1] = a[1];
4098 58584248 : for (i=2; i<L; i++) b[i] = a[i];
4099 5840065 : return Flx_renormalize(b,L);
4100 : }
4101 :
4102 : GEN
4103 5064325 : Flxn_mul_pre(GEN a, GEN b, long n, ulong p, ulong pi)
4104 5064325 : { return Flxn_red(Flx_mul_pre(a, b, p, pi), n); }
4105 : GEN
4106 75391 : Flxn_mul(GEN a, GEN b, long n, ulong p)
4107 75391 : { return Flxn_mul_pre(a, b, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4108 :
4109 : GEN
4110 0 : Flxn_sqr_pre(GEN a, long n, ulong p, ulong pi)
4111 0 : { return Flxn_red(Flx_sqr_pre(a, p, pi), n); }
4112 : GEN
4113 0 : Flxn_sqr(GEN a, long n, ulong p)
4114 0 : { return Flxn_sqr_pre(a, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4115 :
4116 : /* (f*g) \/ x^n */
4117 : static GEN
4118 938603 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p, ulong pi)
4119 938603 : { return Flx_shift(Flx_mul_pre(f, g, p, pi),-n); }
4120 :
4121 : static GEN
4122 516770 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p, ulong pi)
4123 : {
4124 516770 : GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
4125 516503 : return Flx_add(Flx_mulhigh_i(fl, g, n2, p, pi),
4126 : Flxn_mul_pre(fh, g, n - n2, p, pi), p);
4127 : }
4128 :
4129 : /* g==NULL -> assume g==1 */
4130 : GEN
4131 55217 : Flxn_div_pre(GEN g, GEN f, long e, ulong p, ulong pi)
4132 : {
4133 55217 : pari_sp av = avma, av2;
4134 : ulong mask;
4135 : GEN W;
4136 55217 : long n = 1;
4137 55217 : if (lg(f) <= 2) pari_err_INV("Flxn_inv",f);
4138 55217 : W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
4139 55261 : mask = quadratic_prec_mask(e);
4140 55253 : av2 = avma;
4141 259063 : for (;mask>1;)
4142 : {
4143 : GEN u, fr;
4144 203804 : long n2 = n;
4145 203804 : n<<=1; if (mask & 1) n--;
4146 203804 : mask >>= 1;
4147 203804 : fr = Flxn_red(f, n);
4148 203640 : if (mask>1 || !g)
4149 : {
4150 149474 : u = Flxn_mul_pre(W, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
4151 149885 : W = Flx_sub(W, Flx_shift(u, n2), p);
4152 : } else
4153 : {
4154 54166 : GEN y = Flxn_mul_pre(g, W, n, p, pi), yt = Flxn_red(y, n-n2);
4155 54135 : u = Flxn_mul_pre(yt, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
4156 54163 : W = Flx_sub(y, Flx_shift(u, n2), p);
4157 : }
4158 203852 : if (gc_needed(av2,2))
4159 : {
4160 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_div, e = %ld", n);
4161 0 : W = gerepileupto(av2, W);
4162 : }
4163 : }
4164 55259 : return gerepileupto(av, W);
4165 : }
4166 : GEN
4167 55194 : Flxn_div(GEN g, GEN f, long e, ulong p)
4168 55194 : { return Flxn_div_pre(g, f, e, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4169 :
4170 : GEN
4171 1030 : Flxn_inv(GEN f, long e, ulong p)
4172 1030 : { return Flxn_div(NULL, f, e, p); }
4173 :
4174 : GEN
4175 109426 : Flxn_expint(GEN h, long e, ulong p)
4176 : {
4177 109426 : pari_sp av = avma, av2;
4178 109426 : long v = h[1], n=1;
4179 109426 : GEN f = pol1_Flx(v), g = pol1_Flx(v);
4180 109394 : ulong mask = quadratic_prec_mask(e), pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4181 109400 : av2 = avma;
4182 422870 : for (;mask>1;)
4183 : {
4184 : GEN u, w;
4185 422797 : long n2 = n;
4186 422797 : n<<=1; if (mask & 1) n--;
4187 422797 : mask >>= 1;
4188 422797 : u = Flxn_mul_pre(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p,pi), n-n2, p,pi);
4189 422797 : u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
4190 422838 : w = Flxn_mul_pre(f, Flx_integXn(u, n2-1, p), n-n2, p, pi);
4191 422768 : f = Flx_add(f, Flx_shift(w, n2), p);
4192 422918 : if (mask<=1) break;
4193 313491 : u = Flxn_mul_pre(g, Flxn_mulhigh(f, g, n2, n, p, pi), n-n2, p, pi);
4194 313476 : g = Flx_sub(g, Flx_shift(u, n2), p);
4195 313470 : if (gc_needed(av2,2))
4196 : {
4197 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
4198 0 : gerepileall(av2, 2, &f, &g);
4199 : }
4200 : }
4201 109500 : return gerepileupto(av, f);
4202 : }
4203 :
4204 : GEN
4205 0 : Flxn_exp(GEN h, long e, ulong p)
4206 : {
4207 0 : if (degpol(h)<1 || uel(h,2)!=0)
4208 0 : pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
4209 0 : return Flxn_expint(Flx_deriv(h, p), e, p);
4210 : }
4211 :
4212 : INLINE GEN
4213 217447 : Flxn_recip(GEN x, long n)
4214 : {
4215 217447 : GEN z=Flx_recipspec(x+2,lgpol(x),n);
4216 217249 : z[1]=x[1];
4217 217249 : return z;
4218 : }
4219 :
4220 : GEN
4221 54172 : Flx_Newton(GEN P, long n, ulong p)
4222 : {
4223 54172 : pari_sp av = avma;
4224 54172 : long d = degpol(P);
4225 54160 : GEN dP = Flxn_recip(Flx_deriv(P, p), d);
4226 54051 : GEN Q = Flxn_div(dP, Flxn_recip(P, d+1), n, p);
4227 54112 : return gerepileuptoleaf(av, Q);
4228 : }
4229 :
4230 : GEN
4231 109431 : Flx_fromNewton(GEN P, ulong p)
4232 : {
4233 109431 : pari_sp av = avma;
4234 109431 : ulong n = Flx_constant(P)+1;
4235 109430 : GEN z = Flx_neg(Flx_shift(P, -1), p);
4236 109426 : GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
4237 109417 : return gerepileuptoleaf(av, Q);
4238 : }
4239 :
4240 : static void
4241 12514 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
4242 : {
4243 : long i;
4244 : ulong e;
4245 12514 : GEN P = cgetg(d+1, t_VECSMALL);
4246 12514 : GEN V = cgetg(d+1, t_VECSMALL);
4247 1396581 : for (i=1, e=1; i<=d; i++, e++)
4248 : {
4249 1384067 : if (e==p)
4250 : {
4251 459153 : e = 0;
4252 459153 : V[i] = u_lvalrem(i, p, &uel(P,i));
4253 : } else
4254 : {
4255 924914 : V[i] = 0; uel(P,i) = i;
4256 : }
4257 : }
4258 12514 : *pt_P = P; *pt_V = V;
4259 12514 : }
4260 :
4261 : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
4262 : * val large enough to compensate for the power of p in the factorials */
4263 :
4264 : static GEN
4265 497 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long sv)
4266 : {
4267 497 : pari_sp av = avma;
4268 497 : long i, d = n-1, w;
4269 : GEN y, W, E, t;
4270 497 : init_invlaplace(d, p, &E, &W);
4271 497 : t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
4272 497 : w = zv_sum(W);
4273 497 : if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
4274 497 : y = cgetg(d+3,t_POL);
4275 497 : y[1] = evalsigne(1) | sv;
4276 28882 : for (i=d; i>=1; i--)
4277 : {
4278 28385 : gel(y,i+2) = t;
4279 28385 : t = Fp_mulu(t, uel(E,i), q);
4280 28385 : if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
4281 : }
4282 497 : gel(y,2) = t;
4283 497 : return gerepilecopy(av, ZX_renormalize(y, d+3));
4284 : }
4285 :
4286 : GEN
4287 27583 : Flx_composedsum(GEN P, GEN Q, ulong p)
4288 : {
4289 27583 : pari_sp av = avma;
4290 27583 : long n = 1 + degpol(P)*degpol(Q);
4291 27580 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4292 27581 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4293 : GEN R;
4294 27588 : if (p >= (ulong)n)
4295 : {
4296 27091 : GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
4297 27084 : GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
4298 27077 : GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
4299 27087 : R = Flx_fromNewton(L, p);
4300 : } else
4301 : {
4302 497 : long v = factorial_lval(n-1, p);
4303 497 : long w = 1 + ulogint(n-1, p);
4304 497 : GEN pv = powuu(p, v);
4305 497 : GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
4306 497 : GEN iL = ZpX_invLaplace_init(n, q, p, v, P[1]);
4307 497 : GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
4308 497 : GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
4309 497 : GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
4310 497 : GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
4311 497 : R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
4312 : }
4313 27575 : return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
4314 : }
4315 :
4316 : static GEN
4317 3910 : _Flx_composedsum(void *E, GEN a, GEN b)
4318 3910 : { return Flx_composedsum(a, b, (ulong)E); }
4319 :
4320 : GEN
4321 28985 : FlxV_composedsum(GEN V, ulong p)
4322 28985 : { return gen_product(V, (void *)p, &_Flx_composedsum); }
4323 :
4324 : GEN
4325 0 : Flx_composedprod(GEN P, GEN Q, ulong p)
4326 : {
4327 0 : pari_sp av = avma;
4328 0 : long n = 1+ degpol(P)*degpol(Q);
4329 0 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
4330 0 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
4331 : GEN R;
4332 0 : if (p >= (ulong)n)
4333 : {
4334 0 : GEN L = Flx_convol(Flx_Newton(P,n,p), Flx_Newton(Q,n,p), p);
4335 0 : R = Flx_fromNewton(L, p);
4336 : } else
4337 : {
4338 0 : long w = 1 + ulogint(n, p);
4339 0 : GEN qf = powuu(p, w);
4340 0 : GEN Pl = FpX_convol(FpX_Newton(Flx_to_ZX(P), n, qf), FpX_Newton(Flx_to_ZX(Q), n, qf), qf);
4341 0 : R = ZX_to_Flx(FpX_fromNewton(Pl, qf), p);
4342 : }
4343 0 : return gerepileuptoleaf(av, Flx_Fl_mul(R, lead, p));
4344 :
4345 : }
4346 :
4347 : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
4348 : static GEN
4349 0 : Fl_Xp1_powu(ulong n, ulong p, long v)
4350 : {
4351 0 : ulong k, d = (n + 1) >> 1;
4352 0 : GEN C, V = identity_zv(d);
4353 :
4354 0 : Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
4355 0 : C = cgetg(n+3, t_VECSMALL);
4356 0 : C[1] = v;
4357 0 : uel(C,2) = 1UL;
4358 0 : uel(C,3) = n%p;
4359 0 : uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
4360 : /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
4361 0 : if (SMALL_ULONG(p))
4362 0 : for (k = 3; k <= d; k++)
4363 0 : uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
4364 : else
4365 : {
4366 0 : ulong pi = get_Fl_red(p);
4367 0 : for (k = 3; k <= d; k++)
4368 0 : uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
4369 : }
4370 0 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4371 0 : return C; /* normalized */
4372 : }
4373 :
4374 : /* p arbitrary */
4375 : GEN
4376 28236 : Flx_translate1_basecase(GEN P, ulong p)
4377 : {
4378 28236 : GEN R = Flx_copy(P);
4379 28236 : long i, k, n = degpol(P);
4380 654893 : for (i = 1; i <= n; i++)
4381 14846873 : for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
4382 28236 : return R;
4383 : }
4384 :
4385 : static int
4386 41401 : translate_basecase(long n, ulong p)
4387 : {
4388 : #ifdef LONG_IS_64BIT
4389 36102 : if (p <= 19) return n < 40;
4390 29910 : if (p < 1UL<<30) return n < 58;
4391 0 : if (p < 1UL<<59) return n < 100;
4392 0 : if (p < 1UL<<62) return n < 120;
4393 0 : if (p < 1UL<<63) return n < 240;
4394 0 : return n < 250;
4395 : #else
4396 5299 : if (p <= 13) return n < 18;
4397 4136 : if (p <= 17) return n < 22;
4398 4078 : if (p <= 29) return n < 39;
4399 3886 : if (p <= 67) return n < 69;
4400 3667 : if (p < 1UL<< 15) return n < 80;
4401 2047 : if (p < 1UL<< 16) return n < 100;
4402 0 : if (p < 1UL<< 28) return n < 300;
4403 0 : return n < 650;
4404 : #endif
4405 : }
4406 : /* assume p prime */
4407 : GEN
4408 16142 : Flx_translate1(GEN P, ulong p)
4409 : {
4410 16142 : long d, n = degpol(P);
4411 : GEN R, Q, S;
4412 16142 : if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
4413 : /* n > 0 */
4414 1148 : d = n >> 1;
4415 1148 : if ((ulong)n < p)
4416 : {
4417 0 : R = Flx_translate1(Flxn_red(P, d), p);
4418 0 : Q = Flx_translate1(Flx_shift(P, -d), p);
4419 0 : S = Fl_Xp1_powu(d, p, P[1]);
4420 0 : return Flx_add(Flx_mul(Q, S, p), R, p);
4421 : }
4422 : else
4423 : {
4424 : ulong q;
4425 1148 : if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
4426 1148 : R = Flx_translate1(Flxn_red(P, q), p);
4427 1148 : Q = Flx_translate1(Flx_shift(P, -q), p);
4428 1148 : S = Flx_add(Flx_shift(Q, q), Q, p);
4429 1148 : return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
4430 : }
4431 : }
4432 :
4433 : GEN
4434 0 : Flx_translate(GEN P, ulong c, ulong p)
4435 : {
4436 0 : pari_sp av = avma;
4437 : GEN Q;
4438 0 : if (c==0) return Flx_copy(P);
4439 0 : if (c==1) return Flx_translate1(P, p);
4440 0 : Q = Flx_unscale(Flx_translate1(Flx_unscale(P, c, p), p), Fl_inv(c, p), p);
4441 0 : return gerepileuptoleaf(av, Q);
4442 : }
4443 :
4444 : static GEN
4445 12017 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
4446 : {
4447 12017 : ulong k, d = n >> 1, c, v = 0;
4448 12017 : GEN C, V, W, U = upowers(p, e-1);
4449 12017 : init_invlaplace(d, p, &V, &W);
4450 12017 : Flv_inv_inplace(V, q);
4451 12017 : C = cgetg(n+3, t_VECSMALL);
4452 12017 : C[1] = vs;
4453 12017 : uel(C,2) = 1UL;
4454 12017 : uel(C,3) = n%q;
4455 12017 : v = u_lvalrem(n, p, &c);
4456 1355682 : for (k = 2; k <= d; k++)
4457 : {
4458 : ulong w;
4459 1343665 : v += u_lvalrem(n-k+1, p, &w) - W[k];
4460 1343665 : c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
4461 1343665 : uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
4462 : }
4463 1374521 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4464 12017 : return C; /* normalized */
4465 : }
4466 :
4467 : GEN
4468 25259 : zlx_translate1(GEN P, ulong p, long e)
4469 : {
4470 25259 : ulong d, q = upowuu(p,e), n = degpol(P);
4471 : GEN R, Q, S;
4472 25259 : if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
4473 : /* n > 0 */
4474 12017 : d = n >> 1;
4475 12017 : R = zlx_translate1(Flxn_red(P, d), p, e);
4476 12017 : Q = zlx_translate1(Flx_shift(P, -d), p, e);
4477 12017 : S = zl_Xp1_powu(d, p, q, e, P[1]);
4478 12017 : return Flx_add(Flx_mul(Q, S, q), R, q);
4479 : }
4480 :
4481 : /***********************************************************************/
4482 : /** Fl2 **/
4483 : /***********************************************************************/
4484 : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
4485 : * a nonsquare D. */
4486 :
4487 : INLINE GEN
4488 7186642 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
4489 :
4490 : /* allow pi = 0 */
4491 : GEN
4492 1914829 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4493 : {
4494 : ulong xaya, xbyb, Db2, mid, z1, z2;
4495 1914829 : ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
4496 1914829 : if (pi)
4497 : {
4498 1914841 : xaya = Fl_mul_pre(x1,y1,p,pi);
4499 1915371 : if (x2==0 && y2==0) return mkF2(xaya,0);
4500 1846518 : if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
4501 1821795 : if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
4502 1821486 : xbyb = Fl_mul_pre(x2,y2,p,pi);
4503 1821436 : mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
4504 1821637 : Db2 = Fl_mul_pre(D, xbyb, p,pi);
4505 : }
4506 0 : else if (p & HIGHMASK)
4507 : {
4508 0 : xaya = Fl_mul(x1,y1,p);
4509 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4510 0 : if (x2==0) return mkF2(xaya,Fl_mul(x1,y2,p));
4511 0 : if (y2==0) return mkF2(xaya,Fl_mul(x2,y1,p));
4512 0 : xbyb = Fl_mul(x2,y2,p);
4513 0 : mid = Fl_mul(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p);
4514 0 : Db2 = Fl_mul(D, xbyb, p);
4515 : }
4516 : else
4517 : {
4518 0 : xaya = (x1 * y1) % p;
4519 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4520 0 : if (x2==0) return mkF2(xaya, (x1 * y2) % p);
4521 0 : if (y2==0) return mkF2(xaya, (x2 * y1) % p);
4522 0 : xbyb = (x2 * y2) % p;
4523 0 : mid = (Fl_add(x1,x2,p) * Fl_add(y1,y2,p)) % p;
4524 0 : Db2 = (D * xbyb) % p;
4525 : }
4526 1821570 : z1 = Fl_add(xaya,Db2,p);
4527 1821529 : z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
4528 1821533 : return mkF2(z1,z2);
4529 : }
4530 :
4531 : /* allow pi = 0 */
4532 : GEN
4533 4819072 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
4534 : {
4535 4819072 : ulong a = x[1], b = x[2];
4536 : ulong a2, Db2, ab;
4537 4819072 : if (pi)
4538 : {
4539 4819095 : a2 = Fl_sqr_pre(a,p,pi);
4540 4821941 : if (b==0) return mkF2(a2,0);
4541 4609379 : Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
4542 4609390 : ab = Fl_mul_pre(a,b,p,pi);
4543 : }
4544 0 : else if (p & HIGHMASK)
4545 : {
4546 0 : a2 = Fl_sqr(a,p);
4547 0 : if (b==0) return mkF2(a2,0);
4548 0 : Db2= Fl_mul(D, Fl_sqr(b,p), p);
4549 0 : ab = Fl_mul(a,b,p);
4550 : }
4551 : else
4552 : {
4553 0 : a2 = (a * a) % p;
4554 0 : if (b==0) return mkF2(a2,0);
4555 0 : Db2= (D * ((b * b) % p)) % p;
4556 0 : ab = (a * b) % p;
4557 : }
4558 4609397 : return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
4559 : }
4560 :
4561 : /* allow pi = 0 */
4562 : ulong
4563 126016 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
4564 : {
4565 126016 : ulong a = x[1], b = x[2], a2;
4566 126016 : if (pi)
4567 : {
4568 72373 : a2 = Fl_sqr_pre(a,p,pi);
4569 72373 : return b? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p): a2;
4570 : }
4571 53643 : else if (p & HIGHMASK)
4572 : {
4573 0 : a2 = Fl_sqr(a,p);
4574 0 : return b? Fl_sub(a2, Fl_mul(D, Fl_sqr(b, p), p), p): a2;
4575 : }
4576 : else
4577 : {
4578 53643 : a2 = (a * a) % p;
4579 53643 : return b? Fl_sub(a2, (D * ((b * b) % p)) % p, p): a2;
4580 : }
4581 : }
4582 :
4583 : /* allow pi = 0 */
4584 : GEN
4585 192742 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
4586 : {
4587 192742 : ulong a = x[1], b = x[2], n, ni;
4588 192742 : if (b == 0) return mkF2(Fl_inv(a,p), 0);
4589 162040 : b = Fl_neg(b, p);
4590 162043 : if (pi)
4591 : {
4592 162043 : n = Fl_sub(Fl_sqr_pre(a, p,pi),
4593 : Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p);
4594 162042 : ni = Fl_inv(n,p);
4595 162042 : return mkF2(Fl_mul_pre(a, ni, p,pi), Fl_mul_pre(b, ni, p,pi));
4596 : }
4597 0 : else if (p & HIGHMASK)
4598 : {
4599 0 : n = Fl_sub(Fl_sqr(a, p), Fl_mul(D, Fl_sqr(b, p), p), p);
4600 0 : ni = Fl_inv(n,p);
4601 0 : return mkF2(Fl_mul(a, ni, p), Fl_mul(b, ni, p));
4602 : }
4603 : else
4604 : {
4605 0 : n = Fl_sub((a * a) % p, (D * ((b * b) % p)) % p, p);
4606 0 : ni = Fl_inv(n,p);
4607 0 : return mkF2((a * ni) % p, (b * ni) % p);
4608 : }
4609 : }
4610 :
4611 : int
4612 439402 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
4613 :
4614 : struct _Fl2 {
4615 : ulong p, pi, D;
4616 : };
4617 :
4618 : static GEN
4619 4819114 : _Fl2_sqr(void *data, GEN x)
4620 : {
4621 4819114 : struct _Fl2 *D = (struct _Fl2*)data;
4622 4819114 : return Fl2_sqr_pre(x, D->D, D->p, D->pi);
4623 : }
4624 : static GEN
4625 1886390 : _Fl2_mul(void *data, GEN x, GEN y)
4626 : {
4627 1886390 : struct _Fl2 *D = (struct _Fl2*)data;
4628 1886390 : return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
4629 : }
4630 :
4631 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL; allow pi = 0 */
4632 : GEN
4633 656464 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
4634 : {
4635 656464 : pari_sp av = avma;
4636 : struct _Fl2 d;
4637 : GEN y;
4638 656464 : long s = signe(n);
4639 656464 : if (!s) return mkF2(1,0);
4640 581860 : if (s < 0)
4641 192741 : x = Fl2_inv_pre(x,D,p,pi);
4642 581856 : if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
4643 428773 : d.p = p; d.pi = pi; d.D=D;
4644 428773 : y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
4645 428795 : return gerepileuptoleaf(av, y);
4646 : }
4647 :
4648 : static GEN
4649 656455 : _Fl2_pow(void *data, GEN x, GEN n)
4650 : {
4651 656455 : struct _Fl2 *D = (struct _Fl2*)data;
4652 656455 : return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
4653 : }
4654 :
4655 : static GEN
4656 111072 : _Fl2_rand(void *data)
4657 : {
4658 111072 : struct _Fl2 *D = (struct _Fl2*)data;
4659 111072 : ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
4660 111077 : return mkF2(a,b);
4661 : }
4662 :
4663 : GEN
4664 67676 : Fl2_sqrt_pre(GEN z, ulong D, ulong p, ulong pi)
4665 : {
4666 67676 : ulong a = uel(z,1), b = uel(z,2), as2, u, v, s;
4667 67676 : ulong y = Fl_2gener_pre_i(D, p, pi);
4668 67676 : if (b == 0)
4669 19383 : return krouu(a, p)==1 ? mkF2(Fl_sqrt_pre_i(a, y, p, pi), 0)
4670 19383 : : mkF2(0, Fl_sqrt_pre_i(Fl_div(a, D, p), y, p, pi));
4671 54473 : s = Fl_sqrt_pre_i(Fl2_norm_pre(z, D, p, pi), y, p, pi);
4672 54473 : if (s==~0UL) return NULL;
4673 51299 : as2 = Fl_halve(Fl_add(a, s, p), p);
4674 51299 : if (krouu(as2, p)==-1) as2 = Fl_sub(as2, s, p);
4675 51299 : u = Fl_sqrt_pre_i(as2, y, p, pi);
4676 51299 : v = Fl_div(b, Fl_double(u, p), p);
4677 51299 : return mkF2(u,v);
4678 : }
4679 :
4680 : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
4681 : hash_GEN, zv_equal, Fl2_equal1, NULL};
4682 :
4683 : /* allow pi = 0 */
4684 : GEN
4685 74603 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
4686 : {
4687 : struct _Fl2 E;
4688 : GEN o;
4689 74603 : if (a[1]==0 && a[2]==0)
4690 : {
4691 0 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
4692 0 : if (zeta) *zeta=mkF2(1,0);
4693 0 : return zv_copy(a);
4694 : }
4695 74603 : E.p=p; E.pi = pi; E.D = D;
4696 74603 : o = subiu(powuu(p,2), 1);
4697 74600 : return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
4698 : }
4699 :
4700 : /* allow pi = 0 */
4701 : GEN
4702 10528 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4703 : {
4704 : GEN p1;
4705 10528 : long i = lg(x)-1;
4706 10528 : if (i <= 2)
4707 2086 : return mkF2(i == 2? x[2]: 0, 0);
4708 8442 : p1 = mkF2(x[i], 0);
4709 36876 : for (i--; i>=2; i--)
4710 : {
4711 28434 : p1 = Fl2_mul_pre(p1, y, D, p, pi);
4712 28434 : uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
4713 : }
4714 8442 : return p1;
4715 : }
4716 :
4717 : /***********************************************************************/
4718 : /** FlxV **/
4719 : /***********************************************************************/
4720 : /* FlxV are t_VEC with Flx coefficients. */
4721 :
4722 : GEN
4723 34482 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
4724 : {
4725 34482 : pari_sp ltop=avma;
4726 : long i;
4727 34482 : GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
4728 257068 : for(i=2;i<lg(V);i++)
4729 222586 : z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
4730 34482 : return gerepileuptoleaf(ltop,z);
4731 : }
4732 :
4733 : GEN
4734 0 : ZXV_to_FlxV(GEN x, ulong p)
4735 0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
4736 :
4737 : GEN
4738 3798927 : ZXT_to_FlxT(GEN x, ulong p)
4739 : {
4740 3798927 : if (typ(x) == t_POL)
4741 3740443 : return ZX_to_Flx(x, p);
4742 : else
4743 192032 : pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
4744 : }
4745 :
4746 : GEN
4747 171823 : FlxV_to_Flm(GEN x, long n)
4748 927200 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
4749 :
4750 : GEN
4751 0 : FlxV_red(GEN x, ulong p)
4752 0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4753 :
4754 : GEN
4755 292093 : FlxT_red(GEN x, ulong p)
4756 : {
4757 292093 : if (typ(x) == t_VECSMALL)
4758 196544 : return Flx_red(x, p);
4759 : else
4760 320408 : pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4761 : }
4762 :
4763 : GEN
4764 113589 : FlxqV_dotproduct_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
4765 : {
4766 113589 : long i, lx = lg(x);
4767 : pari_sp av;
4768 : GEN c;
4769 113589 : if (lx == 1) return pol0_Flx(get_Flx_var(T));
4770 113589 : av = avma; c = Flx_mul_pre(gel(x,1),gel(y,1), p, pi);
4771 464499 : for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4772 113589 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4773 : }
4774 : GEN
4775 0 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4776 0 : { return FlxqV_dotproduct_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4777 :
4778 : GEN
4779 2000 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4780 : {
4781 2000 : long i, l = minss(lg(x), lg(y));
4782 : ulong pi;
4783 : pari_sp av;
4784 : GEN c;
4785 2000 : if (l == 2) return pol0_Flx(get_Flx_var(T));
4786 1987 : av = avma; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4787 1987 : c = Flx_mul_pre(gel(x,2),gel(y,2), p, pi);
4788 6576 : for (i=3; i<l; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4789 1987 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4790 : }
4791 :
4792 : /* allow pi = 0 */
4793 : GEN
4794 251067 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4795 : {
4796 251067 : long i, l = lg(z);
4797 251067 : GEN y = cgetg(l, t_VECSMALL);
4798 12737351 : for (i=1; i<l; i++) uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4799 251207 : return y;
4800 : }
4801 :
4802 : /***********************************************************************/
4803 : /** FlxM **/
4804 : /***********************************************************************/
4805 : /* allow pi = 0 */
4806 : GEN
4807 19452 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4808 : {
4809 19452 : long i, l = lg(z);
4810 19452 : GEN y = cgetg(l, t_MAT);
4811 270520 : for (i=1; i<l; i++) gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4812 19453 : return y;
4813 : }
4814 :
4815 : GEN
4816 0 : zero_FlxC(long n, long sv)
4817 : {
4818 0 : GEN x = cgetg(n + 1, t_COL), z = zero_Flx(sv);
4819 : long i;
4820 0 : for (i = 1; i <= n; i++) gel(x, i) = z;
4821 0 : return x;
4822 : }
4823 :
4824 : GEN
4825 0 : FlxC_neg(GEN x, ulong p)
4826 0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4827 :
4828 : GEN
4829 0 : FlxC_sub(GEN x, GEN y, ulong p)
4830 0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4831 :
4832 : GEN
4833 0 : zero_FlxM(long r, long c, long sv)
4834 : {
4835 0 : GEN x = cgetg(c + 1, t_MAT), z = zero_FlxC(r, sv);
4836 : long j;
4837 0 : for (j = 1; j <= c; j++) gel(x, j) = z;
4838 0 : return x;
4839 : }
4840 :
4841 : GEN
4842 0 : FlxM_neg(GEN x, ulong p)
4843 0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4844 :
4845 : GEN
4846 0 : FlxM_sub(GEN x, GEN y, ulong p)
4847 0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4848 :
4849 : GEN
4850 0 : FlxC_translate(GEN x, ulong c, ulong p)
4851 0 : { pari_APPLY_type(t_COL, Flx_translate(gel(x,i), c, p)) }
4852 :
4853 : GEN
4854 0 : FlxM_translate(GEN x, ulong c, ulong p)
4855 0 : { pari_APPLY_same(FlxC_translate(gel(x,i), c, p)) }
4856 :
4857 : GEN
4858 235087 : FlxqC_red_pre(GEN x, GEN T, ulong p, ulong pi)
4859 4061629 : { pari_APPLY_type(t_COL, Flx_rem_pre(gel(x,i), T, p, pi)) }
4860 :
4861 : GEN
4862 81833 : FlxqM_red_pre(GEN x, GEN T, ulong p, ulong pi)
4863 316920 : { pari_APPLY_same(FlxqC_red_pre(gel(x,i), T, p, pi)) }
4864 :
4865 : GEN
4866 0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4867 0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4868 :
4869 : GEN
4870 0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4871 0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4872 :
4873 : static GEN
4874 47231 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4875 : long i, j, l, lc;
4876 47231 : GEN N = cgetg_copy(M, &l), x;
4877 47231 : if (l == 1)
4878 0 : return N;
4879 47231 : lc = lgcols(M);
4880 206197 : for (j = 1; j < l; j++) {
4881 158966 : gel(N, j) = cgetg(lc, t_COL);
4882 906164 : for (i = 1; i < lc; i++) {
4883 747198 : x = gcoeff(M, i, j);
4884 747198 : gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4885 : }
4886 : }
4887 47231 : return N;
4888 : }
4889 :
4890 : static GEN
4891 690005 : kron_pack_Flx_spec_half(GEN x, long l) {
4892 690005 : if (l == 0) return gen_0;
4893 459191 : return Flx_to_int_halfspec(x, l);
4894 : }
4895 :
4896 : static GEN
4897 53804 : kron_pack_Flx_spec(GEN x, long l) {
4898 : long i;
4899 : GEN w, y;
4900 53804 : if (l == 0)
4901 10072 : return gen_0;
4902 43732 : y = cgetipos(l + 2);
4903 159397 : for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4904 115665 : *w = x[i];
4905 43732 : return y;
4906 : }
4907 :
4908 : static GEN
4909 3389 : kron_pack_Flx_spec_2(GEN x, long l) { return Flx_eval2BILspec(x, 2, l); }
4910 :
4911 : static GEN
4912 0 : kron_pack_Flx_spec_3(GEN x, long l) { return Flx_eval2BILspec(x, 3, l); }
4913 :
4914 : static GEN
4915 42969 : kron_unpack_Flx(GEN z, ulong p)
4916 : {
4917 42969 : long i, l = lgefint(z);
4918 42969 : GEN x = cgetg(l, t_VECSMALL), w;
4919 202007 : for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4920 159038 : x[i] = ((ulong) *w) % p;
4921 42969 : return Flx_renormalize(x, l);
4922 : }
4923 :
4924 : static GEN
4925 2930 : kron_unpack_Flx_2(GEN x, ulong p) {
4926 2930 : long d = (lgefint(x)-1)/2 - 1;
4927 2930 : return Z_mod2BIL_Flx_2(x, d, p);
4928 : }
4929 :
4930 : static GEN
4931 0 : kron_unpack_Flx_3(GEN x, ulong p) {
4932 0 : long d = lgefint(x)/3 - 1;
4933 0 : return Z_mod2BIL_Flx_3(x, d, p);
4934 : }
4935 :
4936 : static GEN
4937 116347 : FlxM_pack_ZM_bits(GEN M, long b)
4938 : {
4939 : long i, j, l, lc;
4940 116347 : GEN N = cgetg_copy(M, &l), x;
4941 116347 : if (l == 1)
4942 0 : return N;
4943 116347 : lc = lgcols(M);
4944 479970 : for (j = 1; j < l; j++) {
4945 363623 : gel(N, j) = cgetg(lc, t_COL);
4946 5955976 : for (i = 1; i < lc; i++) {
4947 5592353 : x = gcoeff(M, i, j);
4948 5592353 : gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
4949 : }
4950 : }
4951 116347 : return N;
4952 : }
4953 :
4954 : static GEN
4955 23619 : ZM_unpack_FlxM(GEN M, ulong p, ulong sv, GEN (*unpack)(GEN, ulong))
4956 : {
4957 : long i, j, l, lc;
4958 23619 : GEN N = cgetg_copy(M, &l), x;
4959 23619 : if (l == 1)
4960 0 : return N;
4961 23619 : lc = lgcols(M);
4962 111622 : for (j = 1; j < l; j++) {
4963 88003 : gel(N, j) = cgetg(lc, t_COL);
4964 635759 : for (i = 1; i < lc; i++) {
4965 547756 : x = unpack(gcoeff(M, i, j), p);
4966 547756 : x[1] = sv;
4967 547756 : gcoeff(N, i, j) = x;
4968 : }
4969 : }
4970 23619 : return N;
4971 : }
4972 :
4973 : static GEN
4974 58214 : ZM_unpack_FlxM_bits(GEN M, long b, ulong p, ulong pi, long sv)
4975 : {
4976 : long i, j, l, lc;
4977 58214 : GEN N = cgetg_copy(M, &l), x;
4978 58214 : if (l == 1)
4979 0 : return N;
4980 58214 : lc = lgcols(M);
4981 58214 : if (b < BITS_IN_LONG) {
4982 195466 : for (j = 1; j < l; j++) {
4983 138929 : gel(N, j) = cgetg(lc, t_COL);
4984 3250599 : for (i = 1; i < lc; i++) {
4985 3111670 : x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
4986 3111670 : x[1] = sv;
4987 3111670 : gcoeff(N, i, j) = x;
4988 : }
4989 : }
4990 : } else {
4991 1677 : if (!pi) pi = get_Fl_red(p); /* unset if !SMALL_ULONG(p) */
4992 9832 : for (j = 1; j < l; j++) {
4993 8155 : gel(N, j) = cgetg(lc, t_COL);
4994 175271 : for (i = 1; i < lc; i++) {
4995 167116 : x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
4996 167116 : x[1] = sv;
4997 167116 : gcoeff(N, i, j) = x;
4998 : }
4999 : }
5000 : }
5001 58214 : return N;
5002 : }
5003 :
5004 : static GEN
5005 81833 : FlxM_mul_Kronecker_i(GEN A, GEN B, ulong p, ulong pi, long d, long sv)
5006 : {
5007 81833 : long b, n = lg(A) - 1;
5008 : GEN C, z;
5009 : GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
5010 81833 : int is_sqr = A==B;
5011 :
5012 81833 : z = muliu(muliu(sqru(p - 1), d), n);
5013 81833 : b = expi(z) + 1;
5014 : /* only do expensive bit-packing if it saves at least 1 limb */
5015 81833 : if (b <= BITS_IN_HALFULONG)
5016 77428 : { if (nbits2nlong(d*b) == (d + 1)/2) b = BITS_IN_HALFULONG; }
5017 : else
5018 : {
5019 4405 : long l = lgefint(z) - 2;
5020 4405 : if (nbits2nlong(d*b) == d*l) b = l*BITS_IN_LONG;
5021 : }
5022 :
5023 81833 : switch (b) {
5024 22554 : case BITS_IN_HALFULONG:
5025 22554 : pack = kron_pack_Flx_spec_half;
5026 22554 : unpack = int_to_Flx_half;
5027 22554 : break;
5028 1016 : case BITS_IN_LONG:
5029 1016 : pack = kron_pack_Flx_spec;
5030 1016 : unpack = kron_unpack_Flx;
5031 1016 : break;
5032 49 : case 2*BITS_IN_LONG:
5033 49 : pack = kron_pack_Flx_spec_2;
5034 49 : unpack = kron_unpack_Flx_2;
5035 49 : break;
5036 0 : case 3*BITS_IN_LONG:
5037 0 : pack = kron_pack_Flx_spec_3;
5038 0 : unpack = kron_unpack_Flx_3;
5039 0 : break;
5040 58214 : default:
5041 58214 : A = FlxM_pack_ZM_bits(A, b);
5042 58214 : B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
5043 58214 : C = ZM_mul(A, B);
5044 58214 : return ZM_unpack_FlxM_bits(C, b, p, pi, sv);
5045 : }
5046 23619 : A = FlxM_pack_ZM(A, pack);
5047 23619 : B = is_sqr? A: FlxM_pack_ZM(B, pack);
5048 23619 : C = ZM_mul(A, B);
5049 23619 : return ZM_unpack_FlxM(C, p, sv, unpack);
5050 : }
5051 :
5052 : GEN
5053 81833 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
5054 : {
5055 81833 : pari_sp av = avma;
5056 81833 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
5057 81833 : long sv = get_Flx_var(T), d = get_Flx_degree(T);
5058 81833 : GEN C = FlxM_mul_Kronecker_i(A, B, p, pi, d, sv);
5059 81833 : C = FlxqM_red_pre(C, T, p, pi);
5060 81833 : return gerepileupto(av, C);
5061 : }
5062 :
5063 : /* assume m > 1 */
5064 : static long
5065 0 : FlxV_max_degree_i(GEN x, long m)
5066 : {
5067 0 : long i, l = degpol(gel(x,1));
5068 0 : for (i = 2; i < m; i++) l = maxss(l, degpol(gel(x,i)));
5069 0 : return l;
5070 : }
5071 :
5072 : /* assume n > 1 and m > 1 */
5073 : static long
5074 0 : FlxM_max_degree_i(GEN x, long n, long m)
5075 : {
5076 0 : long j, l = FlxV_max_degree_i(gel(x,1), m);
5077 0 : for (j = 2; j < n; j++) l = maxss(l, FlxV_max_degree_i(gel(x,j), m));
5078 0 : return l;
5079 : }
5080 :
5081 : static long
5082 0 : FlxM_max_degree(GEN x)
5083 : {
5084 0 : long n = lg(x), m;
5085 0 : if (n == 1) return -1;
5086 0 : m = lgcols(x); return m == 1? -1: FlxM_max_degree_i(x, n, m);
5087 : }
5088 :
5089 : GEN
5090 0 : FlxM_mul(GEN x, GEN y, ulong p)
5091 : {
5092 0 : pari_sp av = avma;
5093 0 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
5094 : long sv, d;
5095 0 : if (lg(x) == 1) return cgetg(1,t_MAT);
5096 0 : if (lg(gel(x,1))==1) return FlxqM_mul(x, y, NULL, p);
5097 0 : sv = mael3(x,1,1,1);
5098 0 : d = maxss(FlxM_max_degree(x), FlxM_max_degree(y));
5099 0 : return gerepilecopy(av, FlxM_mul_Kronecker_i(x, y, p, pi, d+1, sv));
5100 : }
|
