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. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : /* Not so fast arithmetic with polynomials with small coefficients. */
18 :
19 : static GEN
20 528298768 : get_Flx_red(GEN T, GEN *B)
21 : {
22 528298768 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
23 455936 : *B = gel(T,1); return gel(T,2);
24 : }
25 :
26 : /***********************************************************************/
27 : /** **/
28 : /** Flx **/
29 : /** **/
30 : /***********************************************************************/
31 : /* Flx objects are defined as follows:
32 : Let l an ulong. An Flx is a t_VECSMALL:
33 : x[0] = codeword
34 : x[1] = evalvarn(variable number) (signe is not stored).
35 : x[2] = a_0 x[3] = a_1, etc.
36 : With 0 <= a_i < l
37 :
38 : signe(x) is not valid. Use degpol(x)>=0 instead.
39 : */
40 : /***********************************************************************/
41 : /** **/
42 : /** Conversion from Flx **/
43 : /** **/
44 : /***********************************************************************/
45 :
46 : GEN
47 17261710 : Flx_to_ZX(GEN z)
48 : {
49 17261710 : long i, l = lg(z);
50 17261710 : GEN x = cgetg(l,t_POL);
51 135909801 : for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
52 17261345 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
53 : }
54 :
55 : GEN
56 30499 : Flx_to_FlxX(GEN z, long sv)
57 : {
58 30499 : long i, l = lg(z);
59 30499 : GEN x = cgetg(l,t_POL);
60 108423 : for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
61 30499 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
62 : }
63 :
64 : /* same as Flx_to_ZX, in place */
65 : GEN
66 16124711 : Flx_to_ZX_inplace(GEN z)
67 : {
68 16124711 : long i, l = lg(z);
69 130134731 : for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
70 16123680 : settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
71 : }
72 :
73 : /*Flx_to_Flv=zx_to_zv*/
74 : GEN
75 24235055 : Flx_to_Flv(GEN x, long N)
76 : {
77 : long i, l;
78 24235055 : GEN z = cgetg(N+1,t_VECSMALL);
79 24231497 : if (typ(x) != t_VECSMALL) pari_err_TYPE("Flx_to_Flv",x);
80 24231550 : l = lg(x)-1; x++;
81 354081571 : for (i=1; i<l ; i++) z[i]=x[i];
82 178345992 : for ( ; i<=N; i++) z[i]=0;
83 24231550 : return z;
84 : }
85 :
86 : /*Flv_to_Flx=zv_to_zx*/
87 : GEN
88 10929418 : Flv_to_Flx(GEN x, long sv)
89 : {
90 10929418 : long i, l=lg(x)+1;
91 10929418 : GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
92 10926025 : x--;
93 139637268 : for (i=2; i<l ; i++) z[i]=x[i];
94 10926025 : return Flx_renormalize(z,l);
95 : }
96 :
97 : /*Flm_to_FlxV=zm_to_zxV*/
98 : GEN
99 2226 : Flm_to_FlxV(GEN x, long sv)
100 6027 : { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
101 :
102 : /*FlxC_to_ZXC=zxC_to_ZXC*/
103 : GEN
104 26932 : FlxC_to_ZXC(GEN x)
105 115006 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
106 :
107 : /*FlxC_to_ZXC=zxV_to_ZXV*/
108 : GEN
109 166160 : FlxV_to_ZXV(GEN x)
110 940888 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
111 :
112 : void
113 849272 : FlxV_to_ZXV_inplace(GEN v)
114 : {
115 : long i;
116 2380851 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
117 849203 : }
118 :
119 : /*FlxM_to_ZXM=zxM_to_ZXM*/
120 : GEN
121 4130 : FlxM_to_ZXM(GEN x)
122 14311 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
123 :
124 : GEN
125 341460 : FlxV_to_FlxX(GEN x, long v)
126 : {
127 341460 : long i, l = lg(x)+1;
128 341460 : GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
129 341460 : x--;
130 4210772 : for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
131 341460 : return FlxX_renormalize(z,l);
132 : }
133 :
134 : GEN
135 0 : FlxM_Flx_add_shallow(GEN x, GEN y, ulong p)
136 : {
137 0 : long l = lg(x), i, j;
138 0 : GEN z = cgetg(l,t_MAT);
139 :
140 0 : if (l==1) return z;
141 0 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
142 0 : for (i=1; i<l; i++)
143 : {
144 0 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
145 0 : gel(z,i) = zi;
146 0 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
147 0 : gel(zi,i) = Flx_add(gel(zi,i), y, p);
148 : }
149 0 : return z;
150 : }
151 :
152 : /***********************************************************************/
153 : /** **/
154 : /** Conversion to Flx **/
155 : /** **/
156 : /***********************************************************************/
157 : /* Take an integer and return a scalar polynomial mod p, with evalvarn=vs */
158 : GEN
159 9752975 : Fl_to_Flx(ulong x, long sv)
160 : {
161 9752975 : return x? mkvecsmall2(sv, x): pol0_Flx(sv);
162 : }
163 :
164 : /* a X^d */
165 : GEN
166 283645 : monomial_Flx(ulong a, long d, long vs)
167 : {
168 : GEN P;
169 283645 : if (a==0) return pol0_Flx(vs);
170 283645 : P = const_vecsmall(d+2, 0);
171 283661 : P[1] = vs; P[d+2] = a;
172 283661 : return P;
173 : }
174 :
175 : GEN
176 1024499 : Z_to_Flx(GEN x, ulong p, long sv)
177 : {
178 1024499 : long u = umodiu(x,p);
179 1024500 : return u? mkvecsmall2(sv, u): pol0_Flx(sv);
180 : }
181 :
182 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
183 : GEN
184 73805957 : ZX_to_Flx(GEN x, ulong p)
185 : {
186 73805957 : long i, lx = lg(x);
187 73805957 : GEN a = cgetg(lx, t_VECSMALL);
188 73803708 : a[1]=((ulong)x[1])&VARNBITS;
189 547746501 : for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
190 73804514 : return Flx_renormalize(a,lx);
191 : }
192 :
193 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
194 : GEN
195 3741340 : zx_to_Flx(GEN x, ulong p)
196 : {
197 3741340 : long i, lx = lg(x);
198 3741340 : GEN a = cgetg(lx, t_VECSMALL);
199 3741340 : a[1] = x[1];
200 11153737 : for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
201 3741340 : return Flx_renormalize(a,lx);
202 : }
203 :
204 : ulong
205 48174947 : Rg_to_Fl(GEN x, ulong p)
206 : {
207 48174947 : switch(typ(x))
208 : {
209 24573958 : case t_INT: return umodiu(x, p);
210 61017 : case t_FRAC: {
211 61017 : ulong z = umodiu(gel(x,1), p);
212 61017 : if (!z) return 0;
213 57853 : return Fl_div(z, umodiu(gel(x,2), p), p);
214 : }
215 49 : case t_PADIC: return padic_to_Fl(x, p);
216 23539941 : case t_INTMOD: {
217 23539941 : GEN q = gel(x,1), a = gel(x,2);
218 23539941 : if (absequaliu(q, p)) return itou(a);
219 0 : if (!dvdiu(q,p)) pari_err_MODULUS("Rg_to_Fl", q, utoi(p));
220 0 : return umodiu(a, p);
221 : }
222 0 : default: pari_err_TYPE("Rg_to_Fl",x);
223 : return 0; /* LCOV_EXCL_LINE */
224 : }
225 : }
226 :
227 : ulong
228 1100606 : Rg_to_F2(GEN x)
229 : {
230 1100606 : switch(typ(x))
231 : {
232 258163 : case t_INT: return mpodd(x);
233 0 : case t_FRAC:
234 0 : if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
235 0 : return mpodd(gel(x,1));
236 0 : case t_PADIC:
237 0 : if (!absequaliu(gel(x,2),2)) pari_err_OP("",x, mkintmodu(1,2));
238 0 : if (valp(x) < 0) (void)Fl_inv(0,2);
239 0 : return valp(x) & 1;
240 842443 : case t_INTMOD: {
241 842443 : GEN q = gel(x,1), a = gel(x,2);
242 842443 : if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
243 842443 : return mpodd(a);
244 : }
245 0 : default: pari_err_TYPE("Rg_to_F2",x);
246 : return 0; /* LCOV_EXCL_LINE */
247 : }
248 : }
249 :
250 : GEN
251 1989149 : RgX_to_Flx(GEN x, ulong p)
252 : {
253 1989149 : long i, lx = lg(x);
254 1989149 : GEN a = cgetg(lx, t_VECSMALL);
255 1989149 : a[1]=((ulong)x[1])&VARNBITS;
256 17179384 : for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
257 1989149 : return Flx_renormalize(a,lx);
258 : }
259 :
260 : GEN
261 0 : RgXV_to_FlxV(GEN x, ulong p)
262 0 : { pari_APPLY_type(t_VEC, RgX_to_Flx(gel(x,i), p)) }
263 :
264 : /* If x is a POLMOD, assume modulus is a multiple of T. */
265 : GEN
266 1981882 : Rg_to_Flxq(GEN x, GEN T, ulong p)
267 : {
268 1981882 : long ta, tx = typ(x), v = get_Flx_var(T);
269 : GEN a, b;
270 1981886 : if (is_const_t(tx))
271 : {
272 1909564 : if (tx == t_FFELT) return FF_to_Flxq(x);
273 1183340 : return Fl_to_Flx(Rg_to_Fl(x, p), v);
274 : }
275 72325 : switch(tx)
276 : {
277 8576 : case t_POLMOD:
278 8576 : b = gel(x,1);
279 8576 : a = gel(x,2); ta = typ(a);
280 8576 : if (is_const_t(ta)) return Fl_to_Flx(Rg_to_Fl(a, p), v);
281 8422 : b = RgX_to_Flx(b, p); if (b[1] != v) break;
282 8422 : a = RgX_to_Flx(a, p); if (Flx_equal(b,T)) return a;
283 0 : if (lgpol(Flx_rem(b,T,p))==0) return Flx_rem(a, T, p);
284 0 : break;
285 63749 : case t_POL:
286 63749 : x = RgX_to_Flx(x,p);
287 63749 : if (x[1] != v) break;
288 63749 : return Flx_rem(x, T, p);
289 0 : case t_RFRAC:
290 0 : a = Rg_to_Flxq(gel(x,1), T,p);
291 0 : b = Rg_to_Flxq(gel(x,2), T,p);
292 0 : return Flxq_div(a,b, T,p);
293 : }
294 0 : pari_err_TYPE("Rg_to_Flxq",x);
295 : return NULL; /* LCOV_EXCL_LINE */
296 : }
297 :
298 : /***********************************************************************/
299 : /** **/
300 : /** Basic operation on Flx **/
301 : /** **/
302 : /***********************************************************************/
303 : /* = zx_renormalize. Similar to normalizepol, in place */
304 : GEN
305 1235022348 : Flx_renormalize(GEN /*in place*/ x, long lx)
306 : {
307 : long i;
308 1382727264 : for (i = lx-1; i>1; i--)
309 1321423955 : if (x[i]) break;
310 1235022348 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
311 1234734455 : setlg(x, i+1); return x;
312 : }
313 :
314 : GEN
315 306445 : Flx_red(GEN z, ulong p)
316 : {
317 306445 : long i, l = lg(z);
318 306445 : GEN x = cgetg(l, t_VECSMALL);
319 306363 : x[1] = z[1];
320 17879510 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
321 306363 : return Flx_renormalize(x,l);
322 : }
323 :
324 : int
325 17568233 : Flx_equal(GEN V, GEN W)
326 : {
327 17568233 : long l = lg(V);
328 17568233 : if (lg(W) != l) return 0;
329 17777149 : while (--l > 1) /* do not compare variables, V[1] */
330 17256590 : if (V[l] != W[l]) return 0;
331 520559 : return 1;
332 : }
333 :
334 : GEN
335 1314172 : random_Flx(long d1, long vs, ulong p)
336 : {
337 1314172 : long i, d = d1+2;
338 1314172 : GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
339 10648788 : for (i=2; i<d; i++) y[i] = random_Fl(p);
340 1314514 : return Flx_renormalize(y,d);
341 : }
342 :
343 : static GEN
344 4466804 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
345 : {
346 : long i,lz;
347 : GEN z;
348 :
349 4466804 : if (ly>lx) swapspec(x,y, lx,ly);
350 4466804 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
351 80363902 : for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
352 68169907 : for ( ; i<lx; i++) z[i+2] = x[i];
353 4466804 : z[1] = 0; return Flx_renormalize(z, lz);
354 : }
355 :
356 : GEN
357 36435577 : Flx_add(GEN x, GEN y, ulong p)
358 : {
359 : long i,lz;
360 : GEN z;
361 36435577 : long lx=lg(x);
362 36435577 : long ly=lg(y);
363 36435577 : if (ly>lx) swapspec(x,y, lx,ly);
364 36435577 : lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
365 344969495 : for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
366 71480364 : for ( ; i<lx; i++) z[i] = x[i];
367 36394020 : return Flx_renormalize(z, lz);
368 : }
369 :
370 : GEN
371 7580675 : Flx_Fl_add(GEN y, ulong x, ulong p)
372 : {
373 : GEN z;
374 : long lz, i;
375 7580675 : if (!lgpol(y))
376 225998 : return Fl_to_Flx(x,y[1]);
377 7355919 : lz=lg(y);
378 7355919 : z=cgetg(lz,t_VECSMALL);
379 7354343 : z[1]=y[1];
380 7354343 : z[2] = Fl_add(y[2],x,p);
381 40126342 : for(i=3;i<lz;i++)
382 32772144 : z[i] = y[i];
383 7354198 : if (lz==3) z = Flx_renormalize(z,lz);
384 7354130 : return z;
385 : }
386 :
387 : static GEN
388 701865 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
389 : {
390 : long i,lz;
391 : GEN z;
392 :
393 701865 : if (ly <= lx)
394 : {
395 701933 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
396 40196640 : for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
397 2226170 : for ( ; i<lx; i++) z[i+2] = x[i];
398 : }
399 : else
400 : {
401 0 : lz = ly+2; z = cgetg(lz, t_VECSMALL);
402 0 : for (i=0; i<lx; i++) z[i+2] = Fl_sub(x[i],y[i],p);
403 0 : for ( ; i<ly; i++) z[i+2] = Fl_neg(y[i],p);
404 : }
405 701584 : z[1] = 0; return Flx_renormalize(z, lz);
406 : }
407 :
408 : GEN
409 74993728 : Flx_sub(GEN x, GEN y, ulong p)
410 : {
411 74993728 : long i,lz,lx = lg(x), ly = lg(y);
412 : GEN z;
413 :
414 74993728 : if (ly <= lx)
415 : {
416 43563192 : lz = lx; z = cgetg(lz, t_VECSMALL);
417 259472370 : for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
418 88931507 : for ( ; i<lx; i++) z[i] = x[i];
419 : }
420 : else
421 : {
422 31430536 : lz = ly; z = cgetg(lz, t_VECSMALL);
423 242355753 : for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
424 143744225 : for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
425 : }
426 74984368 : z[1]=x[1]; return Flx_renormalize(z, lz);
427 : }
428 :
429 : GEN
430 28332 : Flx_Fl_sub(GEN y, ulong x, ulong p)
431 : {
432 : GEN z;
433 28332 : long lz = lg(y), i;
434 28332 : if (lz==2)
435 72 : return Fl_to_Flx(Fl_neg(x, p),y[1]);
436 28260 : z = cgetg(lz, t_VECSMALL);
437 28260 : z[1] = y[1];
438 28260 : z[2] = Fl_sub(uel(y,2), x, p);
439 225004 : for(i=3; i<lz; i++)
440 196744 : z[i] = y[i];
441 28260 : if (lz==3) z = Flx_renormalize(z,lz);
442 28260 : return z;
443 : }
444 :
445 : static GEN
446 1638195 : Flx_negspec(GEN x, ulong p, long l)
447 : {
448 : long i;
449 1638195 : GEN z = cgetg(l+2, t_VECSMALL) + 2;
450 10298730 : for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
451 1638311 : return z-2;
452 : }
453 :
454 : GEN
455 1638219 : Flx_neg(GEN x, ulong p)
456 : {
457 1638219 : GEN z = Flx_negspec(x+2, p, lgpol(x));
458 1638583 : z[1] = x[1];
459 1638583 : return z;
460 : }
461 :
462 : GEN
463 1091723 : Flx_neg_inplace(GEN x, ulong p)
464 : {
465 1091723 : long i, l = lg(x);
466 39757626 : for (i=2; i<l; i++)
467 38665903 : if (x[i]) x[i] = p - x[i];
468 1091723 : return x;
469 : }
470 :
471 : GEN
472 1823823 : Flx_double(GEN y, ulong p)
473 : {
474 : long i, l;
475 1823823 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
476 15674622 : for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
477 1823823 : return Flx_renormalize(z, l);
478 : }
479 : GEN
480 648071 : Flx_triple(GEN y, ulong p)
481 : {
482 : long i, l;
483 648071 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
484 5715917 : for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
485 648071 : return Flx_renormalize(z, l);
486 : }
487 : GEN
488 10000985 : Flx_Fl_mul(GEN y, ulong x, ulong p)
489 : {
490 : GEN z;
491 : long i, l;
492 10000985 : if (!x) return pol0_Flx(y[1]);
493 9297728 : z = cgetg_copy(y, &l); z[1] = y[1];
494 9296965 : if (HIGHWORD(x | p))
495 5574307 : for(i=2; i<l; i++) z[i] = Fl_mul(y[i], x, p);
496 : else
497 47270750 : for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
498 9297132 : return Flx_renormalize(z, l);
499 : }
500 : GEN
501 6954439 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
502 : {
503 : GEN z;
504 : long i, l;
505 6954439 : z = cgetg_copy(y, &l); z[1] = y[1];
506 6951074 : if (HIGHWORD(x | p))
507 4623967 : for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
508 : else
509 13379413 : for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
510 6951063 : z[l-1] = 1; return z;
511 : }
512 :
513 : /* Return a*x^n if n>=0 and a\x^(-n) if n<0 */
514 : GEN
515 7792227 : Flx_shift(GEN a, long n)
516 : {
517 7792227 : long i, l = lg(a);
518 : GEN b;
519 7792227 : if (l==2 || !n) return Flx_copy(a);
520 7756766 : if (l+n<=2) return pol0_Flx(a[1]);
521 7745917 : b = cgetg(l+n, t_VECSMALL);
522 7743077 : b[1] = a[1];
523 7743077 : if (n < 0)
524 29697444 : for (i=2-n; i<l; i++) b[i+n] = a[i];
525 : else
526 : {
527 13901066 : for (i=0; i<n; i++) b[2+i] = 0;
528 55551945 : for (i=2; i<l; i++) b[i+n] = a[i];
529 : }
530 7743077 : return b;
531 : }
532 :
533 : GEN
534 39429482 : Flx_normalize(GEN z, ulong p)
535 : {
536 39429482 : long l = lg(z)-1;
537 39429482 : ulong p1 = z[l]; /* leading term */
538 39429482 : if (p1 == 1) return z;
539 6940998 : return Flx_Fl_mul_to_monic(z, Fl_inv(p1,p), p);
540 : }
541 :
542 : /* return (x * X^d) + y. Assume d > 0, shallow if x == 0*/
543 : static GEN
544 2313385 : Flx_addshift(GEN x, GEN y, ulong p, long d)
545 : {
546 2313385 : GEN xd,yd,zd = (GEN)avma;
547 2313385 : long a,lz,ny = lgpol(y), nx = lgpol(x);
548 2313385 : long vs = x[1];
549 2313385 : if (nx == 0) return y;
550 2311533 : x += 2; y += 2; a = ny-d;
551 2311533 : if (a <= 0)
552 : {
553 81416 : lz = (a>nx)? ny+2: nx+d+2;
554 81416 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
555 1611331 : while (xd > x) *--zd = *--xd;
556 81416 : x = zd + a;
557 148444 : while (zd > x) *--zd = 0;
558 : }
559 : else
560 : {
561 2230117 : xd = new_chunk(d); yd = y+d;
562 2230117 : x = Flx_addspec(x,yd,p, nx,a);
563 2230117 : lz = (a>nx)? ny+2: lg(x)+d;
564 101159195 : x += 2; while (xd > x) *--zd = *--xd;
565 : }
566 45535651 : while (yd > y) *--zd = *--yd;
567 2311533 : *--zd = vs;
568 2311533 : *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
569 : }
570 :
571 : /* shift polynomial + gerepile */
572 : /* Do not set evalvarn*/
573 : static GEN
574 396325739 : Flx_shiftip(pari_sp av, GEN x, long v)
575 : {
576 396325739 : long i, lx = lg(x), ly;
577 : GEN y;
578 396325739 : if (!v || lx==2) return gerepileuptoleaf(av, x);
579 101378147 : ly = lx + v; /* result length */
580 101378147 : (void)new_chunk(ly); /* check that result fits */
581 101312381 : x += lx; y = (GEN)av;
582 822112839 : for (i = 2; i<lx; i++) *--y = *--x;
583 437265742 : for (i = 0; i< v; i++) *--y = 0;
584 101312381 : y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
585 101485473 : set_avma((pari_sp)y); return y;
586 : }
587 :
588 : static long
589 1325251914 : get_Fl_threshold(ulong p, long mul, long mul2)
590 : {
591 1325251914 : return SMALL_ULONG(p) ? mul: mul2;
592 : }
593 :
594 : #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
595 : #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
596 : #define LLQUARTWORD(x) ((x) & QUARTMASK)
597 : #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
598 : #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
599 : #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
600 : INLINE long
601 5174235 : maxbitcoeffpol(ulong p, long n)
602 : {
603 5174235 : GEN z = muliu(sqru(p - 1), n);
604 5167506 : long b = expi(z) + 1;
605 : /* only do expensive bit-packing if it saves at least 1 limb */
606 5172854 : if (b <= BITS_IN_QUARTULONG)
607 : {
608 614782 : if (nbits2nlong(n*b) == (n + 3)>>2)
609 45872 : b = BITS_IN_QUARTULONG;
610 : }
611 4558072 : else if (b <= BITS_IN_HALFULONG)
612 : {
613 976902 : if (nbits2nlong(n*b) == (n + 1)>>1)
614 5434 : b = BITS_IN_HALFULONG;
615 : }
616 : else
617 : {
618 3581170 : long l = lgefint(z) - 2;
619 3581170 : if (nbits2nlong(n*b) == n*l)
620 286024 : b = l*BITS_IN_LONG;
621 : }
622 5172922 : return b;
623 : }
624 :
625 : INLINE ulong
626 2249500438 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
627 : { /* Assume OK_ULONG*/
628 2249500438 : ulong p1 = 0;
629 : long i;
630 10419088634 : for (i=a; i<b; i++)
631 8169588196 : if (y[i])
632 : {
633 6858714667 : p1 += y[i] * x[-i];
634 6858714667 : if (p1 & HIGHBIT) p1 %= p;
635 : }
636 2249500438 : return p1 % p;
637 : }
638 :
639 : INLINE ulong
640 921037265 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
641 : {
642 921037265 : ulong p1 = 0;
643 : long i;
644 2880767097 : for (i=a; i<b; i++)
645 1958722255 : if (y[i])
646 1921880167 : p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
647 922044842 : return p1;
648 : }
649 :
650 : /* assume nx >= ny > 0 */
651 : static GEN
652 205814488 : Flx_mulspec_basecase(GEN x, GEN y, ulong p, long nx, long ny)
653 : {
654 : long i,lz,nz;
655 : GEN z;
656 :
657 205814488 : lz = nx+ny+1; nz = lz-2;
658 205814488 : z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
659 205493128 : if (SMALL_ULONG(p))
660 : {
661 682781976 : for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
662 527795279 : for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
663 534062311 : for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
664 : }
665 : else
666 : {
667 55582755 : ulong pi = get_Fl_red(p);
668 200765951 : for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
669 151476803 : for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
670 145620854 : for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
671 : }
672 205581263 : z -= 2; return Flx_renormalize(z, lz);
673 : }
674 :
675 : static GEN
676 12108 : int_to_Flx(GEN z, ulong p)
677 : {
678 12108 : long i, l = lgefint(z);
679 12108 : GEN x = cgetg(l, t_VECSMALL);
680 1008693 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
681 12093 : return Flx_renormalize(x, l);
682 : }
683 :
684 : INLINE GEN
685 9840 : Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
686 : {
687 9840 : GEN z=muliispec(a,b,na,nb);
688 9850 : return int_to_Flx(z,p);
689 : }
690 :
691 : static GEN
692 398002 : Flx_to_int_halfspec(GEN a, long na)
693 : {
694 : long j;
695 398002 : long n = (na+1)>>1UL;
696 398002 : GEN V = cgetipos(2+n);
697 : GEN w;
698 1254994 : for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
699 856992 : *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
700 398002 : if (j<na)
701 268576 : *w = a[j];
702 398002 : return V;
703 : }
704 :
705 : static GEN
706 437010 : int_to_Flx_half(GEN z, ulong p)
707 : {
708 : long i;
709 437010 : long lx = (lgefint(z)-2)*2+2;
710 437010 : GEN w, x = cgetg(lx, t_VECSMALL);
711 1640505 : for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
712 : {
713 1203495 : x[i] = LOWWORD((ulong)*w)%p;
714 1203495 : x[i+1] = HIGHWORD((ulong)*w)%p;
715 : }
716 437010 : return Flx_renormalize(x, lx);
717 : }
718 :
719 : static GEN
720 5298 : Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
721 : {
722 5298 : GEN A = Flx_to_int_halfspec(a,na);
723 5298 : GEN B = Flx_to_int_halfspec(b,nb);
724 5298 : GEN z = mulii(A,B);
725 5298 : return int_to_Flx_half(z,p);
726 : }
727 :
728 : static GEN
729 85030 : Flx_to_int_quartspec(GEN a, long na)
730 : {
731 : long j;
732 85030 : long n = (na+3)>>2UL;
733 85030 : GEN V = cgetipos(2+n);
734 : GEN w;
735 2798209 : for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
736 2713179 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
737 85030 : switch (na-j)
738 : {
739 30168 : case 3:
740 30168 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
741 30168 : break;
742 21196 : case 2:
743 21196 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
744 21196 : break;
745 14066 : case 1:
746 14066 : *w = a[j];
747 14066 : break;
748 19600 : case 0:
749 19600 : break;
750 : }
751 85030 : return V;
752 : }
753 :
754 : static GEN
755 45872 : int_to_Flx_quart(GEN z, ulong p)
756 : {
757 : long i;
758 45872 : long lx = (lgefint(z)-2)*4+2;
759 45872 : GEN w, x = cgetg(lx, t_VECSMALL);
760 3181020 : for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
761 : {
762 3135148 : x[i] = LLQUARTWORD((ulong)*w)%p;
763 3135148 : x[i+1] = HLQUARTWORD((ulong)*w)%p;
764 3135148 : x[i+2] = LHQUARTWORD((ulong)*w)%p;
765 3135148 : x[i+3] = HHQUARTWORD((ulong)*w)%p;
766 : }
767 45872 : return Flx_renormalize(x, lx);
768 : }
769 :
770 : static GEN
771 39158 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
772 : {
773 39158 : GEN A = Flx_to_int_quartspec(a,na);
774 39158 : GEN B = Flx_to_int_quartspec(b,nb);
775 39158 : GEN z = mulii(A,B);
776 39158 : return int_to_Flx_quart(z,p);
777 : }
778 :
779 : /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
780 : static GEN
781 540506 : Flx_eval2BILspec(GEN x, long k, long l)
782 : {
783 540506 : long i, lz = k*l, ki;
784 540506 : GEN pz = cgetipos(2+lz);
785 14966918 : for (i=0; i < lz; i++)
786 14426412 : *int_W(pz,i) = 0UL;
787 7753712 : for (i=0, ki=0; i<l; i++, ki+=k)
788 7213206 : *int_W(pz,ki) = x[i];
789 540506 : return int_normalize(pz,0);
790 : }
791 :
792 : static GEN
793 276849 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
794 : {
795 276849 : long i, offset, lm = lgefint(x)-2, l = d+3;
796 276849 : ulong pi = get_Fl_red(p);
797 276849 : GEN pol = cgetg(l, t_VECSMALL);
798 276849 : pol[1] = 0;
799 7316832 : for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
800 7039983 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
801 276849 : if (offset < lm)
802 211562 : pol[i+2] = (*int_W(x,offset)) % p;
803 276849 : return Flx_renormalize(pol,l);
804 : }
805 :
806 : static GEN
807 0 : Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
808 : {
809 0 : long i, offset, lm = lgefint(x)-2, l = d+3;
810 0 : ulong pi = get_Fl_red(p);
811 0 : GEN pol = cgetg(l, t_VECSMALL);
812 0 : pol[1] = 0;
813 0 : for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
814 0 : pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
815 0 : *int_W(x,offset), p, pi);
816 0 : if (offset+1 < lm)
817 0 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
818 0 : else if (offset < lm)
819 0 : pol[i+2] = (*int_W(x,offset)) % p;
820 0 : return Flx_renormalize(pol,l);
821 : }
822 :
823 : static GEN
824 273919 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
825 : {
826 273919 : return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
827 : }
828 :
829 : static GEN
830 263198 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
831 : {
832 263198 : pari_sp av = avma;
833 263198 : GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
834 263198 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
835 : }
836 :
837 : static GEN
838 14095617 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
839 : GEN y;
840 : long i;
841 14095617 : if (l == 0)
842 3305722 : return gen_0;
843 10789895 : y = cgetg(l + 1, t_VECSMALL);
844 518131338 : for(i = 1; i <= l; i++)
845 507352353 : y[i] = x[l - i];
846 10778985 : return nv_fromdigits_2k(y, b);
847 : }
848 :
849 : /* assume b < BITS_IN_LONG */
850 : static GEN
851 4332974 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
852 4332974 : GEN v = binary_2k_nv(z, b), x;
853 4332909 : long i, l = lg(v) + 1;
854 4332909 : x = cgetg(l, t_VECSMALL);
855 412142516 : for (i = 2; i < l; i++)
856 407809104 : x[i] = v[l - i] % p;
857 4333412 : return Flx_renormalize(x, l);
858 : }
859 :
860 : static GEN
861 3234155 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
862 3234155 : GEN v = binary_2k(z, b), x, y;
863 3233867 : long i, l = lg(v) + 1, ly;
864 3233867 : x = cgetg(l, t_VECSMALL);
865 136795588 : for (i = 2; i < l; i++) {
866 133559650 : y = gel(v, l - i);
867 133559650 : ly = lgefint(y);
868 133559650 : switch (ly) {
869 4948658 : case 2: x[i] = 0; break;
870 14889016 : case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
871 105417924 : case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
872 24912132 : case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
873 8304052 : *int_W_lg(y, 0, ly), p, pi); break;
874 0 : default: x[i] = umodiu(gel(v, l - i), p);
875 : }
876 : }
877 3235938 : return Flx_renormalize(x, l);
878 : }
879 :
880 : static GEN
881 4259960 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
882 : {
883 : GEN C, D;
884 4259960 : pari_sp av = avma;
885 4259960 : A = kron_pack_Flx_spec_bits(A, b, lA);
886 4263074 : B = kron_pack_Flx_spec_bits(B, b, lB);
887 4263196 : C = gerepileuptoint(av, mulii(A, B));
888 4261903 : if (b < BITS_IN_LONG)
889 1349091 : D = kron_unpack_Flx_bits_narrow(C, b, p);
890 : else
891 : {
892 2912812 : ulong pi = get_Fl_red(p);
893 2911197 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
894 : }
895 4261197 : return D;
896 : }
897 :
898 : static GEN
899 574970 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
900 : {
901 : GEN C, D;
902 574970 : A = kron_pack_Flx_spec_bits(A, b, lA);
903 574994 : C = sqri(A);
904 575003 : if (b < BITS_IN_LONG)
905 381213 : D = kron_unpack_Flx_bits_narrow(C, b, p);
906 : else
907 : {
908 193790 : ulong pi = get_Fl_red(p);
909 193787 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
910 : }
911 574982 : return D;
912 : }
913 :
914 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
915 : * b+2 were sent instead. na, nb = number of terms of a, b.
916 : * Only c, c0, c1, c2 are genuine GEN.
917 : */
918 : static GEN
919 228222775 : Flx_mulspec(GEN a, GEN b, ulong p, long na, long nb)
920 : {
921 : GEN a0,c,c0;
922 228222775 : long n0, n0a, i, v = 0;
923 : pari_sp av;
924 :
925 299145224 : while (na && !a[0]) { a++; na--; v++; }
926 362377158 : while (nb && !b[0]) { b++; nb--; v++; }
927 228222775 : if (na < nb) swapspec(a,b, na,nb);
928 228222775 : if (!nb) return pol0_Flx(0);
929 :
930 211287321 : av = avma;
931 211287321 : if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
932 : {
933 4579568 : long m = maxbitcoeffpol(p,nb);
934 4577418 : switch (m)
935 : {
936 39158 : case BITS_IN_QUARTULONG:
937 39158 : return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
938 5298 : case BITS_IN_HALFULONG:
939 5298 : return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
940 9840 : case BITS_IN_LONG:
941 9840 : return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
942 263198 : case 2*BITS_IN_LONG:
943 263198 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
944 0 : case 3*BITS_IN_LONG:
945 0 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
946 4259924 : default:
947 4259924 : return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
948 : }
949 : }
950 206861917 : if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
951 205796285 : return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,na,nb), v);
952 1115161 : i=(na>>1); n0=na-i; na=i;
953 1115161 : a0=a+n0; n0a=n0;
954 1784907 : while (n0a && !a[n0a-1]) n0a--;
955 :
956 1115161 : if (nb > n0)
957 : {
958 : GEN b0,c1,c2;
959 : long n0b;
960 :
961 1091723 : nb -= n0; b0 = b+n0; n0b = n0;
962 2096084 : while (n0b && !b[n0b-1]) n0b--;
963 1091723 : c = Flx_mulspec(a,b,p,n0a,n0b);
964 1091723 : c0 = Flx_mulspec(a0,b0,p,na,nb);
965 :
966 1091723 : c2 = Flx_addspec(a0,a,p,na,n0a);
967 1091723 : c1 = Flx_addspec(b0,b,p,nb,n0b);
968 :
969 1091723 : c1 = Flx_mul(c1,c2,p);
970 1091723 : c2 = Flx_add(c0,c,p);
971 :
972 1091723 : c2 = Flx_neg_inplace(c2,p);
973 1091723 : c2 = Flx_add(c1,c2,p);
974 1091723 : c0 = Flx_addshift(c0,c2 ,p, n0);
975 : }
976 : else
977 : {
978 23438 : c = Flx_mulspec(a,b,p,n0a,nb);
979 23438 : c0 = Flx_mulspec(a0,b,p,na,nb);
980 : }
981 1115161 : c0 = Flx_addshift(c0,c,p,n0);
982 1115161 : return Flx_shiftip(av,c0, v);
983 : }
984 :
985 : GEN
986 224319804 : Flx_mul(GEN x, GEN y, ulong p)
987 : {
988 224319804 : GEN z = Flx_mulspec(x+2,y+2,p, lgpol(x),lgpol(y));
989 224504760 : z[1] = x[1]; return z;
990 : }
991 :
992 : static GEN
993 186096039 : Flx_sqrspec_basecase(GEN x, ulong p, long nx)
994 : {
995 : long i, lz, nz;
996 : ulong p1;
997 : GEN z;
998 :
999 186096039 : if (!nx) return pol0_Flx(0);
1000 186096039 : lz = (nx << 1) + 1, nz = lz-2;
1001 186096039 : z = cgetg(lz, t_VECSMALL) + 2;
1002 185539792 : if (SMALL_ULONG(p))
1003 : {
1004 128820569 : z[0] = x[0]*x[0]%p;
1005 617708762 : for (i=1; i<nx; i++)
1006 : {
1007 489125845 : p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1008 488888193 : p1 <<= 1;
1009 488888193 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1010 488888193 : z[i] = p1 % p;
1011 : }
1012 620947736 : for ( ; i<nz; i++)
1013 : {
1014 491958032 : p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1015 492364819 : p1 <<= 1;
1016 492364819 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1017 492364819 : z[i] = p1 % p;
1018 : }
1019 : }
1020 : else
1021 : {
1022 56719223 : ulong pi = get_Fl_red(p);
1023 56726386 : z[0] = Fl_sqr_pre(x[0], p, pi);
1024 358361825 : for (i=1; i<nx; i++)
1025 : {
1026 301643879 : p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1027 302023529 : p1 = Fl_add(p1, p1, p);
1028 301690862 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1029 301458331 : z[i] = p1;
1030 : }
1031 358405731 : for ( ; i<nz; i++)
1032 : {
1033 301578847 : p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1034 302500724 : p1 = Fl_add(p1, p1, p);
1035 302227853 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1036 301687785 : z[i] = p1;
1037 : }
1038 : }
1039 185816588 : z -= 2; return Flx_renormalize(z, lz);
1040 : }
1041 :
1042 : static GEN
1043 2264 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
1044 : {
1045 2264 : GEN z=sqrispec(a,na);
1046 2265 : return int_to_Flx(z,p);
1047 : }
1048 :
1049 : static GEN
1050 136 : Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
1051 : {
1052 136 : GEN z = sqri(Flx_to_int_halfspec(a,na));
1053 136 : return int_to_Flx_half(z,p);
1054 : }
1055 :
1056 : static GEN
1057 6714 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1058 : {
1059 6714 : GEN z = sqri(Flx_to_int_quartspec(a,na));
1060 6714 : return int_to_Flx_quart(z,p);
1061 : }
1062 :
1063 : static GEN
1064 10721 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1065 : {
1066 10721 : pari_sp av = avma;
1067 10721 : GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1068 10721 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1069 : }
1070 :
1071 : static GEN
1072 185995610 : Flx_sqrspec(GEN a, ulong p, long na)
1073 : {
1074 : GEN a0, c, c0;
1075 185995610 : long n0, n0a, i, v = 0, m;
1076 : pari_sp av;
1077 :
1078 252187963 : while (na && !a[0]) { a++; na--; v += 2; }
1079 185995610 : if (!na) return pol0_Flx(0);
1080 :
1081 185895891 : av = avma;
1082 185895891 : if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1083 : {
1084 594807 : m = maxbitcoeffpol(p,na);
1085 594800 : switch(m)
1086 : {
1087 6714 : case BITS_IN_QUARTULONG:
1088 6714 : return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
1089 136 : case BITS_IN_HALFULONG:
1090 136 : return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
1091 2264 : case BITS_IN_LONG:
1092 2264 : return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1093 10721 : case 2*BITS_IN_LONG:
1094 10721 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
1095 0 : case 3*BITS_IN_LONG:
1096 0 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
1097 574965 : default:
1098 574965 : return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1099 : }
1100 : }
1101 185738421 : if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1102 185943762 : return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,na), v);
1103 53260 : i=(na>>1); n0=na-i; na=i;
1104 53260 : a0=a+n0; n0a=n0;
1105 66845 : while (n0a && !a[n0a-1]) n0a--;
1106 :
1107 53260 : c = Flx_sqrspec(a,p,n0a);
1108 53260 : c0= Flx_sqrspec(a0,p,na);
1109 53260 : if (p == 2) n0 *= 2;
1110 : else
1111 : {
1112 53241 : GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1113 53241 : t = Flx_sqr(t,p);
1114 53241 : c1= Flx_add(c0,c, p);
1115 53241 : c1= Flx_sub(t, c1, p);
1116 53241 : c0 = Flx_addshift(c0,c1,p,n0);
1117 : }
1118 53260 : c0 = Flx_addshift(c0,c,p,n0);
1119 53260 : return Flx_shiftip(av,c0,v);
1120 : }
1121 :
1122 : GEN
1123 185715874 : Flx_sqr(GEN x, ulong p)
1124 : {
1125 185715874 : GEN z = Flx_sqrspec(x+2,p, lgpol(x));
1126 186710810 : z[1] = x[1]; return z;
1127 : }
1128 :
1129 : GEN
1130 5179 : Flx_powu(GEN x, ulong n, ulong p)
1131 : {
1132 5179 : GEN y = pol1_Flx(x[1]), z;
1133 : ulong m;
1134 5169 : if (n == 0) return y;
1135 5169 : m = n; z = x;
1136 : for (;;)
1137 : {
1138 19835 : if (m&1UL) y = Flx_mul(y,z, p);
1139 19834 : m >>= 1; if (!m) return y;
1140 14671 : z = Flx_sqr(z, p);
1141 : }
1142 : }
1143 :
1144 : GEN
1145 13145 : Flx_halve(GEN y, ulong p)
1146 : {
1147 : GEN z;
1148 : long i, l;
1149 13145 : z = cgetg_copy(y, &l); z[1] = y[1];
1150 54932 : for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
1151 13145 : return z;
1152 : }
1153 :
1154 : static GEN
1155 3086214 : Flx_recipspec(GEN x, long l, long n)
1156 : {
1157 : long i;
1158 3086214 : GEN z=cgetg(n+2,t_VECSMALL)+2;
1159 80252080 : for(i=0; i<l; i++)
1160 77167276 : z[n-i-1] = x[i];
1161 7122064 : for( ; i<n; i++)
1162 4037260 : z[n-i-1] = 0;
1163 3084804 : return Flx_renormalize(z-2,n+2);
1164 : }
1165 :
1166 : GEN
1167 0 : Flx_recip(GEN x)
1168 : {
1169 0 : GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
1170 0 : z[1]=x[1];
1171 0 : return z;
1172 : }
1173 :
1174 : /* Return h^degpol(P) P(x / h) */
1175 : GEN
1176 1117 : Flx_rescale(GEN P, ulong h, ulong p)
1177 : {
1178 1117 : long i, l = lg(P);
1179 1117 : GEN Q = cgetg(l,t_VECSMALL);
1180 1117 : ulong hi = h;
1181 1117 : Q[l-1] = P[l-1];
1182 12524 : for (i=l-2; i>=2; i--)
1183 : {
1184 12523 : Q[i] = Fl_mul(P[i], hi, p);
1185 12523 : if (i == 2) break;
1186 11406 : hi = Fl_mul(hi,h, p);
1187 : }
1188 1118 : Q[1] = P[1]; return Q;
1189 : }
1190 :
1191 : /*
1192 : * x/polrecip(P)+O(x^n)
1193 : */
1194 : static GEN
1195 118265 : Flx_invBarrett_basecase(GEN T, ulong p)
1196 : {
1197 118265 : long i, l=lg(T)-1, lr=l-1, k;
1198 118265 : GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1199 118265 : r[2] = 1;
1200 118265 : if (SMALL_ULONG(p))
1201 452948 : for (i=3;i<lr;i++)
1202 : {
1203 449100 : ulong u = uel(T, l-i+2);
1204 28843540 : for (k=3; k<i; k++)
1205 28394440 : { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1206 449100 : r[i] = Fl_neg(u % p, p);
1207 : }
1208 : else
1209 1896099 : for (i=3;i<lr;i++)
1210 : {
1211 1781683 : ulong u = Fl_neg(uel(T,l-i+2), p);
1212 46707869 : for (k=3; k<i; k++)
1213 44926187 : u = Fl_sub(u, Fl_mul(uel(T,l-i+k), uel(r,k), p), p);
1214 1781682 : r[i] = u;
1215 : }
1216 118264 : return Flx_renormalize(r,lr);
1217 : }
1218 :
1219 : /* Return new lgpol */
1220 : static long
1221 1698827 : Flx_lgrenormalizespec(GEN x, long lx)
1222 : {
1223 : long i;
1224 11073545 : for (i = lx-1; i>=0; i--)
1225 11073670 : if (x[i]) break;
1226 1698827 : return i+1;
1227 : }
1228 : static GEN
1229 20552 : Flx_invBarrett_Newton(GEN T, ulong p)
1230 : {
1231 20552 : long nold, lx, lz, lq, l = degpol(T), lQ;
1232 20552 : GEN q, y, z, x = zero_zv(l+1) + 2;
1233 20554 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1234 : pari_sp av;
1235 :
1236 20554 : y = T+2;
1237 20554 : q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1238 20553 : av = avma;
1239 : /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1240 :
1241 : /* initialize */
1242 20553 : x[0] = Fl_inv(q[0], p);
1243 20554 : if (lQ>1 && q[1])
1244 5016 : {
1245 5016 : ulong u = q[1];
1246 5016 : if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1247 5016 : x[1] = p - u; lx = 2;
1248 : }
1249 : else
1250 15538 : lx = 1;
1251 20554 : nold = 1;
1252 144327 : for (; mask > 1; set_avma(av))
1253 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1254 123780 : long i, lnew, nnew = nold << 1;
1255 :
1256 123780 : if (mask & 1) nnew--;
1257 123780 : mask >>= 1;
1258 :
1259 123780 : lnew = nnew + 1;
1260 123780 : lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1261 123832 : z = Flx_mulspec(x, q, p, lx, lq); /* FIXME: high product */
1262 123771 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1263 123775 : z += 2;
1264 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1265 276506 : for (i = nold; i < lz; i++) if (z[i]) break;
1266 123775 : nold = nnew;
1267 123775 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1268 :
1269 : /* z + i represents (x*q - 1) / t^i */
1270 85967 : lz = Flx_lgrenormalizespec (z+i, lz-i);
1271 85966 : z = Flx_mulspec(x, z+i, p, lx, lz); /* FIXME: low product */
1272 85967 : lz = lgpol(z); z += 2;
1273 85967 : if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1274 :
1275 85966 : lx = lz+ i;
1276 85966 : y = x + i; /* x -= z * t^i, in place */
1277 807059 : for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1278 : }
1279 20554 : x -= 2; setlg(x, lx + 2); x[1] = T[1];
1280 20554 : return x;
1281 : }
1282 :
1283 : GEN
1284 139518 : Flx_invBarrett(GEN T, ulong p)
1285 : {
1286 139518 : pari_sp ltop = avma;
1287 139518 : long l = lgpol(T);
1288 : GEN r;
1289 139518 : if (l < 3) return pol0_Flx(T[1]);
1290 138817 : if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1291 : {
1292 118265 : ulong c = T[l+1];
1293 118265 : if (c!=1)
1294 : {
1295 97452 : ulong ci = Fl_inv(c,p);
1296 97452 : T=Flx_Fl_mul(T, ci, p);
1297 97452 : r=Flx_invBarrett_basecase(T,p);
1298 97452 : r=Flx_Fl_mul(r,ci,p);
1299 : }
1300 : else
1301 20813 : r=Flx_invBarrett_basecase(T,p);
1302 : }
1303 : else
1304 20552 : r = Flx_invBarrett_Newton(T,p);
1305 138819 : return gerepileuptoleaf(ltop, r);
1306 : }
1307 :
1308 : GEN
1309 45007868 : Flx_get_red(GEN T, ulong p)
1310 : {
1311 45007868 : if (typ(T)!=t_VECSMALL
1312 44980559 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1313 : Flx_BARRETT2_LIMIT))
1314 44997285 : return T;
1315 8020 : retmkvec2(Flx_invBarrett(T,p),T);
1316 : }
1317 :
1318 : /* separate from Flx_divrem for maximal speed. */
1319 : static GEN
1320 432367651 : Flx_rem_basecase(GEN x, GEN y, ulong p)
1321 : {
1322 : pari_sp av;
1323 : GEN z, c;
1324 : long dx,dy,dy1,dz,i,j;
1325 : ulong p1,inv;
1326 432367651 : long vs=x[1];
1327 :
1328 432367651 : dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1329 422887878 : dx = degpol(x);
1330 422643826 : dz = dx-dy; if (dz < 0) return Flx_copy(x);
1331 422643826 : x += 2; y += 2;
1332 422643826 : inv = y[dy];
1333 422643826 : if (inv != 1UL) inv = Fl_inv(inv,p);
1334 479021497 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1335 :
1336 423482734 : c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1337 422217380 : z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1338 :
1339 421254893 : if (SMALL_ULONG(p))
1340 : {
1341 270423972 : z[dz] = (inv*x[dx]) % p;
1342 1112754676 : for (i=dx-1; i>=dy; --i)
1343 : {
1344 842330704 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1345 7456565659 : for (j=i-dy1; j<=i && j<=dz; j++)
1346 : {
1347 6614234955 : p1 += z[j]*y[i-j];
1348 6614234955 : if (p1 & HIGHBIT) p1 %= p;
1349 : }
1350 842330704 : p1 %= p;
1351 842330704 : z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1352 : }
1353 2017416095 : for (i=0; i<dy; i++)
1354 : {
1355 1747072893 : p1 = z[0]*y[i];
1356 10052663459 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1357 : {
1358 8305590566 : p1 += z[j]*y[i-j];
1359 8305590566 : if (p1 & HIGHBIT) p1 %= p;
1360 : }
1361 1747095777 : c[i] = Fl_sub(x[i], p1%p, p);
1362 : }
1363 : }
1364 : else
1365 : {
1366 150830921 : ulong pi = get_Fl_red(p);
1367 150766690 : z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1368 541093872 : for (i=dx-1; i>=dy; --i)
1369 : {
1370 390121874 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1371 1753411594 : for (j=i-dy1; j<=i && j<=dz; j++)
1372 1362741836 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1373 390669758 : z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1374 : }
1375 1134008239 : for (i=0; i<dy; i++)
1376 : {
1377 983335742 : p1 = Fl_mul_pre(z[0],y[i],p,pi);
1378 3006126595 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1379 2019771446 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1380 972911412 : c[i] = Fl_sub(x[i], p1, p);
1381 : }
1382 : }
1383 513291744 : i = dy-1; while (i>=0 && !c[i]) i--;
1384 421015699 : set_avma(av);
1385 421602865 : return Flx_renormalize(c-2, i+3);
1386 : }
1387 :
1388 : /* as FpX_divrem but working only on ulong types.
1389 : * if relevant, *pr is the last object on stack */
1390 : static GEN
1391 27602024 : Flx_divrem_basecase(GEN x, GEN y, ulong p, GEN *pr)
1392 : {
1393 : GEN z,q,c;
1394 : long dx,dy,dy1,dz,i,j;
1395 : ulong p1,inv;
1396 27602024 : long sv=x[1];
1397 :
1398 27602024 : dy = degpol(y);
1399 27601471 : if (dy<0) pari_err_INV("Flx_divrem",y);
1400 27601518 : if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p);
1401 27601057 : if (!dy)
1402 : {
1403 4526998 : if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1404 4526981 : if (y[2] == 1UL) return Flx_copy(x);
1405 2933207 : return Flx_Fl_mul(x, Fl_inv(y[2], p), p);
1406 : }
1407 23074059 : dx = degpol(x);
1408 23077954 : dz = dx-dy;
1409 23077954 : if (dz < 0)
1410 : {
1411 331033 : q = pol0_Flx(sv);
1412 331026 : if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1413 331026 : return q;
1414 : }
1415 22746921 : x += 2;
1416 22746921 : y += 2;
1417 22746921 : z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1418 22746061 : inv = uel(y, dy);
1419 22746061 : if (inv != 1UL) inv = Fl_inv(inv,p);
1420 37917563 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1421 :
1422 22748101 : if (SMALL_ULONG(p))
1423 : {
1424 21637005 : z[dz] = (inv*x[dx]) % p;
1425 65107307 : for (i=dx-1; i>=dy; --i)
1426 : {
1427 43470302 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1428 156916692 : for (j=i-dy1; j<=i && j<=dz; j++)
1429 : {
1430 113446390 : p1 += z[j]*y[i-j];
1431 113446390 : if (p1 & HIGHBIT) p1 %= p;
1432 : }
1433 43470302 : p1 %= p;
1434 43470302 : z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1435 : }
1436 : }
1437 : else
1438 : {
1439 1111096 : z[dz] = Fl_mul(inv, x[dx], p);
1440 7139466 : for (i=dx-1; i>=dy; --i)
1441 : { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1442 6028377 : p1 = p - uel(x,i);
1443 21521224 : for (j=i-dy1; j<=i && j<=dz; j++)
1444 15492846 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1445 6028378 : z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1446 : }
1447 : }
1448 22748094 : q = Flx_renormalize(z-2, dz+3);
1449 22748215 : if (!pr) return q;
1450 :
1451 16807158 : c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1452 16807137 : if (SMALL_ULONG(p))
1453 : {
1454 176612212 : for (i=0; i<dy; i++)
1455 : {
1456 160832346 : p1 = (ulong)z[0]*y[i];
1457 393790308 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1458 : {
1459 232957962 : p1 += (ulong)z[j]*y[i-j];
1460 232957962 : if (p1 & HIGHBIT) p1 %= p;
1461 : }
1462 160831998 : c[i] = Fl_sub(x[i], p1%p, p);
1463 : }
1464 : }
1465 : else
1466 : {
1467 14001417 : for (i=0; i<dy; i++)
1468 : {
1469 12974567 : p1 = Fl_mul(z[0],y[i],p);
1470 47471087 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1471 34496522 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1472 12974572 : c[i] = Fl_sub(x[i], p1, p);
1473 : }
1474 : }
1475 22705395 : i=dy-1; while (i>=0 && !c[i]) i--;
1476 16806716 : c = Flx_renormalize(c-2, i+3);
1477 16807593 : if (pr == ONLY_DIVIDES)
1478 328 : { if (lg(c) != 2) return NULL; }
1479 : else
1480 16807265 : *pr = c;
1481 16807495 : return q;
1482 : }
1483 :
1484 : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1485 : * and mg is the Barrett inverse of T. */
1486 : static GEN
1487 707971 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, GEN *pr)
1488 : {
1489 : GEN q, r;
1490 707971 : long lt = degpol(T); /*We discard the leading term*/
1491 : long ld, lm, lT, lmg;
1492 707957 : ld = l-lt;
1493 707957 : lm = minss(ld, lgpol(mg));
1494 708186 : lT = Flx_lgrenormalizespec(T+2,lt);
1495 708302 : lmg = Flx_lgrenormalizespec(mg+2,lm);
1496 708182 : q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1497 707603 : q = Flx_mulspec(q+2,mg+2,p,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1498 708215 : q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1499 707558 : if (!pr) return q;
1500 701332 : r = Flx_mulspec(q+2,T+2,p,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1501 702056 : r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1502 701732 : if (pr == ONLY_REM) return r;
1503 375283 : *pr = r; return q;
1504 : }
1505 :
1506 : static GEN
1507 446814 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, GEN *pr)
1508 : {
1509 446814 : GEN q = NULL, r = Flx_copy(x);
1510 446858 : long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1511 : long i;
1512 446843 : if (l <= lt)
1513 : {
1514 0 : if (pr == ONLY_REM) return Flx_copy(x);
1515 0 : if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1516 0 : if (pr) *pr = Flx_copy(x);
1517 0 : return pol0_Flx(v);
1518 : }
1519 446843 : if (lt <= 1)
1520 701 : return Flx_divrem_basecase(x,T,p,pr);
1521 446142 : if (pr != ONLY_REM && l>lm)
1522 16395 : { q = zero_zv(l-lt+1); q[1] = T[1]; }
1523 709045 : while (l>lm)
1524 : {
1525 263046 : GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,&zr);
1526 262963 : long lz = lgpol(zr);
1527 262903 : if (pr != ONLY_REM)
1528 : {
1529 26327 : long lq = lgpol(zq);
1530 547118 : for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1531 : }
1532 3919235 : for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1533 262903 : l = l-lm+lz;
1534 : }
1535 445999 : if (pr == ONLY_REM)
1536 : {
1537 326482 : if (l > lt)
1538 326467 : r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,ONLY_REM);
1539 : else
1540 15 : r = Flx_renormalize(r, l+2);
1541 326463 : r[1] = v; return r;
1542 : }
1543 119517 : if (l > lt)
1544 : {
1545 118512 : GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p, pr ? &r: NULL);
1546 118512 : if (!q) q = zq;
1547 : else
1548 : {
1549 15256 : long lq = lgpol(zq);
1550 100940 : for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1551 : }
1552 : }
1553 1005 : else if (pr)
1554 1097 : r = Flx_renormalize(r, l+2);
1555 119517 : q[1] = v; q = Flx_renormalize(q, lg(q));
1556 119651 : if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1557 119651 : if (pr) { r[1] = v; *pr = r; }
1558 119651 : return q;
1559 : }
1560 :
1561 : GEN
1562 34607688 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1563 : {
1564 : GEN B, y;
1565 : long dy, dx, d;
1566 34607688 : if (pr==ONLY_REM) return Flx_rem(x, T, p);
1567 27721940 : y = get_Flx_red(T, &B);
1568 27724311 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1569 27720377 : if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1570 27600707 : return Flx_divrem_basecase(x,y,p,pr);
1571 : else
1572 : {
1573 119891 : pari_sp av = avma;
1574 119891 : GEN mg = B? B: Flx_invBarrett(y, p);
1575 119891 : GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pr);
1576 119891 : if (!q1) return gc_NULL(av);
1577 119891 : if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
1578 113521 : gerepileall(av,2,&q1,pr);
1579 113521 : return q1;
1580 : }
1581 : }
1582 :
1583 : GEN
1584 500057075 : Flx_rem(GEN x, GEN T, ulong p)
1585 : {
1586 500057075 : GEN B, y = get_Flx_red(T, &B);
1587 500197042 : long dy = degpol(y), dx = degpol(x), d = dx-dy;
1588 500116361 : if (d < 0) return Flx_copy(x);
1589 432923695 : if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1590 432497693 : return Flx_rem_basecase(x,y,p);
1591 : else
1592 : {
1593 326926 : pari_sp av=avma;
1594 326926 : GEN mg = B ? B: Flx_invBarrett(y, p);
1595 326928 : GEN r = Flx_divrem_Barrett(x, mg, y, p, ONLY_REM);
1596 326915 : return gerepileuptoleaf(av, r);
1597 : }
1598 : }
1599 :
1600 : /* reduce T mod (X^n - 1, p). Shallow function */
1601 : GEN
1602 4932767 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1603 : {
1604 4932767 : long i, j, L = lg(T), l = n+2;
1605 : GEN S;
1606 4932767 : if (L <= l || n & ~LGBITS) return T;
1607 455 : S = cgetg(l, t_VECSMALL);
1608 455 : S[1] = T[1];
1609 1736 : for (i = 2; i < l; i++) S[i] = T[i];
1610 1316 : for (j = 2; i < L; i++) {
1611 861 : S[j] = Fl_add(S[j], T[i], p);
1612 861 : if (++j == l) j = 2;
1613 : }
1614 455 : return Flx_renormalize(S, l);
1615 : }
1616 : /* reduce T mod (X^n + 1, p). Shallow function */
1617 : GEN
1618 5668 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
1619 : {
1620 5668 : long i, j, L = lg(T), l = n+2;
1621 : GEN S;
1622 5668 : if (L <= l || n & ~LGBITS) return T;
1623 406 : S = cgetg(l, t_VECSMALL);
1624 406 : S[1] = T[1];
1625 1547 : for (i = 2; i < l; i++) S[i] = T[i];
1626 1183 : for (j = 2; i < L; i++) {
1627 777 : S[j] = Fl_sub(S[j], T[i], p);
1628 777 : if (++j == l) j = 2;
1629 : }
1630 406 : return Flx_renormalize(S, l);
1631 : }
1632 :
1633 : struct _Flxq {
1634 : GEN aut;
1635 : GEN T;
1636 : ulong p;
1637 : };
1638 :
1639 : static GEN
1640 0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1641 : {
1642 0 : struct _Flxq *D = (struct _Flxq*) E;
1643 0 : return Flx_divrem(x, y, D->p, r);
1644 : }
1645 : static GEN
1646 504329 : _Flx_add(void * E, GEN x, GEN y) {
1647 504329 : struct _Flxq *D = (struct _Flxq*) E;
1648 504329 : return Flx_add(x, y, D->p);
1649 : }
1650 : static GEN
1651 9033001 : _Flx_mul(void *E, GEN x, GEN y) {
1652 9033001 : struct _Flxq *D = (struct _Flxq*) E;
1653 9033001 : return Flx_mul(x, y, D->p);
1654 : }
1655 : static GEN
1656 0 : _Flx_sqr(void *E, GEN x) {
1657 0 : struct _Flxq *D = (struct _Flxq*) E;
1658 0 : return Flx_sqr(x, D->p);
1659 : }
1660 :
1661 : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1662 :
1663 : GEN
1664 0 : Flx_digits(GEN x, GEN T, ulong p)
1665 : {
1666 0 : pari_sp av = avma;
1667 : struct _Flxq D;
1668 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
1669 : GEN z;
1670 0 : D.p = p;
1671 0 : z = gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1672 0 : return gerepileupto(av, z);
1673 : }
1674 :
1675 : GEN
1676 0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1677 : {
1678 0 : pari_sp av = avma;
1679 : struct _Flxq D;
1680 : GEN z;
1681 0 : D.p = p;
1682 0 : z = gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1683 0 : return gerepileupto(av, z);
1684 : }
1685 :
1686 : long
1687 3140856 : Flx_val(GEN x)
1688 : {
1689 3140856 : long i, l=lg(x);
1690 3140856 : if (l==2) return LONG_MAX;
1691 3146587 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1692 3140856 : return i-2;
1693 : }
1694 : long
1695 22061299 : Flx_valrem(GEN x, GEN *Z)
1696 : {
1697 22061299 : long v, i, l=lg(x);
1698 : GEN y;
1699 22061299 : if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1700 22160074 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1701 22061299 : v = i-2;
1702 22061299 : if (v == 0) { *Z = x; return 0; }
1703 34247 : l -= v;
1704 34247 : y = cgetg(l, t_VECSMALL); y[1] = x[1];
1705 129031 : for (i=2; i<l; i++) y[i] = x[i+v];
1706 39952 : *Z = y; return v;
1707 : }
1708 :
1709 : GEN
1710 6404394 : Flx_deriv(GEN z, ulong p)
1711 : {
1712 6404394 : long i,l = lg(z)-1;
1713 : GEN x;
1714 6404394 : if (l < 2) l = 2;
1715 6404394 : x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1716 6404266 : if (HIGHWORD(l | p))
1717 30012782 : for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1718 : else
1719 27370549 : for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1720 6404824 : return Flx_renormalize(x,l);
1721 : }
1722 :
1723 : static GEN
1724 8981 : Flx_integXn(GEN x, long n, ulong p)
1725 : {
1726 8981 : long i, lx = lg(x);
1727 : GEN y;
1728 8981 : if (lx == 2) return Flx_copy(x);
1729 7854 : y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1730 80683 : for (i=2; i<lx; i++)
1731 : {
1732 72794 : ulong xi = uel(x,i);
1733 72794 : if (xi == 0)
1734 6694 : uel(y,i) = 0;
1735 : else
1736 : {
1737 66100 : ulong j = n+i-1;
1738 66100 : ulong d = ugcd(j, xi);
1739 66081 : if (d==1)
1740 35804 : uel(y,i) = Fl_div(xi, j, p);
1741 : else
1742 30277 : uel(y,i) = Fl_div(xi/d, j/d, p);
1743 : }
1744 : }
1745 7889 : return Flx_renormalize(y, lx);;
1746 : }
1747 :
1748 : GEN
1749 0 : Flx_integ(GEN x, ulong p)
1750 : {
1751 0 : long i, lx = lg(x);
1752 : GEN y;
1753 0 : if (lx == 2) return Flx_copy(x);
1754 0 : y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1755 0 : uel(y,2) = 0;
1756 0 : for (i=3; i<=lx; i++)
1757 0 : uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1758 0 : return Flx_renormalize(y, lx+1);;
1759 : }
1760 :
1761 : /* assume p prime */
1762 : GEN
1763 11848 : Flx_diff1(GEN P, ulong p)
1764 : {
1765 11848 : return Flx_sub(Flx_translate1(P, p), P, p);
1766 : }
1767 :
1768 : GEN
1769 79699 : Flx_deflate(GEN x0, long d)
1770 : {
1771 : GEN z, y, x;
1772 79699 : long i,id, dy, dx = degpol(x0);
1773 79699 : if (d == 1 || dx <= 0) return Flx_copy(x0);
1774 72456 : dy = dx/d;
1775 72456 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1776 72456 : z = y + 2;
1777 72456 : x = x0+ 2;
1778 242017 : for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1779 72456 : return y;
1780 : }
1781 :
1782 : GEN
1783 48905 : Flx_inflate(GEN x0, long d)
1784 : {
1785 48905 : long i, id, dy, dx = degpol(x0);
1786 48905 : GEN x = x0 + 2, z, y;
1787 48905 : if (dx <= 0) return Flx_copy(x0);
1788 47769 : dy = dx*d;
1789 47769 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1790 47748 : z = y + 2;
1791 5368362 : for (i=0; i<=dy; i++) z[i] = 0;
1792 2696589 : for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1793 47748 : return y;
1794 : }
1795 :
1796 : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1797 : GEN
1798 137767 : Flx_splitting(GEN p, long k)
1799 : {
1800 137767 : long n = degpol(p), v = p[1], m, i, j, l;
1801 : GEN r;
1802 :
1803 137765 : m = n/k;
1804 137765 : r = cgetg(k+1,t_VEC);
1805 649998 : for(i=1; i<=k; i++)
1806 : {
1807 512249 : gel(r,i) = cgetg(m+3, t_VECSMALL);
1808 512252 : mael(r,i,1) = v;
1809 : }
1810 3233372 : for (j=1, i=0, l=2; i<=n; i++)
1811 : {
1812 3095623 : mael(r,j,l) = p[2+i];
1813 3095623 : if (j==k) { j=1; l++; } else j++;
1814 : }
1815 650019 : for(i=1; i<=k; i++)
1816 512297 : gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1817 137722 : return r;
1818 : }
1819 : static GEN
1820 277591 : Flx_halfgcd_basecase(GEN a, GEN b, ulong p)
1821 : {
1822 277591 : pari_sp av=avma;
1823 : GEN u,u1,v,v1;
1824 277591 : long vx = a[1];
1825 277591 : long n = lgpol(a)>>1;
1826 277586 : u1 = v = pol0_Flx(vx);
1827 277573 : u = v1 = pol1_Flx(vx);
1828 1203957 : while (lgpol(b)>n)
1829 : {
1830 926385 : GEN r, q = Flx_divrem(a,b,p, &r);
1831 926398 : a = b; b = r; swap(u,u1); swap(v,v1);
1832 926398 : u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
1833 926335 : v1 = Flx_sub(v1, Flx_mul(v, q ,p), p);
1834 926397 : if (gc_needed(av,2))
1835 : {
1836 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
1837 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
1838 : }
1839 : }
1840 277410 : return gerepilecopy(av, mkmat2(mkcol2(u,u1), mkcol2(v,v1)));
1841 : }
1842 : /* ux + vy */
1843 : static GEN
1844 15498 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p)
1845 15498 : { return Flx_add(Flx_mul(u,x, p), Flx_mul(v,y, p), p); }
1846 :
1847 : static GEN
1848 7746 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p)
1849 : {
1850 7746 : GEN res = cgetg(3, t_COL);
1851 7746 : gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p);
1852 7746 : gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p);
1853 7746 : return res;
1854 : }
1855 :
1856 : #if 0
1857 : static GEN
1858 : FlxM_mul2_old(GEN M, GEN N, ulong p)
1859 : {
1860 : GEN res = cgetg(3, t_MAT);
1861 : gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1862 : gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1863 : return res;
1864 : }
1865 : #endif
1866 : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1867 : static GEN
1868 1344 : FlxM_mul2(GEN A, GEN B, ulong p)
1869 : {
1870 1344 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1871 1344 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1872 1344 : GEN M1 = Flx_mul(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p);
1873 1344 : GEN M2 = Flx_mul(Flx_add(A21,A22, p), B11, p);
1874 1344 : GEN M3 = Flx_mul(A11, Flx_sub(B12,B22, p), p);
1875 1344 : GEN M4 = Flx_mul(A22, Flx_sub(B21,B11, p), p);
1876 1344 : GEN M5 = Flx_mul(Flx_add(A11,A12, p), B22, p);
1877 1344 : GEN M6 = Flx_mul(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p);
1878 1344 : GEN M7 = Flx_mul(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p);
1879 1344 : GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1880 1344 : GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1881 1344 : retmkmat2(mkcol2(Flx_add(T1,T2, p), Flx_add(M2,M4, p)),
1882 : mkcol2(Flx_add(M3,M5, p), Flx_add(T3,T4, p)));
1883 : }
1884 :
1885 : /* Return [0,1;1,-q]*M */
1886 : static GEN
1887 1341 : Flx_FlxM_qmul(GEN q, GEN M, ulong p)
1888 : {
1889 1341 : GEN u, v, res = cgetg(3, t_MAT);
1890 1341 : u = Flx_sub(gcoeff(M,1,1), Flx_mul(gcoeff(M,2,1), q, p), p);
1891 1341 : gel(res,1) = mkcol2(gcoeff(M,2,1), u);
1892 1341 : v = Flx_sub(gcoeff(M,1,2), Flx_mul(gcoeff(M,2,2), q, p), p);
1893 1341 : gel(res,2) = mkcol2(gcoeff(M,2,2), v);
1894 1341 : return res;
1895 : }
1896 :
1897 : static GEN
1898 3 : matid2_FlxM(long v)
1899 : {
1900 3 : return mkmat2(mkcol2(pol1_Flx(v),pol0_Flx(v)),
1901 : mkcol2(pol0_Flx(v),pol1_Flx(v)));
1902 : }
1903 :
1904 : static GEN
1905 7719 : Flx_halfgcd_split(GEN x, GEN y, ulong p)
1906 : {
1907 7719 : pari_sp av=avma;
1908 : GEN R, S, V;
1909 : GEN y1, r, q;
1910 7719 : long l = lgpol(x), n = l>>1, k;
1911 7719 : if (lgpol(y)<=n) return matid2_FlxM(x[1]);
1912 7719 : R = Flx_halfgcd(Flx_shift(x,-n),Flx_shift(y,-n),p);
1913 7719 : V = FlxM_Flx_mul2(R,x,y,p); y1 = gel(V,2);
1914 7719 : if (lgpol(y1)<=n) return gerepilecopy(av, R);
1915 1341 : q = Flx_divrem(gel(V,1), y1, p, &r);
1916 1341 : k = 2*n-degpol(y1);
1917 1341 : S = Flx_halfgcd(Flx_shift(y1,-k), Flx_shift(r,-k),p);
1918 1341 : return gerepileupto(av, FlxM_mul2(S,Flx_FlxM_qmul(q,R,p),p));
1919 : }
1920 :
1921 : /* Return M in GL_2(Fl[X]) such that:
1922 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
1923 : */
1924 :
1925 : static GEN
1926 285316 : Flx_halfgcd_i(GEN x, GEN y, ulong p)
1927 : {
1928 285316 : if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
1929 277590 : return Flx_halfgcd_basecase(x,y,p);
1930 7719 : return Flx_halfgcd_split(x,y,p);
1931 : }
1932 :
1933 : GEN
1934 285316 : Flx_halfgcd(GEN x, GEN y, ulong p)
1935 : {
1936 : pari_sp av;
1937 : GEN M,q,r;
1938 285316 : long lx=lgpol(x), ly=lgpol(y);
1939 285316 : if (!lx)
1940 : {
1941 0 : long v = x[1];
1942 0 : retmkmat2(mkcol2(pol0_Flx(v),pol1_Flx(v)),
1943 : mkcol2(pol1_Flx(v),pol0_Flx(v)));
1944 : }
1945 285316 : if (ly < lx) return Flx_halfgcd_i(x,y,p);
1946 2030 : av = avma;
1947 2030 : q = Flx_divrem(y,x,p,&r);
1948 2030 : M = Flx_halfgcd_i(x,r,p);
1949 2030 : gcoeff(M,1,1) = Flx_sub(gcoeff(M,1,1), Flx_mul(q, gcoeff(M,1,2), p), p);
1950 2030 : gcoeff(M,2,1) = Flx_sub(gcoeff(M,2,1), Flx_mul(q, gcoeff(M,2,2), p), p);
1951 2030 : return gerepilecopy(av, M);
1952 : }
1953 :
1954 : /*Do not garbage collect*/
1955 : static GEN
1956 35243899 : Flx_gcd_basecase(GEN a, GEN b, ulong p)
1957 : {
1958 35243899 : pari_sp av = avma;
1959 35243899 : ulong iter = 0;
1960 35243899 : if (lg(b) > lg(a)) swap(a, b);
1961 136988272 : while (lgpol(b))
1962 : {
1963 101508136 : GEN c = Flx_rem(a,b,p);
1964 101744373 : iter++; a = b; b = c;
1965 101744373 : if (gc_needed(av,2))
1966 : {
1967 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
1968 0 : gerepileall(av,2, &a,&b);
1969 : }
1970 : }
1971 35170235 : return iter < 2 ? Flx_copy(a) : a;
1972 : }
1973 :
1974 : GEN
1975 36030679 : Flx_gcd(GEN x, GEN y, ulong p)
1976 : {
1977 36030679 : pari_sp av = avma;
1978 : long lim;
1979 36030679 : if (!lgpol(x)) return Flx_copy(y);
1980 35248664 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
1981 35249683 : while (lgpol(y) >= lim)
1982 : {
1983 : GEN c;
1984 24 : if (lgpol(y)<=(lgpol(x)>>1))
1985 : {
1986 0 : GEN r = Flx_rem(x, y, p);
1987 0 : x = y; y = r;
1988 : }
1989 24 : c = FlxM_Flx_mul2(Flx_halfgcd(x,y, p), x, y, p);
1990 24 : x = gel(c,1); y = gel(c,2);
1991 24 : if (gc_needed(av,2))
1992 : {
1993 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
1994 0 : gerepileall(av,2,&x,&y);
1995 : }
1996 : }
1997 35240761 : return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p));
1998 : }
1999 :
2000 : int
2001 3754279 : Flx_is_squarefree(GEN z, ulong p)
2002 : {
2003 3754279 : pari_sp av = avma;
2004 3754279 : GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
2005 3754271 : return gc_bool(av, degpol(d) == 0);
2006 : }
2007 :
2008 : static long
2009 137521 : Flx_is_smooth_squarefree(GEN f, long r, ulong p)
2010 : {
2011 137521 : pari_sp av = avma;
2012 : long i;
2013 137521 : GEN sx = polx_Flx(f[1]), a = sx;
2014 137309 : for(i=1;;i++)
2015 : {
2016 588311 : if (degpol(f)<=r) return gc_long(av,1);
2017 562647 : a = Flxq_powu(Flx_rem(a,f,p), p, f, p);
2018 567069 : if (Flx_equal(a, sx)) return gc_long(av,1);
2019 563003 : if (i==r) return gc_long(av,0);
2020 452777 : f = Flx_div(f, Flx_gcd(Flx_sub(a,sx,p),f,p),p);
2021 : }
2022 : }
2023 :
2024 : static long
2025 8953 : Flx_is_l_pow(GEN x, ulong p)
2026 : {
2027 8953 : ulong i, lx = lgpol(x);
2028 18025 : for (i=1; i<lx; i++)
2029 16107 : if (x[i+2] && i%p) return 0;
2030 1918 : return 1;
2031 : }
2032 :
2033 : int
2034 137549 : Flx_is_smooth(GEN g, long r, ulong p)
2035 : {
2036 : while (1)
2037 8958 : {
2038 137549 : GEN f = Flx_gcd(g, Flx_deriv(g, p), p);
2039 137427 : if (!Flx_is_smooth_squarefree(Flx_div(g, f, p), r, p))
2040 110197 : return 0;
2041 27479 : if (degpol(f)==0) return 1;
2042 8932 : g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2043 : }
2044 : }
2045 :
2046 : static GEN
2047 4073864 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
2048 : {
2049 4073864 : pari_sp av=avma;
2050 : GEN u,v,d,d1,v1;
2051 4073864 : long vx = a[1];
2052 4073864 : d = a; d1 = b;
2053 4073864 : v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2054 20931173 : while (lgpol(d1))
2055 : {
2056 16857528 : GEN r, q = Flx_divrem(d,d1,p, &r);
2057 16857667 : v = Flx_sub(v,Flx_mul(q,v1,p),p);
2058 16857420 : u=v; v=v1; v1=u;
2059 16857420 : u=r; d=d1; d1=u;
2060 16857420 : if (gc_needed(av,2))
2061 : {
2062 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(d));
2063 0 : gerepileall(av,5, &d,&d1,&u,&v,&v1);
2064 : }
2065 : }
2066 4073123 : if (ptu) *ptu = Flx_div(Flx_sub(d, Flx_mul(b,v,p), p), a, p);
2067 4073740 : *ptv = v; return d;
2068 : }
2069 :
2070 : static GEN
2071 3 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2072 : {
2073 3 : pari_sp av=avma;
2074 3 : GEN u,v,R = matid2_FlxM(x[1]);
2075 3 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2076 6 : while (lgpol(y) >= lim)
2077 : {
2078 : GEN M, c;
2079 3 : if (lgpol(y)<=(lgpol(x)>>1))
2080 : {
2081 0 : GEN r, q = Flx_divrem(x, y, p, &r);
2082 0 : x = y; y = r;
2083 0 : R = Flx_FlxM_qmul(q, R, p);
2084 : }
2085 3 : M = Flx_halfgcd(x,y, p);
2086 3 : c = FlxM_Flx_mul2(M, x,y, p);
2087 3 : R = FlxM_mul2(M, R, p);
2088 3 : x = gel(c,1); y = gel(c,2);
2089 3 : gerepileall(av,3,&x,&y,&R);
2090 : }
2091 3 : y = Flx_extgcd_basecase(x,y,p,&u,&v);
2092 3 : if (ptu) *ptu = Flx_addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1),p);
2093 3 : *ptv = Flx_addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2),p);
2094 3 : return y;
2095 : }
2096 :
2097 : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2098 : * ux + vy = gcd (mod p) */
2099 : GEN
2100 4073841 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2101 : {
2102 : GEN d;
2103 4073841 : pari_sp ltop=avma;
2104 4073841 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2105 4073857 : if (lgpol(y) >= lim)
2106 3 : d = Flx_extgcd_halfgcd(x, y, p, ptu, ptv);
2107 : else
2108 4073847 : d = Flx_extgcd_basecase(x, y, p, ptu, ptv);
2109 4073746 : gerepileall(ltop,ptu?3:2,&d,ptv,ptu);
2110 4074020 : return d;
2111 : }
2112 :
2113 : ulong
2114 2089082 : Flx_resultant(GEN a, GEN b, ulong p)
2115 : {
2116 : long da,db,dc,cnt;
2117 2089082 : ulong lb, res = 1UL;
2118 : pari_sp av;
2119 : GEN c;
2120 :
2121 2089082 : if (lgpol(a)==0 || lgpol(b)==0) return 0;
2122 2087750 : da = degpol(a);
2123 2087759 : db = degpol(b);
2124 2087787 : if (db > da)
2125 : {
2126 100387 : swapspec(a,b, da,db);
2127 100387 : if (both_odd(da,db)) res = p-res;
2128 : }
2129 1987400 : else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2130 2087786 : cnt = 0; av = avma;
2131 30526909 : while (db)
2132 : {
2133 28448654 : lb = b[db+2];
2134 28448654 : c = Flx_rem(a,b, p);
2135 28360699 : a = b; b = c; dc = degpol(c);
2136 28357552 : if (dc < 0) return gc_long(av,0);
2137 :
2138 28357324 : if (both_odd(da,db)) res = p - res;
2139 28364637 : if (lb != 1) res = Fl_mul(res, Fl_powu(lb, da - dc, p), p);
2140 28418482 : if (++cnt == 100) { cnt = 0; gerepileall(av, 2, &a, &b); }
2141 28439123 : da = db; /* = degpol(a) */
2142 28439123 : db = dc; /* = degpol(b) */
2143 : }
2144 2078255 : return gc_ulong(av, Fl_mul(res, Fl_powu(b[2], da, p), p));
2145 : }
2146 :
2147 : /* If resultant is 0, *ptU and *ptV are not set */
2148 : ulong
2149 0 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2150 : {
2151 0 : GEN z,q,u,v, x = a, y = b;
2152 0 : ulong lb, res = 1UL;
2153 0 : pari_sp av = avma;
2154 : long dx, dy, dz;
2155 0 : long vs=a[1];
2156 :
2157 0 : dx = degpol(x);
2158 0 : dy = degpol(y);
2159 0 : if (dy > dx)
2160 : {
2161 0 : swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
2162 0 : a = x; b = y;
2163 0 : if (both_odd(dx,dy)) res = p-res;
2164 : }
2165 : /* dy <= dx */
2166 0 : if (dy < 0) return 0;
2167 :
2168 0 : u = pol0_Flx(vs);
2169 0 : v = pol1_Flx(vs); /* v = 1 */
2170 0 : while (dy)
2171 : { /* b u = x (a), b v = y (a) */
2172 0 : lb = y[dy+2];
2173 0 : q = Flx_divrem(x,y, p, &z);
2174 0 : x = y; y = z; /* (x,y) = (y, x - q y) */
2175 0 : dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2176 0 : z = Flx_sub(u, Flx_mul(q,v, p), p);
2177 0 : u = v; v = z; /* (u,v) = (v, u - q v) */
2178 :
2179 0 : if (both_odd(dx,dy)) res = p - res;
2180 0 : if (lb != 1) res = Fl_mul(res, Fl_powu(lb, dx-dz, p), p);
2181 0 : dx = dy; /* = degpol(x) */
2182 0 : dy = dz; /* = degpol(y) */
2183 : }
2184 0 : res = Fl_mul(res, Fl_powu(y[2], dx, p), p);
2185 0 : lb = Fl_mul(res, Fl_inv(y[2],p), p);
2186 0 : v = gerepileuptoleaf(av, Flx_Fl_mul(v, lb, p));
2187 0 : av = avma;
2188 0 : u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul(b,v,p), p);
2189 0 : u = gerepileuptoleaf(av, Flx_div(u,a,p)); /* = (res - b v) / a */
2190 0 : *ptU = u;
2191 0 : *ptV = v; return res;
2192 : }
2193 :
2194 : ulong
2195 31463013 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2196 : {
2197 31463013 : ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2198 : LOCAL_OVERFLOW;
2199 : LOCAL_HIREMAINDER;
2200 31463013 : x++;
2201 :
2202 31463013 : if (lx == 1)
2203 3479910 : return 0;
2204 27983103 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2205 64560027 : while (++i < lx) {
2206 36576924 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2207 36576924 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2208 : }
2209 27983103 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
2210 60798 : else return remlll_pre(v1, h1, l1, p, pi);
2211 : }
2212 :
2213 : INLINE ulong
2214 3975776 : Flx_eval_pre_i(GEN x, ulong y, ulong p, ulong pi)
2215 : {
2216 : ulong p1;
2217 3975776 : long i=lg(x)-1;
2218 3975776 : if (i<=2)
2219 2341521 : return (i==2)? x[2]: 0;
2220 1634255 : p1 = x[i];
2221 6495147 : for (i--; i>=2; i--)
2222 4860155 : p1 = Fl_addmul_pre(uel(x, i), p1, y, p, pi);
2223 1634992 : return p1;
2224 : }
2225 :
2226 : ulong
2227 4055610 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2228 : {
2229 4055610 : if (degpol(x) > 15)
2230 : {
2231 79871 : pari_sp av = avma;
2232 79871 : GEN v = Fl_powers_pre(y, degpol(x), p, pi);
2233 79897 : ulong r = Flx_eval_powers_pre(x, v, p, pi);
2234 79907 : return gc_ulong(av,r);
2235 : }
2236 : else
2237 3975617 : return Flx_eval_pre_i(x, y, p, pi);
2238 : }
2239 :
2240 : ulong
2241 4051204 : Flx_eval(GEN x, ulong y, ulong p)
2242 : {
2243 4051204 : return Flx_eval_pre(x, y, p, get_Fl_red(p));
2244 : }
2245 :
2246 : ulong
2247 3073 : Flv_prod_pre(GEN x, ulong p, ulong pi)
2248 : {
2249 3073 : pari_sp ltop = avma;
2250 : GEN v;
2251 3073 : long i,k,lx = lg(x);
2252 3073 : if (lx == 1) return 1UL;
2253 3073 : if (lx == 2) return uel(x,1);
2254 2884 : v = cgetg(1+(lx << 1), t_VECSMALL);
2255 2884 : k = 1;
2256 27244 : for (i=1; i<lx-1; i+=2)
2257 24360 : uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2258 2884 : if (i < lx) uel(v,k++) = uel(x,i);
2259 12964 : while (k > 2)
2260 : {
2261 10080 : lx = k; k = 1;
2262 34440 : for (i=1; i<lx-1; i+=2)
2263 24360 : uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2264 10080 : if (i < lx) uel(v,k++) = uel(v,i);
2265 : }
2266 2884 : return gc_ulong(ltop, uel(v,1));
2267 : }
2268 :
2269 : ulong
2270 0 : Flv_prod(GEN v, ulong p)
2271 : {
2272 0 : return Flv_prod_pre(v, p, get_Fl_red(p));
2273 : }
2274 :
2275 : GEN
2276 0 : FlxV_prod(GEN V, ulong p)
2277 : {
2278 : struct _Flxq D;
2279 0 : D.T = NULL; D.aut = NULL; D.p = p;
2280 0 : return gen_product(V, (void *)&D, &_Flx_mul);
2281 : }
2282 :
2283 : /* compute prod (x - a[i]) */
2284 : GEN
2285 615741 : Flv_roots_to_pol(GEN a, ulong p, long vs)
2286 : {
2287 : struct _Flxq D;
2288 615741 : long i,k,lx = lg(a);
2289 : GEN p1;
2290 615741 : if (lx == 1) return pol1_Flx(vs);
2291 615741 : p1 = cgetg(lx, t_VEC);
2292 10201494 : for (k=1,i=1; i<lx-1; i+=2)
2293 9585674 : gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2294 9586320 : Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2295 615174 : if (i < lx)
2296 52421 : gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2297 615174 : D.T = NULL; D.aut = NULL; D.p = p;
2298 615174 : setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2299 : }
2300 :
2301 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2302 : INLINE void
2303 10333735 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2304 : {
2305 10333735 : pari_sp av = avma;
2306 10333735 : long n = lg(w), i;
2307 : ulong u;
2308 : GEN c;
2309 :
2310 10333735 : if (n == 1) return;
2311 10333735 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2312 56529444 : for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2313 10333740 : i = n-1; u = Fl_inv(c[i], p);
2314 56529417 : for ( ; i > 1; --i)
2315 : {
2316 46195682 : ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2317 46195687 : u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2318 : }
2319 10333735 : v[1] = u; set_avma(av);
2320 : }
2321 :
2322 : void
2323 10294276 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2324 :
2325 : GEN
2326 10702 : Flv_inv_pre(GEN w, ulong p, ulong pi)
2327 10702 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2328 :
2329 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2330 : INLINE void
2331 32121 : Flv_inv_indir(GEN w, GEN v, ulong p)
2332 : {
2333 32121 : pari_sp av = avma;
2334 32121 : long n = lg(w), i;
2335 : ulong u;
2336 : GEN c;
2337 :
2338 32121 : if (n == 1) return;
2339 32121 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2340 418710 : for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2341 32119 : i = n-1; u = Fl_inv(c[i], p);
2342 418762 : for ( ; i > 1; --i)
2343 : {
2344 386633 : ulong t = Fl_mul(u, c[i-1], p);
2345 386629 : u = Fl_mul(u, w[i], p); v[i] = t;
2346 : }
2347 32129 : v[1] = u; set_avma(av);
2348 : }
2349 : static void
2350 60878 : Flv_inv_i(GEN v, GEN w, ulong p)
2351 : {
2352 60878 : if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2353 28757 : else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2354 60885 : }
2355 : void
2356 221 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2357 : GEN
2358 60664 : Flv_inv(GEN w, ulong p)
2359 60664 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2360 :
2361 : GEN
2362 28601543 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2363 : {
2364 28601543 : long l = lg(a), i;
2365 : GEN a0, z0, z;
2366 28601543 : if (l <= 3)
2367 : {
2368 0 : if (rem) *rem = l == 2? 0: a[2];
2369 0 : return zero_Flx(a[1]);
2370 : }
2371 28601543 : z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2372 28456600 : a0 = a + l-1;
2373 28456600 : z0 = z + l-2; *z0 = *a0--;
2374 28456600 : if (SMALL_ULONG(p))
2375 : {
2376 69725685 : for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2377 : {
2378 52302174 : ulong t = (*a0-- + x * *z0--) % p;
2379 52302174 : *z0 = (long)t;
2380 : }
2381 17423511 : if (rem) *rem = (*a0 + x * *z0) % p;
2382 : }
2383 : else
2384 : {
2385 43348296 : for (i=l-3; i>1; i--)
2386 : {
2387 32262943 : ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2388 32315207 : *z0 = (long)t;
2389 : }
2390 11085353 : if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2391 : }
2392 28537142 : return z;
2393 : }
2394 :
2395 : /* xa, ya = t_VECSMALL */
2396 : static GEN
2397 58930 : Flv_producttree(GEN xa, GEN s, ulong p, long vs)
2398 : {
2399 58930 : long n = lg(xa)-1;
2400 58930 : long m = n==1 ? 1: expu(n-1)+1;
2401 58930 : long i, j, k, ls = lg(s);
2402 58930 : GEN T = cgetg(m+1, t_VEC);
2403 58921 : GEN t = cgetg(ls, t_VEC);
2404 689262 : for (j=1, k=1; j<ls; k+=s[j++])
2405 630263 : gel(t, j) = s[j] == 1 ?
2406 630344 : mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2407 215881 : mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2408 215880 : Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2409 58918 : gel(T,1) = t;
2410 221201 : for (i=2; i<=m; i++)
2411 : {
2412 162306 : GEN u = gel(T, i-1);
2413 162306 : long n = lg(u)-1;
2414 162306 : GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2415 733615 : for (j=1, k=1; k<n; j++, k+=2)
2416 571332 : gel(t, j) = Flx_mul(gel(u, k), gel(u, k+1), p);
2417 162283 : gel(T, i) = t;
2418 : }
2419 58895 : return T;
2420 : }
2421 :
2422 : static GEN
2423 58928 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p)
2424 : {
2425 : long i,j,k;
2426 58928 : long m = lg(T)-1;
2427 : GEN t;
2428 58928 : GEN R = cgetg(lg(xa), t_VECSMALL);
2429 58924 : GEN Tp = cgetg(m+1, t_VEC);
2430 58915 : gel(Tp, m) = mkvec(P);
2431 221212 : for (i=m-1; i>=1; i--)
2432 : {
2433 162289 : GEN u = gel(T, i);
2434 162289 : GEN v = gel(Tp, i+1);
2435 162289 : long n = lg(u)-1;
2436 162289 : t = cgetg(n+1, t_VEC);
2437 733477 : for (j=1, k=1; k<n; j++, k+=2)
2438 : {
2439 571188 : gel(t, k) = Flx_rem(gel(v, j), gel(u, k), p);
2440 571213 : gel(t, k+1) = Flx_rem(gel(v, j), gel(u, k+1), p);
2441 : }
2442 162289 : gel(Tp, i) = t;
2443 : }
2444 : {
2445 58923 : GEN u = gel(T, i+1);
2446 58923 : GEN v = gel(Tp, i+1);
2447 58923 : long n = lg(u)-1;
2448 689479 : for (j=1, k=1; j<=n; j++)
2449 : {
2450 630515 : long c, d = degpol(gel(u,j));
2451 1476935 : for (c=1; c<=d; c++, k++)
2452 846379 : R[k] = Flx_eval(gel(v, j), xa[k], p);
2453 : }
2454 58964 : set_avma((pari_sp)R); return R;
2455 : }
2456 : }
2457 :
2458 : static GEN
2459 752313 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, long vs)
2460 : {
2461 752313 : pari_sp av = avma;
2462 752313 : long m = lg(T)-1;
2463 752313 : long i, j, k, ls = lg(s);
2464 752313 : GEN Tp = cgetg(m+1, t_VEC);
2465 751510 : GEN t = cgetg(ls, t_VEC);
2466 13085815 : for (j=1, k=1; j<ls; k+=s[j++])
2467 12334298 : if (s[j]==2)
2468 : {
2469 4216317 : ulong a = Fl_mul(ya[k], R[k], p);
2470 4217383 : ulong b = Fl_mul(ya[k+1], R[k+1], p);
2471 4225236 : gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2472 4221296 : Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2473 4218131 : gel(t, j) = Flx_renormalize(gel(t, j), 4);
2474 : }
2475 : else
2476 8117981 : gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2477 751517 : gel(Tp, 1) = t;
2478 3336643 : for (i=2; i<=m; i++)
2479 : {
2480 2585457 : GEN u = gel(T, i-1);
2481 2585457 : GEN t = cgetg(lg(gel(T,i)), t_VEC);
2482 2582324 : GEN v = gel(Tp, i-1);
2483 2582324 : long n = lg(v)-1;
2484 14110544 : for (j=1, k=1; k<n; j++, k+=2)
2485 11542263 : gel(t, j) = Flx_add(Flx_mul(gel(u, k), gel(v, k+1), p),
2486 11525418 : Flx_mul(gel(u, k+1), gel(v, k), p), p);
2487 2585126 : gel(Tp, i) = t;
2488 : }
2489 751186 : return gerepileuptoleaf(av, gmael(Tp,m,1));
2490 : }
2491 :
2492 : GEN
2493 0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2494 : {
2495 0 : pari_sp av = avma;
2496 0 : GEN s = producttree_scheme(lg(xa)-1);
2497 0 : GEN T = Flv_producttree(xa, s, p, P[1]);
2498 0 : return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p));
2499 : }
2500 :
2501 : static GEN
2502 0 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p)
2503 0 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p)) }
2504 :
2505 : GEN
2506 0 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2507 : {
2508 0 : pari_sp av = avma;
2509 0 : GEN s = producttree_scheme(lg(xa)-1);
2510 0 : GEN T = Flv_producttree(xa, s, p, P[1]);
2511 0 : return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p));
2512 : }
2513 :
2514 : GEN
2515 15776 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2516 : {
2517 15776 : pari_sp av = avma;
2518 15776 : GEN s = producttree_scheme(lg(xa)-1);
2519 15776 : GEN T = Flv_producttree(xa, s, p, vs);
2520 15776 : long m = lg(T)-1;
2521 15776 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2522 15777 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2523 15776 : return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, vs));
2524 : }
2525 :
2526 : GEN
2527 38430 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2528 : {
2529 38430 : pari_sp av = avma;
2530 38430 : GEN s = producttree_scheme(lg(xa)-1);
2531 38430 : GEN T = Flv_producttree(xa, s, p, vs);
2532 38421 : long i, m = lg(T)-1, l = lg(ya)-1;
2533 38421 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2534 38425 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2535 38426 : GEN M = cgetg(l+1, t_VEC);
2536 774618 : for (i=1; i<=l; i++)
2537 736210 : gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, vs);
2538 38408 : return gerepileupto(av, M);
2539 : }
2540 :
2541 : GEN
2542 4726 : Flv_invVandermonde(GEN L, ulong den, ulong p)
2543 : {
2544 4726 : pari_sp av = avma;
2545 4726 : long i, n = lg(L);
2546 : GEN M, R;
2547 4726 : GEN s = producttree_scheme(n-1);
2548 4726 : GEN tree = Flv_producttree(L, s, p, 0);
2549 4726 : long m = lg(tree)-1;
2550 4726 : GEN T = gmael(tree, m, 1);
2551 4726 : R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p), p);
2552 4726 : if (den!=1) R = Flv_Fl_mul(R, den, p);
2553 4726 : M = cgetg(n, t_MAT);
2554 19874 : for (i = 1; i < n; i++)
2555 : {
2556 15148 : GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2557 15148 : gel(M,i) = Flx_to_Flv(P, n-1);
2558 : }
2559 4726 : return gerepilecopy(av, M);
2560 : }
2561 :
2562 : /***********************************************************************/
2563 : /** **/
2564 : /** Flxq **/
2565 : /** **/
2566 : /***********************************************************************/
2567 : /* Flxq objects are defined as follows:
2568 : They are Flx modulo another Flx called q.
2569 : */
2570 :
2571 : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2572 : GEN
2573 121017678 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2574 : {
2575 121017678 : return Flx_rem(Flx_mul(x,y,p),T,p);
2576 : }
2577 :
2578 : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2579 : GEN
2580 185081539 : Flxq_sqr(GEN x,GEN T,ulong p)
2581 : {
2582 185081539 : return Flx_rem(Flx_sqr(x,p),T,p);
2583 : }
2584 :
2585 : static GEN
2586 1658644 : _Flxq_red(void *E, GEN x)
2587 1658644 : { struct _Flxq *s = (struct _Flxq *)E;
2588 1658644 : return Flx_rem(x, s->T, s->p); }
2589 : #if 0
2590 : static GEN
2591 : _Flx_sub(void *E, GEN x, GEN y)
2592 : { struct _Flxq *s = (struct _Flxq *)E;
2593 : return Flx_sub(x,y,s->p); }
2594 : #endif
2595 : static GEN
2596 179244697 : _Flxq_sqr(void *data, GEN x)
2597 : {
2598 179244697 : struct _Flxq *D = (struct _Flxq*)data;
2599 179244697 : return Flxq_sqr(x, D->T, D->p);
2600 : }
2601 : static GEN
2602 92729714 : _Flxq_mul(void *data, GEN x, GEN y)
2603 : {
2604 92729714 : struct _Flxq *D = (struct _Flxq*)data;
2605 92729714 : return Flxq_mul(x,y, D->T, D->p);
2606 : }
2607 : static GEN
2608 5374021 : _Flxq_one(void *data)
2609 : {
2610 5374021 : struct _Flxq *D = (struct _Flxq*)data;
2611 5374021 : return pol1_Flx(get_Flx_var(D->T));
2612 : }
2613 : #if 0
2614 : static GEN
2615 : _Flxq_zero(void *data)
2616 : {
2617 : struct _Flxq *D = (struct _Flxq*)data;
2618 : return pol0_Flx(get_Flx_var(D->T));
2619 : }
2620 : static GEN
2621 : _Flxq_cmul(void *data, GEN P, long a, GEN x)
2622 : {
2623 : struct _Flxq *D = (struct _Flxq*)data;
2624 : return Flx_Fl_mul(x, P[a+2], D->p);
2625 : }
2626 : #endif
2627 :
2628 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2629 : GEN
2630 11963111 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2631 : {
2632 11963111 : pari_sp av = avma;
2633 : struct _Flxq D;
2634 : GEN y;
2635 11963111 : switch(n)
2636 : {
2637 0 : case 0: return pol1_Flx(get_Flx_var(T));
2638 62679 : case 1: return Flx_copy(x);
2639 189770 : case 2: return Flxq_sqr(x, T, p);
2640 : }
2641 11710662 : D.T = Flx_get_red(T, p); D.p = p;
2642 11710123 : y = gen_powu_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2643 11690543 : return gerepileuptoleaf(av, y);
2644 : }
2645 :
2646 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2647 : GEN
2648 16258041 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
2649 : {
2650 16258041 : pari_sp av = avma;
2651 : struct _Flxq D;
2652 : GEN y;
2653 16258041 : long s = signe(n);
2654 16258041 : if (!s) return pol1_Flx(get_Flx_var(T));
2655 16082950 : if (s < 0)
2656 566538 : x = Flxq_inv(x,T,p);
2657 16082950 : if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2658 15367867 : D.T = Flx_get_red(T, p); D.p = p;
2659 15367865 : y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2660 15367723 : return gerepileuptoleaf(av, y);
2661 : }
2662 :
2663 : GEN
2664 28 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
2665 : {
2666 : struct _Flxq D;
2667 28 : D.T = Flx_get_red(T, p); D.p = p;
2668 28 : return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2669 : }
2670 :
2671 : GEN
2672 4397 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
2673 : {
2674 : struct _Flxq D;
2675 4397 : D.T = Flx_get_red(T, p); D.p = p;
2676 4397 : return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
2677 : }
2678 :
2679 : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
2680 : * not stack clean.
2681 : */
2682 : GEN
2683 3920912 : Flxq_invsafe(GEN x, GEN T, ulong p)
2684 : {
2685 3920912 : GEN V, z = Flx_extgcd(get_Flx_mod(T), x, p, NULL, &V);
2686 : ulong iz;
2687 3921078 : if (degpol(z)) return NULL;
2688 3920405 : iz = Fl_inv (uel(z,2), p);
2689 3920365 : return Flx_Fl_mul(V, iz, p);
2690 : }
2691 :
2692 : GEN
2693 3690987 : Flxq_inv(GEN x,GEN T,ulong p)
2694 : {
2695 3690987 : pari_sp av=avma;
2696 3690987 : GEN U = Flxq_invsafe(x, T, p);
2697 3690987 : if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
2698 3690959 : return gerepileuptoleaf(av, U);
2699 : }
2700 :
2701 : GEN
2702 1906072 : Flxq_div(GEN x,GEN y,GEN T,ulong p)
2703 : {
2704 1906072 : pari_sp av = avma;
2705 1906072 : return gerepileuptoleaf(av, Flxq_mul(x,Flxq_inv(y,T,p),T,p));
2706 : }
2707 :
2708 : GEN
2709 5369737 : Flxq_powers(GEN x, long l, GEN T, ulong p)
2710 : {
2711 : struct _Flxq D;
2712 5369737 : int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
2713 5369633 : D.T = Flx_get_red(T, p); D.p = p;
2714 5369548 : return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
2715 : }
2716 :
2717 : GEN
2718 34905 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
2719 : {
2720 34905 : return FlxV_to_Flm(Flxq_powers(y,m-1,P,l),n);
2721 : }
2722 :
2723 : GEN
2724 4279222 : Flx_Frobenius(GEN T, ulong p)
2725 : {
2726 4279222 : return Flxq_powu(polx_Flx(get_Flx_var(T)), p, T, p);
2727 : }
2728 :
2729 : GEN
2730 17978 : Flx_matFrobenius(GEN T, ulong p)
2731 : {
2732 17978 : long n = get_Flx_degree(T);
2733 17978 : return Flxq_matrix_pow(Flx_Frobenius(T, p), n, n, T, p);
2734 : }
2735 :
2736 : static GEN
2737 5071753 : Flx_blocks_Flm(GEN P, long n, long m)
2738 : {
2739 5071753 : GEN z = cgetg(m+1,t_MAT);
2740 5071459 : long i,j, k=2, l = lg(P);
2741 15829808 : for(i=1; i<=m; i++)
2742 : {
2743 10760135 : GEN zi = cgetg(n+1,t_VECSMALL);
2744 10758349 : gel(z,i) = zi;
2745 50895994 : for(j=1; j<=n; j++)
2746 40137645 : uel(zi, j) = k==l ? 0 : uel(P,k++);
2747 : }
2748 5069673 : return z;
2749 : }
2750 :
2751 : GEN
2752 21654 : Flx_blocks(GEN P, long n, long m)
2753 : {
2754 21654 : GEN z = cgetg(m+1,t_VEC);
2755 21656 : long i,j, k=2, l = lg(P);
2756 64956 : for(i=1; i<=m; i++)
2757 : {
2758 43299 : GEN zi = cgetg(n+2,t_VECSMALL);
2759 43299 : zi[1] = P[1];
2760 43299 : gel(z,i) = zi;
2761 380925 : for(j=2; j<n+2; j++)
2762 337626 : uel(zi, j) = k==l ? 0 : uel(P,k++);
2763 43299 : zi = Flx_renormalize(zi, n+2);
2764 : }
2765 21657 : return z;
2766 : }
2767 :
2768 : static GEN
2769 5071824 : FlxV_to_Flm_lg(GEN x, long m, long n)
2770 : {
2771 : long i;
2772 5071824 : GEN y = cgetg(n+1, t_MAT);
2773 24692257 : for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
2774 5070995 : return y;
2775 : }
2776 :
2777 : GEN
2778 5269542 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
2779 : {
2780 5269542 : pari_sp btop, av = avma;
2781 5269542 : long sv = get_Flx_var(T), m = get_Flx_degree(T);
2782 5269605 : long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
2783 : GEN A, B, C, S, g;
2784 5269969 : if (lQ == 0) return pol0_Flx(sv);
2785 5072137 : if (lQ <= l)
2786 : {
2787 2219959 : n = l;
2788 2219959 : d = 1;
2789 : }
2790 : else
2791 : {
2792 2852178 : n = l-1;
2793 2852178 : d = (lQ+n-1)/n;
2794 : }
2795 5072137 : A = FlxV_to_Flm_lg(x, m, n);
2796 5071651 : B = Flx_blocks_Flm(Q, n, d);
2797 5070833 : C = gerepileupto(av, Flm_mul(A, B, p));
2798 5072475 : g = gel(x, l);
2799 5072475 : T = Flx_get_red(T, p);
2800 5072276 : btop = avma;
2801 5072276 : S = Flv_to_Flx(gel(C, d), sv);
2802 10762374 : for (i = d-1; i>0; i--)
2803 : {
2804 5691116 : S = Flx_add(Flxq_mul(S, g, T, p), Flv_to_Flx(gel(C,i), sv), p);
2805 5689865 : if (gc_needed(btop,1))
2806 0 : S = gerepileuptoleaf(btop, S);
2807 : }
2808 5071258 : return gerepileuptoleaf(av, S);
2809 : }
2810 :
2811 : GEN
2812 1153729 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
2813 : {
2814 1153729 : pari_sp av = avma;
2815 : GEN z, V;
2816 1153729 : long d = degpol(Q), rtd;
2817 1153725 : if (d < 0) return pol0_Flx(get_Flx_var(T));
2818 1153634 : rtd = (long) sqrt((double)d);
2819 1153634 : T = Flx_get_red(T, p);
2820 1153648 : V = Flxq_powers(x, rtd, T, p);
2821 1153634 : z = Flx_FlxqV_eval(Q, V, T, p);
2822 1153635 : return gerepileupto(av, z);
2823 : }
2824 :
2825 : GEN
2826 0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
2827 0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval(gel(x,i), v, T, p))
2828 : }
2829 :
2830 : GEN
2831 0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
2832 : {
2833 0 : long d = brent_kung_optpow(degpol(T)-1,lg(x)-1,1);
2834 0 : GEN Fp = Flxq_powers(F, d, T, p);
2835 0 : return FlxC_FlxqV_eval(x, Fp, T, p);
2836 : }
2837 :
2838 : #if 0
2839 : static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
2840 : _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
2841 : #endif
2842 :
2843 : static GEN
2844 368573 : Flxq_autpow_sqr(void *E, GEN x)
2845 : {
2846 368573 : struct _Flxq *D = (struct _Flxq*)E;
2847 368573 : return Flx_Flxq_eval(x, x, D->T, D->p);
2848 : }
2849 : static GEN
2850 20223 : Flxq_autpow_msqr(void *E, GEN x)
2851 : {
2852 20223 : struct _Flxq *D = (struct _Flxq*)E;
2853 20223 : return Flx_FlxqV_eval(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p);
2854 : }
2855 :
2856 : GEN
2857 295553 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
2858 : {
2859 295553 : pari_sp av = avma;
2860 : struct _Flxq D;
2861 : long d;
2862 295553 : if (n==0) return Flx_rem(polx_Flx(x[1]), T, p);
2863 295546 : if (n==1) return Flx_rem(x, T, p);
2864 295049 : D.T = Flx_get_red(T, p); D.p = p;
2865 295049 : d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
2866 295049 : D.aut = Flxq_powers(x, d, T, p);
2867 295049 : x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
2868 295049 : return gerepilecopy(av, x);
2869 : }
2870 :
2871 : GEN
2872 1659 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
2873 : {
2874 1659 : long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
2875 : ulong i;
2876 1659 : pari_sp av = avma;
2877 1659 : GEN xp, V = cgetg(l+2,t_VEC);
2878 1659 : gel(V,1) = polx_Flx(vT); if (l==0) return V;
2879 1659 : gel(V,2) = gcopy(x); if (l==1) return V;
2880 1659 : T = Flx_get_red(T, p);
2881 1659 : d = brent_kung_optpow(dT-1, l-1, 1);
2882 1659 : xp = Flxq_powers(x, d, T, p);
2883 6958 : for(i = 3; i < l+2; i++)
2884 5299 : gel(V,i) = Flx_FlxqV_eval(gel(V,i-1), xp, T, p);
2885 1659 : return gerepilecopy(av, V);
2886 : }
2887 :
2888 : static GEN
2889 477829 : Flxq_autsum_mul(void *E, GEN x, GEN y)
2890 : {
2891 477829 : struct _Flxq *D = (struct _Flxq*)E;
2892 477829 : GEN T = D->T;
2893 477829 : ulong p = D->p;
2894 477829 : GEN phi1 = gel(x,1), a1 = gel(x,2);
2895 477829 : GEN phi2 = gel(y,1), a2 = gel(y,2);
2896 477829 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
2897 477829 : GEN V2 = Flxq_powers(phi2, d, T, p);
2898 477829 : GEN phi3 = Flx_FlxqV_eval(phi1, V2, T, p);
2899 477829 : GEN aphi = Flx_FlxqV_eval(a1, V2, T, p);
2900 477829 : GEN a3 = Flxq_mul(aphi, a2, T, p);
2901 477829 : return mkvec2(phi3, a3);
2902 : }
2903 : static GEN
2904 280218 : Flxq_autsum_sqr(void *E, GEN x)
2905 280218 : { return Flxq_autsum_mul(E, x, x); }
2906 :
2907 : GEN
2908 224118 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
2909 : {
2910 224118 : pari_sp av = avma;
2911 : struct _Flxq D;
2912 224118 : D.T = Flx_get_red(T, p); D.p = p;
2913 224118 : x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
2914 224118 : return gerepilecopy(av, x);
2915 : }
2916 :
2917 : static GEN
2918 310634 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
2919 : {
2920 310634 : struct _Flxq *D = (struct _Flxq*)E;
2921 310634 : GEN T = D->T;
2922 310634 : ulong p = D->p;
2923 310634 : GEN phi1 = gel(x,1), a1 = gel(x,2);
2924 310634 : GEN phi2 = gel(y,1), a2 = gel(y,2);
2925 310634 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
2926 310659 : GEN V1 = Flxq_powers(phi1, d, T, p);
2927 310598 : GEN phi3 = Flx_FlxqV_eval(phi2, V1, T, p);
2928 310594 : GEN aphi = Flx_FlxqV_eval(a2, V1, T, p);
2929 310616 : GEN a3 = Flx_add(a1, aphi, p);
2930 310604 : return mkvec2(phi3, a3);
2931 : }
2932 :
2933 : static GEN
2934 247315 : Flxq_auttrace_sqr(void *E, GEN x)
2935 247315 : { return Flxq_auttrace_mul(E, x, x); }
2936 :
2937 : GEN
2938 311003 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
2939 : {
2940 311003 : pari_sp av = avma;
2941 : struct _Flxq D;
2942 311003 : D.T = Flx_get_red(T, p); D.p = p;
2943 311006 : x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
2944 310993 : return gerepilecopy(av, x);
2945 : }
2946 :
2947 : static long
2948 496468 : bounded_order(ulong p, GEN b, long k)
2949 : {
2950 : long i;
2951 496468 : GEN a=modii(utoi(p),b);
2952 1373189 : for(i=1;i<k;i++)
2953 : {
2954 1168076 : if (equali1(a))
2955 291355 : return i;
2956 876721 : a = modii(muliu(a,p),b);
2957 : }
2958 205113 : return 0;
2959 : }
2960 :
2961 : /*
2962 : n = (p^d-a)\b
2963 : b = bb*p^vb
2964 : p^k = 1 [bb]
2965 : d = m*k+r+vb
2966 : u = (p^k-1)/bb;
2967 : v = (p^(r+vb)-a)/b;
2968 : w = (p^(m*k)-1)/(p^k-1)
2969 : n = p^r*w*u+v
2970 : w*u = p^vb*(p^(m*k)-1)/b
2971 : n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b
2972 : */
2973 :
2974 : static GEN
2975 15805365 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p)
2976 : {
2977 15805365 : pari_sp av=avma;
2978 15805365 : long d = get_Flx_degree(T);
2979 15805365 : GEN an = absi_shallow(n), z, q;
2980 15805365 : if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow(x, n, T, p);
2981 496629 : q = powuu(p, d);
2982 496629 : if (dvdii(q, n))
2983 : {
2984 139 : long vn = logint(an,utoi(p));
2985 139 : GEN autvn = vn==1 ? aut: Flxq_autpow(aut,vn,T,p);
2986 139 : z = Flx_Flxq_eval(x,autvn,T,p);
2987 : } else
2988 : {
2989 496490 : GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
2990 : GEN bb, u, v, autk;
2991 496490 : long vb = Z_lvalrem(b,p,&bb);
2992 496490 : long m, r, k = is_pm1(bb) ? 1 : bounded_order(p,bb,d);
2993 496490 : if (!k || d-vb<k) return Flxq_pow(x,n, T, p);
2994 291370 : m = (d-vb)/k; r = (d-vb)%k;
2995 291370 : u = diviiexact(subiu(powuu(p,k),1),bb);
2996 291370 : v = diviiexact(subii(powuu(p,r+vb),a),b);
2997 291370 : autk = k==1 ? aut: Flxq_autpow(aut,k,T,p);
2998 291370 : if (r)
2999 : {
3000 92970 : GEN autr = r==1 ? aut: Flxq_autpow(aut,r,T,p);
3001 92970 : z = Flx_Flxq_eval(x,autr,T,p);
3002 198400 : } else z = x;
3003 291370 : if (m > 1) z = gel(Flxq_autsum(mkvec2(autk, z), m, T, p), 2);
3004 291370 : if (!is_pm1(u)) z = Flxq_pow(z, u, T, p);
3005 291370 : if (signe(v)) z = Flxq_mul(z, Flxq_pow(x, v, T, p), T, p);
3006 : }
3007 291509 : return gerepileupto(av,signe(n)>0 ? z : Flxq_inv(z,T,p));
3008 : }
3009 :
3010 : static GEN
3011 15798454 : _Flxq_pow(void *data, GEN x, GEN n)
3012 : {
3013 15798454 : struct _Flxq *D = (struct _Flxq*)data;
3014 15798454 : return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p);
3015 : }
3016 :
3017 : static GEN
3018 319130 : _Flxq_rand(void *data)
3019 : {
3020 319130 : pari_sp av=avma;
3021 319130 : struct _Flxq *D = (struct _Flxq*)data;
3022 : GEN z;
3023 : do
3024 : {
3025 322182 : set_avma(av);
3026 322182 : z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3027 322182 : } while (lgpol(z)==0);
3028 319130 : return z;
3029 : }
3030 :
3031 : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3032 : static GEN
3033 14262 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3034 : {
3035 14262 : pari_sp av = avma;
3036 : GEN q,n_q,ord,ordp, op;
3037 :
3038 14262 : if (a == 1UL) return gen_0;
3039 : /* p > 2 */
3040 :
3041 14262 : ordp = utoi(p - 1);
3042 14262 : ord = get_arith_Z(o);
3043 14262 : if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3044 14262 : if (a == p - 1) /* -1 */
3045 1551 : return gerepileuptoint(av, shifti(ord,-1));
3046 12711 : ordp = gcdii(ordp, ord);
3047 12711 : op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3048 :
3049 12711 : q = NULL;
3050 12711 : if (T)
3051 : { /* we want < g > = Fp^* */
3052 12711 : if (!equalii(ord,ordp)) {
3053 4455 : q = diviiexact(ord,ordp);
3054 4455 : g = Flxq_pow(g,q,T,p);
3055 : }
3056 : }
3057 12711 : n_q = Fp_log(utoi(a), utoi(uel(g,2)), op, utoi(p));
3058 12711 : if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
3059 12711 : if (q) n_q = mulii(q, n_q);
3060 12711 : return gerepileuptoint(av, n_q);
3061 : }
3062 :
3063 : static GEN
3064 246976 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3065 : {
3066 246976 : struct _Flxq *f = (struct _Flxq *)E;
3067 246976 : GEN T = f->T;
3068 246976 : ulong p = f->p;
3069 246976 : long d = get_Flx_degree(T);
3070 246976 : if (Flx_equal1(a)) return gen_0;
3071 225583 : if (Flx_equal(a,g)) return gen_1;
3072 36747 : if (!degpol(a))
3073 14262 : return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3074 22485 : if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3075 22457 : return NULL;
3076 28 : return Flxq_log_index(a, g, ord, T, p);
3077 : }
3078 :
3079 : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3080 :
3081 : const struct bb_group *
3082 210295 : get_Flxq_star(void **E, GEN T, ulong p)
3083 : {
3084 210295 : struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3085 210295 : e->T = T; e->p = p; e->aut = Flx_Frobenius(T, p);
3086 210296 : *E = (void*)e; return &Flxq_star;
3087 : }
3088 :
3089 : GEN
3090 12755 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3091 : {
3092 : void *E;
3093 12755 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3094 12755 : return gen_order(a,ord,E,S);
3095 : }
3096 :
3097 : GEN
3098 30300 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3099 : {
3100 : void *E;
3101 30300 : pari_sp av = avma;
3102 30300 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3103 30300 : GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3104 30300 : if (Flxq_log_use_index(gel(F,lg(F)-1), T, p))
3105 5635 : v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3106 30300 : return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
3107 : }
3108 :
3109 : GEN
3110 170362 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3111 : {
3112 170362 : if (!lgpol(a))
3113 : {
3114 3122 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3115 3115 : if (zeta)
3116 0 : *zeta=pol1_Flx(get_Flx_var(T));
3117 3115 : return pol0_Flx(get_Flx_var(T));
3118 : }
3119 : else
3120 : {
3121 : void *E;
3122 167240 : pari_sp av = avma;
3123 167240 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3124 167241 : GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3125 167241 : GEN s = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
3126 167241 : if (s) gerepileall(av, zeta?2:1, &s, zeta);
3127 167241 : return s;
3128 : }
3129 : }
3130 :
3131 : GEN
3132 162077 : Flxq_sqrt(GEN a, GEN T, ulong p)
3133 : {
3134 162077 : return Flxq_sqrtn(a, gen_2, T, p, NULL);
3135 : }
3136 :
3137 : /* assume T irreducible mod p */
3138 : int
3139 358704 : Flxq_issquare(GEN x, GEN T, ulong p)
3140 : {
3141 358704 : if (lgpol(x) == 0 || p == 2) return 1;
3142 355533 : return krouu(Flxq_norm(x,T,p), p) == 1;
3143 : }
3144 :
3145 : /* assume T irreducible mod p */
3146 : int
3147 280 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3148 : {
3149 : pari_sp av;
3150 : GEN m;
3151 280 : if (n==1) return Flxq_issquare(x, T, p);
3152 280 : if (lgpol(x) == 0 || p == 2) return 1;
3153 280 : av = avma;
3154 280 : m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3155 280 : return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3156 : }
3157 :
3158 : GEN
3159 113505 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
3160 : {
3161 113505 : pari_sp av=avma;
3162 113505 : GEN A = Flx_splitting(a,p);
3163 113505 : return gerepileuptoleaf(av, FlxqV_dotproduct(A,sqx,T,p));
3164 : }
3165 :
3166 : GEN
3167 25032 : Flxq_lroot(GEN a, GEN T, long p)
3168 : {
3169 25032 : pari_sp av=avma;
3170 25032 : long n = get_Flx_degree(T), d = degpol(a);
3171 : GEN sqx, V;
3172 25032 : if (n==1) return leafcopy(a);
3173 25032 : if (n==2) return Flxq_powu(a, p, T, p);
3174 25032 : sqx = Flxq_autpow(Flx_Frobenius(T, p), n-1, T, p);
3175 25032 : if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
3176 0 : if (d>=p)
3177 : {
3178 0 : V = Flxq_powers(sqx,p-1,T,p);
3179 0 : return gerepileuptoleaf(av, Flxq_lroot_fast(a,V,T,p));
3180 : } else
3181 0 : return gerepileuptoleaf(av, Flx_Flxq_eval(a,sqx,T,p));
3182 : }
3183 :
3184 : ulong
3185 387724 : Flxq_norm(GEN x, GEN TB, ulong p)
3186 : {
3187 387724 : GEN T = get_Flx_mod(TB);
3188 387724 : ulong y = Flx_resultant(T, x, p);
3189 387724 : ulong L = Flx_lead(T);
3190 387724 : if ( L==1 || lgpol(x)==0) return y;
3191 0 : return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3192 : }
3193 :
3194 : ulong
3195 3212 : Flxq_trace(GEN x, GEN TB, ulong p)
3196 : {
3197 3212 : pari_sp av = avma;
3198 : ulong t;
3199 3212 : GEN T = get_Flx_mod(TB);
3200 3212 : long n = degpol(T)-1;
3201 3212 : GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3202 3212 : t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3203 3212 : return gc_ulong(av, t);
3204 : }
3205 :
3206 : /*x must be reduced*/
3207 : GEN
3208 2771 : Flxq_charpoly(GEN x, GEN TB, ulong p)
3209 : {
3210 2771 : pari_sp ltop=avma;
3211 2771 : GEN T = get_Flx_mod(TB);
3212 2771 : long vs = evalvarn(fetch_var());
3213 2771 : GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3214 2771 : GEN r = Flx_FlxY_resultant(T, xm1, p);
3215 2771 : r[1] = x[1];
3216 2771 : (void)delete_var(); return gerepileupto(ltop, r);
3217 : }
3218 :
3219 : /* Computing minimal polynomial : */
3220 : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3221 : /* in Algebraic Extensions of Finite Fields' */
3222 :
3223 : /* Let v a linear form, return the linear form z->v(tau*z)
3224 : that is, v*(M_tau) */
3225 :
3226 : static GEN
3227 548359 : Flxq_transmul_init(GEN tau, GEN T, ulong p)
3228 : {
3229 : GEN bht;
3230 548359 : GEN h, Tp = get_Flx_red(T, &h);
3231 548360 : long n = degpol(Tp), vT = Tp[1];
3232 548345 : GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3233 548345 : GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3234 548343 : ft[1] = vT; bt[1] = vT;
3235 548343 : if (h)
3236 2850 : bht = Flxn_mul(bt, h, n-1, p);
3237 : else
3238 : {
3239 545493 : GEN bh = Flx_div(Flx_shift(tau, n-1), T, p);
3240 545488 : bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3241 545489 : bht[1] = vT;
3242 : }
3243 548339 : return mkvec3(bt, bht, ft);
3244 : }
3245 :
3246 : static GEN
3247 1448913 : Flxq_transmul(GEN tau, GEN a, long n, ulong p)
3248 : {
3249 1448913 : pari_sp ltop = avma;
3250 : GEN t1, t2, t3, vec;
3251 1448913 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3252 1448913 : if (lgpol(a)==0) return pol0_Flx(a[1]);
3253 1439567 : t2 = Flx_shift(Flx_mul(bt, a, p),1-n);
3254 1439276 : if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
3255 1118950 : t1 = Flx_shift(Flx_mul(ft, a, p),-n);
3256 1118985 : t3 = Flxn_mul(t1, bht, n-1, p);
3257 1119075 : vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3258 1118961 : return gerepileuptoleaf(ltop, vec);
3259 : }
3260 :
3261 : GEN
3262 257360 : Flxq_minpoly(GEN x, GEN T, ulong p)
3263 : {
3264 257360 : pari_sp ltop = avma;
3265 257360 : long vT = get_Flx_var(T), n = get_Flx_degree(T);
3266 : GEN v_x;
3267 257351 : GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3268 257327 : T = Flx_get_red(T, p);
3269 257329 : v_x = Flxq_powers(x, usqrt(2*n), T, p);
3270 531513 : while (lgpol(tau) != 0)
3271 : {
3272 : long i, j, m, k1;
3273 : GEN M, v, tr;
3274 : GEN g_prime, c;
3275 274168 : if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3276 274169 : v = random_Flx(n, vT, p);
3277 274186 : tr = Flxq_transmul_init(tau, T, p);
3278 274171 : v = Flxq_transmul(tr, v, n, p);
3279 274183 : m = 2*(n-degpol(g));
3280 274182 : k1 = usqrt(m);
3281 274183 : tr = Flxq_transmul_init(gel(v_x,k1+1), T, p);
3282 274173 : c = cgetg(m+2,t_VECSMALL);
3283 274205 : c[1] = vT;
3284 1448725 : for (i=0; i<m; i+=k1)
3285 : {
3286 1174544 : long mj = minss(m-i, k1);
3287 5318120 : for (j=0; j<mj; j++)
3288 4143204 : uel(c,m+1-(i+j)) = Flx_dotproduct(v, gel(v_x,j+1), p);
3289 1174916 : v = Flxq_transmul(tr, v, n, p);
3290 : }
3291 274181 : c = Flx_renormalize(c, m+2);
3292 : /* now c contains <v,x^i> , i = 0..m-1 */
3293 274183 : M = Flx_halfgcd(monomial_Flx(1, m, vT), c, p);
3294 274201 : g_prime = gmael(M, 2, 2);
3295 274201 : if (degpol(g_prime) < 1) continue;
3296 270575 : g = Flx_mul(g, g_prime, p);
3297 270558 : tau = Flxq_mul(tau, Flx_FlxqV_eval(g_prime, v_x, T, p), T, p);
3298 : }
3299 257308 : g = Flx_normalize(g,p);
3300 257346 : return gerepileuptoleaf(ltop,g);
3301 : }
3302 :
3303 : GEN
3304 20 : Flxq_conjvec(GEN x, GEN T, ulong p)
3305 : {
3306 20 : long i, l = 1+get_Flx_degree(T);
3307 20 : GEN z = cgetg(l,t_COL);
3308 20 : T = Flx_get_red(T,p);
3309 20 : gel(z,1) = Flx_copy(x);
3310 88 : for (i=2; i<l; i++) gel(z,i) = Flxq_powu(gel(z,i-1), p, T, p);
3311 20 : return z;
3312 : }
3313 :
3314 : GEN
3315 7194 : gener_Flxq(GEN T, ulong p, GEN *po)
3316 : {
3317 7194 : long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
3318 : ulong p_1;
3319 : GEN g, L, L2, o, q, F;
3320 : pari_sp av0, av;
3321 :
3322 7194 : if (f == 1) {
3323 : GEN fa;
3324 28 : o = utoipos(p-1);
3325 28 : fa = Z_factor(o);
3326 28 : L = gel(fa,1);
3327 28 : L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
3328 28 : g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
3329 28 : if (po) *po = mkvec2(o, fa);
3330 28 : return g;
3331 : }
3332 :
3333 7166 : av0 = avma; p_1 = p - 1;
3334 7166 : q = diviuexact(subiu(powuu(p,f), 1), p_1);
3335 :
3336 7166 : L = cgetg(1, t_VECSMALL);
3337 7166 : if (p > 3)
3338 : {
3339 1797 : ulong t = p_1 >> vals(p_1);
3340 1797 : GEN P = gel(factoru(t), 1);
3341 1797 : L = cgetg_copy(P, &i);
3342 2513 : while (--i) L[i] = p_1 / P[i];
3343 : }
3344 7166 : o = factor_pn_1(utoipos(p),f);
3345 7166 : L2 = leafcopy( gel(o, 1) );
3346 17553 : for (i = j = 1; i < lg(L2); i++)
3347 : {
3348 10387 : if (umodui(p_1, gel(L2,i)) == 0) continue;
3349 6012 : gel(L2,j++) = diviiexact(q, gel(L2,i));
3350 : }
3351 7166 : setlg(L2, j);
3352 7166 : F = Flx_Frobenius(T, p);
3353 17120 : for (av = avma;; set_avma(av))
3354 9954 : {
3355 : GEN tt;
3356 17120 : g = random_Flx(f, vT, p);
3357 17120 : if (degpol(g) < 1) continue;
3358 10914 : if (p == 2) tt = g;
3359 : else
3360 : {
3361 7139 : ulong t = Flxq_norm(g, T, p);
3362 7139 : if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
3363 4160 : tt = Flxq_powu(g, p_1>>1, T, p);
3364 : }
3365 14077 : for (i = 1; i < j; i++)
3366 : {
3367 6911 : GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p);
3368 6911 : if (!degpol(a) && uel(a,2) == p_1) break;
3369 : }
3370 7935 : if (i == j) break;
3371 : }
3372 7166 : if (!po)
3373 : {
3374 180 : set_avma((pari_sp)g);
3375 180 : g = gerepileuptoleaf(av0, g);
3376 : }
3377 : else {
3378 6986 : *po = mkvec2(subiu(powuu(p,f), 1), o);
3379 6986 : gerepileall(av0, 2, &g, po);
3380 : }
3381 7166 : return g;
3382 : }
3383 :
3384 : static GEN
3385 483266 : _Flxq_neg(void *E, GEN x)
3386 483266 : { struct _Flxq *s = (struct _Flxq *)E;
3387 483266 : return Flx_neg(x,s->p); }
3388 :
3389 : static GEN
3390 1483609 : _Flxq_rmul(void *E, GEN x, GEN y)
3391 1483609 : { struct _Flxq *s = (struct _Flxq *)E;
3392 1483609 : return Flx_mul(x,y,s->p); }
3393 :
3394 : static GEN
3395 17199 : _Flxq_inv(void *E, GEN x)
3396 17199 : { struct _Flxq *s = (struct _Flxq *)E;
3397 17199 : return Flxq_inv(x,s->T,s->p); }
3398 :
3399 : static int
3400 147511 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
3401 :
3402 : static GEN
3403 22841 : _Flxq_s(void *E, long x)
3404 22841 : { struct _Flxq *s = (struct _Flxq *)E;
3405 22841 : ulong u = x<0 ? s->p+x: (ulong)x;
3406 22841 : return Fl_to_Flx(u, get_Flx_var(s->T));
3407 : }
3408 :
3409 : static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
3410 : _Flxq_inv,_Flxq_equal0,_Flxq_s};
3411 :
3412 66420 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
3413 : {
3414 66420 : GEN z = new_chunk(sizeof(struct _Flxq));
3415 66420 : struct _Flxq *e = (struct _Flxq *) z;
3416 66420 : e->T = Flx_get_red(T, p); e->p = p; *E = (void*)e;
3417 66420 : return &Flxq_field;
3418 : }
3419 :
3420 : /***********************************************************************/
3421 : /** **/
3422 : /** Flxn **/
3423 : /** **/
3424 : /***********************************************************************/
3425 :
3426 : GEN
3427 3194 : Flx_invLaplace(GEN x, ulong p)
3428 : {
3429 3194 : long i, d = degpol(x);
3430 : ulong t;
3431 : GEN y;
3432 3194 : if (d <= 1) return Flx_copy(x);
3433 3194 : t = Fl_inv(factorial_Fl(d, p), p);
3434 3196 : y = cgetg(d+3, t_VECSMALL);
3435 3196 : y[1] = x[1];
3436 111087 : for (i=d; i>=2; i--)
3437 : {
3438 107891 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3439 107890 : t = Fl_mul(t, i, p);
3440 : }
3441 3196 : uel(y,3) = uel(x,3);
3442 3196 : uel(y,2) = uel(x,2);
3443 3196 : return y;
3444 : }
3445 :
3446 : GEN
3447 1598 : Flx_Laplace(GEN x, ulong p)
3448 : {
3449 1598 : long i, d = degpol(x);
3450 1598 : ulong t = 1;
3451 : GEN y;
3452 1598 : if (d <= 1) return Flx_copy(x);
3453 1598 : y = cgetg(d+3, t_VECSMALL);
3454 1598 : y[1] = x[1];
3455 1598 : uel(y,2) = uel(x,2);
3456 1598 : uel(y,3) = uel(x,3);
3457 55640 : for (i=2; i<=d; i++)
3458 : {
3459 54042 : t = Fl_mul(t, i%p, p);
3460 54042 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3461 : }
3462 1598 : return y;
3463 : }
3464 :
3465 : GEN
3466 1235001 : Flxn_red(GEN a, long n)
3467 : {
3468 1235001 : long i, L, l = lg(a);
3469 : GEN b;
3470 1235001 : if (l == 2 || !n) return zero_Flx(a[1]);
3471 1215347 : L = n+2; if (L > l) L = l;
3472 1215347 : b = cgetg(L, t_VECSMALL); b[1] = a[1];
3473 13550249 : for (i=2; i<L; i++) b[i] = a[i];
3474 1215075 : return Flx_renormalize(b,L);
3475 : }
3476 :
3477 : GEN
3478 1201956 : Flxn_mul(GEN a, GEN b, long n, ulong p)
3479 1201956 : { return Flxn_red(Flx_mul(a, b, p), n); }
3480 :
3481 : GEN
3482 0 : Flxn_sqr(GEN a, long n, ulong p)
3483 0 : { return Flxn_red(Flx_sqr(a, p), n); }
3484 :
3485 : /* (f*g) \/ x^n */
3486 : static GEN
3487 30627 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p)
3488 : {
3489 30627 : return Flx_shift(Flx_mul(f,g, p),-n);
3490 : }
3491 :
3492 : static GEN
3493 21654 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p)
3494 : {
3495 21654 : GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
3496 21657 : return Flx_add(Flx_mulhigh_i(fl, g, n2, p), Flxn_mul(fh, g, n - n2, p), p);
3497 : }
3498 :
3499 : GEN
3500 3384 : Flxn_inv(GEN f, long e, ulong p)
3501 : {
3502 3384 : pari_sp av = avma, av2;
3503 : ulong mask;
3504 : GEN W;
3505 3384 : long n=1;
3506 3384 : if (lg(f)==2) pari_err_INV("Flxn_inv",f);
3507 3384 : W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
3508 3384 : mask = quadratic_prec_mask(e);
3509 3385 : av2 = avma;
3510 18093 : for (;mask>1;)
3511 : {
3512 : GEN u, fr;
3513 14709 : long n2 = n;
3514 14709 : n<<=1; if (mask & 1) n--;
3515 14709 : mask >>= 1;
3516 14709 : fr = Flxn_red(f, n);
3517 14707 : u = Flxn_mul(W, Flxn_mulhigh(fr, W, n2, n, p), n-n2, p);
3518 14708 : W = Flx_sub(W, Flx_shift(u, n2), p);
3519 14709 : if (gc_needed(av2,2))
3520 : {
3521 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_inv, e = %ld", n);
3522 0 : W = gerepileupto(av2, W);
3523 : }
3524 : }
3525 3384 : return gerepileupto(av, W);
3526 : }
3527 :
3528 : GEN
3529 2030 : Flxn_expint(GEN h, long e, ulong p)
3530 : {
3531 2030 : pari_sp av = avma, av2;
3532 2030 : long v = h[1], n=1;
3533 2030 : GEN f = pol1_Flx(v), g = pol1_Flx(v);
3534 2030 : ulong mask = quadratic_prec_mask(e);
3535 2030 : av2 = avma;
3536 8980 : for (;mask>1;)
3537 : {
3538 : GEN u, w;
3539 8980 : long n2 = n;
3540 8980 : n<<=1; if (mask & 1) n--;
3541 8980 : mask >>= 1;
3542 8980 : u = Flxn_mul(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p), n-n2, p);
3543 8981 : u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
3544 8981 : w = Flxn_mul(f, Flx_integXn(u, n2-1, p), n-n2, p);
3545 8981 : f = Flx_add(f, Flx_shift(w, n2), p);
3546 8981 : if (mask<=1) break;
3547 6951 : u = Flxn_mul(g, Flxn_mulhigh(f, g, n2, n, p), n-n2, p);
3548 6950 : g = Flx_sub(g, Flx_shift(u, n2), p);
3549 6950 : if (gc_needed(av2,2))
3550 : {
3551 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
3552 0 : gerepileall(av2, 2, &f, &g);
3553 : }
3554 : }
3555 2030 : return gerepileupto(av, f);
3556 : }
3557 :
3558 : GEN
3559 0 : Flxn_exp(GEN h, long e, ulong p)
3560 : {
3561 0 : if (degpol(h)<1 || uel(h,2)!=0)
3562 0 : pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
3563 0 : return Flxn_expint(Flx_deriv(h, p), e, p);
3564 : }
3565 :
3566 : INLINE GEN
3567 8419 : Flxn_recip(GEN x, long n)
3568 : {
3569 8419 : GEN z=Flx_recipspec(x+2,lgpol(x),n);
3570 8419 : z[1]=x[1];
3571 8419 : return z;
3572 : }
3573 :
3574 : GEN
3575 3195 : Flx_Newton(GEN P, long n, ulong p)
3576 : {
3577 3195 : pari_sp av = avma;
3578 3195 : long d = degpol(P);
3579 3195 : GEN dP = Flxn_recip(Flx_deriv(P, p), d);
3580 3195 : GEN Q = Flxn_mul(Flxn_inv(Flxn_recip(P, d+1), n, p), dP, n, p);
3581 3194 : return gerepileuptoleaf(av, Q);
3582 : }
3583 :
3584 : GEN
3585 2030 : Flx_fromNewton(GEN P, ulong p)
3586 : {
3587 2030 : pari_sp av = avma;
3588 2030 : ulong n = Flx_constant(P)+1;
3589 2030 : GEN z = Flx_neg(Flx_shift(P, -1), p);
3590 2030 : GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
3591 2030 : return gerepileuptoleaf(av, Q);
3592 : }
3593 :
3594 : static long
3595 434 : newtonlogint(ulong n, ulong pp)
3596 : {
3597 434 : long s = 0;
3598 1960 : while (n > pp)
3599 : {
3600 1526 : s += ulogint(n, pp);
3601 1526 : n = (n+1)>>1;
3602 : }
3603 434 : return s;
3604 : }
3605 :
3606 : static void
3607 655 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
3608 : {
3609 : long i;
3610 : ulong e;
3611 655 : GEN P = cgetg(d+1, t_VECSMALL);
3612 655 : GEN V = cgetg(d+1, t_VECSMALL);
3613 27292 : for (i=1, e=1; i<=d; i++, e++)
3614 : {
3615 26637 : if (e==p)
3616 : {
3617 9512 : e = 0;
3618 9512 : V[i] = u_lvalrem(i, p, &uel(P,i));
3619 : } else
3620 : {
3621 17125 : V[i] = 0; uel(P,i) = i;
3622 : }
3623 : }
3624 655 : *pt_P = P; *pt_V = V;
3625 655 : }
3626 :
3627 : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
3628 : * val large enough to compensate for the power of p in the factorials */
3629 :
3630 : static GEN
3631 434 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long var)
3632 : {
3633 434 : pari_sp av = avma;
3634 434 : long i, d = n-1, w;
3635 : GEN y, W, E, t;
3636 434 : init_invlaplace(d, p, &E, &W);
3637 434 : t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
3638 434 : w = zv_sum(W);
3639 434 : if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
3640 434 : y = cgetg(d+3,t_POL);
3641 434 : y[1] = evalsigne(1) | evalvarn(var);
3642 21203 : for (i=d; i>=1; i--)
3643 : {
3644 20769 : gel(y,i+2) = t;
3645 20769 : t = Fp_mulu(t, uel(E,i), q);
3646 20769 : if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
3647 : }
3648 434 : gel(y,2) = t;
3649 434 : return gerepilecopy(av, ZX_renormalize(y, d+3));
3650 : }
3651 :
3652 : static GEN
3653 2032 : Flx_composedsum(GEN P, GEN Q, ulong p)
3654 : {
3655 2032 : long n = 1 + degpol(P)*degpol(Q);
3656 2032 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
3657 2032 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
3658 : GEN R;
3659 2032 : if (p >= (ulong)n)
3660 : {
3661 1598 : GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
3662 1598 : GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
3663 1598 : GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
3664 1598 : R = Flx_fromNewton(L, p);
3665 : } else
3666 : {
3667 434 : long v = factorial_lval(n-1, p);
3668 434 : long w = 1 + newtonlogint(n-1, p);
3669 434 : GEN pv = powuu(p, v);
3670 434 : GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
3671 434 : GEN iL = ZpX_invLaplace_init(n, q, p, v, varn(P));
3672 434 : GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
3673 434 : GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
3674 434 : GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
3675 434 : GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
3676 434 : R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
3677 : }
3678 2032 : return Flx_Fl_mul(R, lead, p);
3679 : }
3680 :
3681 : GEN
3682 2032 : Flx_direct_compositum(GEN a, GEN b, ulong p)
3683 : {
3684 2032 : return Flx_composedsum(a, b, p);
3685 : }
3686 :
3687 : static GEN
3688 628 : _Flx_direct_compositum(void *E, GEN a, GEN b)
3689 628 : { return Flx_direct_compositum(a, b, (ulong)E); }
3690 :
3691 : GEN
3692 4982 : FlxV_direct_compositum(GEN V, ulong p)
3693 : {
3694 4982 : return gen_product(V, (void *)p, &_Flx_direct_compositum);
3695 : }
3696 :
3697 : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
3698 : static GEN
3699 0 : Fl_Xp1_powu(ulong n, ulong p, long v)
3700 : {
3701 0 : ulong k, d = (n + 1) >> 1;
3702 0 : GEN C, V = identity_zv(d);
3703 :
3704 0 : Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
3705 0 : C = cgetg(n+3, t_VECSMALL);
3706 0 : C[1] = v;
3707 0 : uel(C,2) = 1UL;
3708 0 : uel(C,3) = n%p;
3709 0 : uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
3710 : /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
3711 0 : if (SMALL_ULONG(p))
3712 0 : for (k = 3; k <= d; k++)
3713 0 : uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
3714 : else
3715 : {
3716 0 : ulong pi = get_Fl_red(p);
3717 0 : for (k = 3; k <= d; k++)
3718 0 : uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
3719 : }
3720 0 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
3721 0 : return C; /* normalized */
3722 : }
3723 :
3724 : /* p arbitrary */
3725 : GEN
3726 13279 : Flx_translate1_basecase(GEN P, ulong p)
3727 : {
3728 13279 : GEN R = Flx_copy(P);
3729 13277 : long i, k, n = degpol(P);
3730 72853 : for (i = 1; i <= n; i++)
3731 498556 : for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
3732 13279 : return R;
3733 : }
3734 :
3735 : static int
3736 13500 : translate_basecase(long n, ulong p)
3737 : {
3738 : #ifdef LONG_IS_64BIT
3739 11588 : if (p <= 19) return n < 40;
3740 8030 : if (p < 1UL<<30) return n < 58;
3741 0 : if (p < 1UL<<59) return n < 100;
3742 0 : if (p < 1UL<<62) return n < 120;
3743 0 : if (p < 1UL<<63) return n < 240;
3744 0 : return n < 250;
3745 : #else
3746 1912 : if (p <= 13) return n < 18;
3747 1428 : if (p <= 17) return n < 22;
3748 1376 : if (p <= 29) return n < 39;
3749 1202 : if (p <= 67) return n < 69;
3750 983 : if (p < 1UL<< 15) return n < 80;
3751 0 : if (p < 1UL<< 16) return n < 100;
3752 0 : if (p < 1UL<< 28) return n < 300;
3753 0 : return n < 650;
3754 : #endif
3755 : }
3756 : /* assume p prime */
3757 : GEN
3758 11848 : Flx_translate1(GEN P, ulong p)
3759 : {
3760 11848 : long d, n = degpol(P);
3761 : GEN R, Q, S;
3762 11847 : if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
3763 : /* n > 0 */
3764 0 : d = n >> 1;
3765 0 : if ((ulong)n < p)
3766 : {
3767 0 : R = Flx_translate1(Flxn_red(P, d), p);
3768 0 : Q = Flx_translate1(Flx_shift(P, -d), p);
3769 0 : S = Fl_Xp1_powu(d, p, P[1]);
3770 0 : return Flx_add(Flx_mul(Q, S, p), R, p);
3771 : }
3772 : else
3773 : {
3774 : ulong q;
3775 0 : if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
3776 0 : R = Flx_translate1(Flxn_red(P, q), p);
3777 0 : Q = Flx_translate1(Flx_shift(P, -q), p);
3778 0 : S = Flx_add(Flx_shift(Q, q), Q, p);
3779 0 : return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
3780 : }
3781 : }
3782 :
3783 : static GEN
3784 221 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
3785 : {
3786 221 : ulong k, d = n >> 1, c, v = 0;
3787 221 : GEN C, V, W, U = upowers(p, e-1);
3788 221 : init_invlaplace(d, p, &V, &W);
3789 221 : Flv_inv_inplace(V, q);
3790 221 : C = cgetg(n+3, t_VECSMALL);
3791 221 : C[1] = vs;
3792 221 : uel(C,2) = 1UL;
3793 221 : uel(C,3) = n%q;
3794 221 : v = u_lvalrem(n, p, &c);
3795 5868 : for (k = 2; k <= d; k++)
3796 : {
3797 : ulong w;
3798 5647 : v += u_lvalrem(n-k+1, p, &w) - W[k];
3799 5647 : c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
3800 5647 : uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
3801 : }
3802 6160 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
3803 221 : return C; /* normalized */
3804 : }
3805 :
3806 : GEN
3807 1653 : zlx_translate1(GEN P, ulong p, long e)
3808 : {
3809 1653 : ulong d, q = upowuu(p,e), n = degpol(P);
3810 : GEN R, Q, S;
3811 1653 : if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
3812 : /* n > 0 */
3813 221 : d = n >> 1;
3814 221 : R = zlx_translate1(Flxn_red(P, d), p, e);
3815 221 : Q = zlx_translate1(Flx_shift(P, -d), p, e);
3816 221 : S = zl_Xp1_powu(d, p, q, e, P[1]);
3817 221 : return Flx_add(Flx_mul(Q, S, q), R, q);
3818 : }
3819 :
3820 : /***********************************************************************/
3821 : /** **/
3822 : /** Fl2 **/
3823 : /** **/
3824 : /***********************************************************************/
3825 : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
3826 : a non-square D.
3827 : */
3828 :
3829 : INLINE GEN
3830 6249715 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
3831 :
3832 : GEN
3833 1679315 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
3834 : {
3835 : ulong xaya, xbyb, Db2, mid;
3836 : ulong z1, z2;
3837 1679315 : ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
3838 1679315 : xaya = Fl_mul_pre(x1,y1,p,pi);
3839 1679389 : if (x2==0 && y2==0) return mkF2(xaya,0);
3840 1625172 : if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
3841 1604347 : if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
3842 1604091 : xbyb = Fl_mul_pre(x2,y2,p,pi);
3843 1604088 : mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
3844 1604090 : Db2 = Fl_mul_pre(D, xbyb, p,pi);
3845 1604097 : z1 = Fl_add(xaya,Db2,p);
3846 1604094 : z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
3847 1604084 : return mkF2(z1,z2);
3848 : }
3849 :
3850 : GEN
3851 4232502 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
3852 : {
3853 4232502 : ulong a = x[1], b = x[2];
3854 : ulong a2, Db2, ab;
3855 4232502 : a2 = Fl_sqr_pre(a,p,pi);
3856 4232991 : if (b==0) return mkF2(a2,0);
3857 4058419 : Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
3858 4058350 : ab = Fl_mul_pre(a,b,p,pi);
3859 4058383 : return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
3860 : }
3861 :
3862 : ulong
3863 66104 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
3864 : {
3865 66104 : ulong a2 = Fl_sqr_pre(x[1],p,pi);
3866 66104 : return x[2]? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(x[2], p,pi), p,pi), p): a2;
3867 : }
3868 :
3869 : GEN
3870 167357 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
3871 : {
3872 : ulong n, ni;
3873 167357 : if (x[2] == 0) return mkF2(Fl_inv(x[1],p),0);
3874 143265 : n = Fl_sub(Fl_sqr_pre(x[1], p,pi),
3875 143265 : Fl_mul_pre(D, Fl_sqr_pre(x[2], p,pi), p,pi), p);
3876 143264 : ni = Fl_inv(n,p);
3877 143265 : return mkF2(Fl_mul_pre(x[1], ni, p,pi),
3878 143265 : Fl_neg(Fl_mul_pre(x[2], ni, p,pi), p));
3879 : }
3880 :
3881 : int
3882 379769 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
3883 :
3884 : struct _Fl2 {
3885 : ulong p, pi, D;
3886 : };
3887 :
3888 : static GEN
3889 4232391 : _Fl2_sqr(void *data, GEN x)
3890 : {
3891 4232391 : struct _Fl2 *D = (struct _Fl2*)data;
3892 4232391 : return Fl2_sqr_pre(x, D->D, D->p, D->pi);
3893 : }
3894 : static GEN
3895 1651578 : _Fl2_mul(void *data, GEN x, GEN y)
3896 : {
3897 1651578 : struct _Fl2 *D = (struct _Fl2*)data;
3898 1651578 : return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
3899 : }
3900 :
3901 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
3902 : GEN
3903 568011 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
3904 : {
3905 568011 : pari_sp av = avma;
3906 : struct _Fl2 d;
3907 : GEN y;
3908 568011 : long s = signe(n);
3909 568011 : if (!s) return mkF2(1,0);
3910 503556 : if (s < 0)
3911 167357 : x = Fl2_inv_pre(x,D,p,pi);
3912 503555 : if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
3913 370309 : d.p = p; d.pi = pi; d.D=D;
3914 370309 : y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
3915 370305 : return gerepileuptoleaf(av, y);
3916 : }
3917 :
3918 : static GEN
3919 568010 : _Fl2_pow(void *data, GEN x, GEN n)
3920 : {
3921 568010 : struct _Fl2 *D = (struct _Fl2*)data;
3922 568010 : return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
3923 : }
3924 :
3925 : static GEN
3926 96040 : _Fl2_rand(void *data)
3927 : {
3928 96040 : struct _Fl2 *D = (struct _Fl2*)data;
3929 96040 : ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
3930 96039 : return mkF2(a,b);
3931 : }
3932 :
3933 : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
3934 : hash_GEN, zv_equal, Fl2_equal1, NULL};
3935 :
3936 : GEN
3937 64455 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
3938 : {
3939 : struct _Fl2 E;
3940 : GEN o;
3941 64455 : if (a[1]==0 && a[2]==0)
3942 : {
3943 0 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3944 0 : if (zeta) *zeta=mkF2(1,0);
3945 0 : return zv_copy(a);
3946 : }
3947 64455 : E.p=p; E.pi = pi; E.D = D;
3948 64455 : o = subiu(powuu(p,2), 1);
3949 64455 : return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
3950 : }
3951 :
3952 : GEN
3953 10108 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
3954 : {
3955 : GEN p1;
3956 10108 : long i = lg(x)-1;
3957 10108 : if (i <= 2)
3958 1883 : return mkF2(i == 2? x[2]: 0, 0);
3959 8225 : p1 = mkF2(x[i], 0);
3960 35952 : for (i--; i>=2; i--)
3961 : {
3962 27727 : p1 = Fl2_mul_pre(p1, y, D, p, pi);
3963 27727 : uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
3964 : }
3965 8225 : return p1;
3966 : }
3967 :
3968 : /***********************************************************************/
3969 : /** **/
3970 : /** FlxV **/
3971 : /** **/
3972 : /***********************************************************************/
3973 : /* FlxV are t_VEC with Flx coefficients. */
3974 :
3975 : GEN
3976 0 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
3977 : {
3978 0 : pari_sp ltop=avma;
3979 : long i;
3980 0 : GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
3981 0 : for(i=2;i<lg(V);i++)
3982 0 : z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
3983 0 : return gerepileuptoleaf(ltop,z);
3984 : }
3985 :
3986 : GEN
3987 0 : ZXV_to_FlxV(GEN x, ulong p)
3988 0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
3989 :
3990 : GEN
3991 1486558 : ZXT_to_FlxT(GEN x, ulong p)
3992 : {
3993 1486558 : if (typ(x) == t_POL)
3994 1435847 : return ZX_to_Flx(x, p);
3995 : else
3996 167744 : pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
3997 : }
3998 :
3999 : GEN
4000 34905 : FlxV_to_Flm(GEN x, long n)
4001 170420 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
4002 :
4003 : GEN
4004 0 : FlxV_red(GEN x, ulong p)
4005 0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4006 :
4007 : GEN
4008 208566 : FlxT_red(GEN x, ulong p)
4009 : {
4010 208566 : if (typ(x) == t_VECSMALL)
4011 140967 : return Flx_red(x, p);
4012 : else
4013 227559 : pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4014 : }
4015 :
4016 : GEN
4017 113505 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4018 : {
4019 113505 : long i, lx = lg(x);
4020 : pari_sp av;
4021 : GEN c;
4022 113505 : if (lx == 1) return pol0_Flx(get_Flx_var(T));
4023 113505 : av = avma; c = Flx_mul(gel(x,1),gel(y,1), p);
4024 463799 : for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul(gel(x,i),gel(y,i), p), p);
4025 113505 : return gerepileuptoleaf(av, Flx_rem(c,T,p));
4026 : }
4027 :
4028 : GEN
4029 1620 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4030 : {
4031 1620 : long i, l = minss(lg(x), lg(y));
4032 : pari_sp av;
4033 : GEN c;
4034 1620 : if (l == 2) return pol0_Flx(get_Flx_var(T));
4035 1620 : av = avma; c = Flx_mul(gel(x,2),gel(y,2), p);
4036 4686 : for (i=3; i<l; i++) c = Flx_add(c, Flx_mul(gel(x,i),gel(y,i), p), p);
4037 1620 : return gerepileuptoleaf(av, Flx_rem(c,T,p));
4038 : }
4039 :
4040 : GEN
4041 218135 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4042 : {
4043 218135 : long i, l = lg(z);
4044 218135 : GEN y = cgetg(l, t_VECSMALL);
4045 7157621 : for (i=1; i<l; i++)
4046 6939482 : uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4047 218139 : return y;
4048 : }
4049 :
4050 : /***********************************************************************/
4051 : /** **/
4052 : /** FlxM **/
4053 : /** **/
4054 : /***********************************************************************/
4055 :
4056 : GEN
4057 18130 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4058 : {
4059 18130 : long i, l = lg(z);
4060 18130 : GEN y = cgetg(l, t_MAT);
4061 236262 : for (i=1; i<l; i++)
4062 218134 : gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4063 18128 : return y;
4064 : }
4065 :
4066 : GEN
4067 0 : zero_FlxC(long n, long sv)
4068 : {
4069 : long i;
4070 0 : GEN x = cgetg(n + 1, t_COL);
4071 0 : GEN z = zero_Flx(sv);
4072 0 : for (i = 1; i <= n; i++)
4073 0 : gel(x, i) = z;
4074 0 : return x;
4075 : }
4076 :
4077 : GEN
4078 0 : FlxC_neg(GEN x, ulong p)
4079 0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4080 :
4081 : GEN
4082 0 : FlxC_sub(GEN x, GEN y, ulong p)
4083 0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4084 :
4085 : GEN
4086 0 : zero_FlxM(long r, long c, long sv)
4087 : {
4088 : long j;
4089 0 : GEN x = cgetg(c + 1, t_MAT);
4090 0 : GEN z = zero_FlxC(r, sv);
4091 0 : for (j = 1; j <= c; j++)
4092 0 : gel(x, j) = z;
4093 0 : return x;
4094 : }
4095 :
4096 : GEN
4097 0 : FlxM_neg(GEN x, ulong p)
4098 0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4099 :
4100 : GEN
4101 0 : FlxM_sub(GEN x, GEN y, ulong p)
4102 0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4103 :
4104 : GEN
4105 0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4106 0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4107 :
4108 : GEN
4109 0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4110 0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4111 :
4112 : static GEN
4113 47494 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4114 : long i, j, l, lc;
4115 47494 : GEN N = cgetg_copy(M, &l), x;
4116 47494 : if (l == 1)
4117 0 : return N;
4118 47494 : lc = lgcols(M);
4119 221155 : for (j = 1; j < l; j++) {
4120 173661 : gel(N, j) = cgetg(lc, t_COL);
4121 862754 : for (i = 1; i < lc; i++) {
4122 689093 : x = gcoeff(M, i, j);
4123 689093 : gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4124 : }
4125 : }
4126 47494 : return N;
4127 : }
4128 :
4129 : static GEN
4130 639654 : kron_pack_Flx_spec_half(GEN x, long l) {
4131 639654 : if (l == 0)
4132 252384 : return gen_0;
4133 387270 : return Flx_to_int_halfspec(x, l);
4134 : }
4135 :
4136 : static GEN
4137 46050 : kron_pack_Flx_spec(GEN x, long l) {
4138 : long i;
4139 : GEN w, y;
4140 46050 : if (l == 0)
4141 8595 : return gen_0;
4142 37455 : y = cgetipos(l + 2);
4143 143761 : for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4144 106306 : *w = x[i];
4145 37455 : return y;
4146 : }
4147 :
4148 : static GEN
4149 3389 : kron_pack_Flx_spec_2(GEN x, long l) {
4150 3389 : return Flx_eval2BILspec(x, 2, l);
4151 : }
4152 :
4153 : static GEN
4154 0 : kron_pack_Flx_spec_3(GEN x, long l) {
4155 0 : return Flx_eval2BILspec(x, 3, l);
4156 : }
4157 :
4158 : static GEN
4159 36194 : kron_unpack_Flx(GEN z, ulong p)
4160 : {
4161 36194 : long i, l = lgefint(z);
4162 36194 : GEN x = cgetg(l, t_VECSMALL), w;
4163 184151 : for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4164 147957 : x[i] = ((ulong) *w) % p;
4165 36194 : return Flx_renormalize(x, l);
4166 : }
4167 :
4168 : static GEN
4169 2930 : kron_unpack_Flx_2(GEN x, ulong p) {
4170 2930 : long d = (lgefint(x)-1)/2 - 1;
4171 2930 : return Z_mod2BIL_Flx_2(x, d, p);
4172 : }
4173 :
4174 : static GEN
4175 0 : kron_unpack_Flx_3(GEN x, ulong p) {
4176 0 : long d = lgefint(x)/3 - 1;
4177 0 : return Z_mod2BIL_Flx_3(x, d, p);
4178 : }
4179 :
4180 : static GEN
4181 105745 : FlxM_pack_ZM_bits(GEN M, long b)
4182 : {
4183 : long i, j, l, lc;
4184 105745 : GEN N = cgetg_copy(M, &l), x;
4185 105745 : if (l == 1)
4186 0 : return N;
4187 105745 : lc = lgcols(M);
4188 425524 : for (j = 1; j < l; j++) {
4189 319779 : gel(N, j) = cgetg(lc, t_COL);
4190 5320624 : for (i = 1; i < lc; i++) {
4191 5000845 : x = gcoeff(M, i, j);
4192 5000845 : gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
4193 : }
4194 : }
4195 105745 : return N;
4196 : }
4197 :
4198 : static GEN
4199 23747 : ZM_unpack_FlxqM(GEN M, GEN T, ulong p, GEN (*unpack)(GEN, ulong))
4200 : {
4201 23747 : long i, j, l, lc, sv = get_Flx_var(T);
4202 23747 : GEN N = cgetg_copy(M, &l), x;
4203 23747 : if (l == 1)
4204 0 : return N;
4205 23747 : lc = lgcols(M);
4206 128456 : for (j = 1; j < l; j++) {
4207 104709 : gel(N, j) = cgetg(lc, t_COL);
4208 575409 : for (i = 1; i < lc; i++) {
4209 470700 : x = unpack(gcoeff(M, i, j), p);
4210 470700 : x[1] = sv;
4211 470700 : gcoeff(N, i, j) = Flx_rem(x, T, p);
4212 : }
4213 : }
4214 23747 : return N;
4215 : }
4216 :
4217 : static GEN
4218 52913 : ZM_unpack_FlxqM_bits(GEN M, long b, GEN T, ulong p)
4219 : {
4220 52913 : long i, j, l, lc, sv = get_Flx_var(T);
4221 52913 : GEN N = cgetg_copy(M, &l), x;
4222 52913 : if (l == 1)
4223 0 : return N;
4224 52913 : lc = lgcols(M);
4225 52913 : if (b < BITS_IN_LONG) {
4226 170070 : for (j = 1; j < l; j++) {
4227 118255 : gel(N, j) = cgetg(lc, t_COL);
4228 2720984 : for (i = 1; i < lc; i++) {
4229 2602729 : x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
4230 2602729 : x[1] = sv;
4231 2602729 : gcoeff(N, i, j) = Flx_rem(x, T, p);
4232 : }
4233 : }
4234 : } else {
4235 1098 : ulong pi = get_Fl_red(p);
4236 7410 : for (j = 1; j < l; j++) {
4237 6312 : gel(N, j) = cgetg(lc, t_COL);
4238 135663 : for (i = 1; i < lc; i++) {
4239 129351 : x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
4240 129351 : x[1] = sv;
4241 129351 : gcoeff(N, i, j) = Flx_rem(x, T, p);
4242 : }
4243 : }
4244 : }
4245 52913 : return N;
4246 : }
4247 :
4248 : GEN
4249 76660 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
4250 : {
4251 76660 : pari_sp av = avma;
4252 76660 : long b, d = get_Flx_degree(T), n = lg(A) - 1;
4253 : GEN C, D, z;
4254 : GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
4255 76660 : int is_sqr = A==B;
4256 :
4257 76660 : z = muliu(muliu(sqru(p - 1), d), n);
4258 76660 : b = expi(z) + 1;
4259 : /* only do expensive bit-packing if it saves at least 1 limb */
4260 76660 : if (b <= BITS_IN_HALFULONG) {
4261 73686 : if (nbits2nlong(d*b) == (d + 1)/2)
4262 22897 : b = BITS_IN_HALFULONG;
4263 : } else {
4264 2974 : long l = lgefint(z) - 2;
4265 2974 : if (nbits2nlong(d*b) == d*l)
4266 850 : b = l*BITS_IN_LONG;
4267 : }
4268 76660 : set_avma(av);
4269 :
4270 76660 : switch (b) {
4271 22897 : case BITS_IN_HALFULONG:
4272 22897 : pack = kron_pack_Flx_spec_half;
4273 22897 : unpack = int_to_Flx_half;
4274 22897 : break;
4275 801 : case BITS_IN_LONG:
4276 801 : pack = kron_pack_Flx_spec;
4277 801 : unpack = kron_unpack_Flx;
4278 801 : break;
4279 49 : case 2*BITS_IN_LONG:
4280 49 : pack = kron_pack_Flx_spec_2;
4281 49 : unpack = kron_unpack_Flx_2;
4282 49 : break;
4283 0 : case 3*BITS_IN_LONG:
4284 0 : pack = kron_pack_Flx_spec_3;
4285 0 : unpack = kron_unpack_Flx_3;
4286 0 : break;
4287 52913 : default:
4288 52913 : A = FlxM_pack_ZM_bits(A, b);
4289 52913 : B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
4290 52913 : C = ZM_mul(A, B);
4291 52913 : D = ZM_unpack_FlxqM_bits(C, b, T, p);
4292 52913 : return gerepilecopy(av, D);
4293 : }
4294 23747 : A = FlxM_pack_ZM(A, pack);
4295 23747 : B = is_sqr? A: FlxM_pack_ZM(B, pack);
4296 23747 : C = ZM_mul(A, B);
4297 23747 : D = ZM_unpack_FlxqM(C, T, p, unpack);
4298 23747 : return gerepilecopy(av, D);
4299 : }
|
