Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_mat
19 :
20 : /********************************************************************/
21 : /** **/
22 : /** REDUCTION **/
23 : /** **/
24 : /********************************************************************/
25 : /* z in Z^n, return lift(Col(z) * Mod(1,p)) */
26 : GEN
27 30021202 : FpC_red(GEN x, GEN p)
28 209968660 : { pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
29 :
30 : /* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
31 : GEN
32 420 : FpV_red(GEN x, GEN p)
33 2695 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
34 :
35 : GEN
36 7004897 : FpC_center(GEN x, GEN p, GEN pov2)
37 42865433 : { pari_APPLY_type(t_COL, Fp_center(gel(x,i), p, pov2)) }
38 :
39 : GEN
40 164710 : Flv_center(GEN x, ulong p, ulong ps2)
41 2621724 : { pari_APPLY_ulong(Fl_center(uel(x,i),p,ps2)) }
42 :
43 : /* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
44 : GEN
45 5616736 : FpM_red(GEN x, GEN p)
46 24448226 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2792834 : FpM_center(GEN x, GEN p, GEN pov2)
50 8370014 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
51 :
52 : /* HACK: assume p > 3, which ensures u = gen_[0-2] is never written to */
53 : static void
54 168 : Fp_center_inplace(GEN u, GEN p, GEN ps2)
55 : {
56 168 : if (abscmpii(u, ps2) > 0)
57 : {
58 28 : pari_sp av = avma;
59 28 : affii(subii(u, p), u);
60 28 : set_avma(av);
61 : }
62 168 : }
63 :
64 : /* p > 3; assume entries in [0,p[ and ps2 = p>>1. */
65 : static void
66 42 : _FpC_center_inplace(GEN z, GEN p, GEN ps2)
67 : {
68 42 : long i, l = lg(z);
69 210 : for (i = 1; i < l; i++) Fp_center_inplace(gel(z,i), p, ps2);
70 42 : }
71 : static void
72 7 : _F3C_center_inplace(GEN z)
73 : {
74 7 : long i, l = lg(z);
75 35 : for (i = 1; i < l; i++) /* z[i] = 0, 1 : no-op */
76 28 : if (equaliu(gel(z,i), 2)) gel(z,i) = gen_m1;
77 7 : }
78 : void
79 0 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
80 : {
81 0 : switch (itou_or_0(p))
82 : {
83 0 : case 2: break;
84 0 : case 3:_F3C_center_inplace(z); break;
85 0 : default: _FpC_center_inplace(z, p, ps2);
86 : }
87 0 : }
88 :
89 : void
90 42 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
91 : {
92 42 : long i, l = lg(z);
93 42 : switch (itou_or_0(p))
94 : {
95 21 : case 2: break;
96 7 : case 3:
97 14 : for (i = 1; i < l; i++) _F3C_center_inplace(gel(z,i));
98 7 : break;
99 14 : default:
100 56 : for (i = 1; i < l; i++) _FpC_center_inplace(gel(z,i), p, pov2);
101 : }
102 42 : }
103 : GEN
104 11347 : Flm_center(GEN x, ulong p, ulong ps2)
105 175875 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
106 :
107 : GEN
108 147 : random_FpV(long d, GEN p)
109 : {
110 : long i;
111 147 : GEN y = cgetg(d+1,t_VEC);
112 23183 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
113 147 : return y;
114 : }
115 :
116 : GEN
117 924 : random_FpC(long d, GEN p)
118 : {
119 : long i;
120 924 : GEN y = cgetg(d+1,t_COL);
121 6188 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
122 924 : return y;
123 : }
124 :
125 : GEN
126 41404 : random_Flv(long n, ulong p)
127 : {
128 41404 : GEN y = cgetg(n+1, t_VECSMALL);
129 : long i;
130 262074 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
131 41406 : return y;
132 : }
133 :
134 : /********************************************************************/
135 : /** **/
136 : /** ADD, SUB **/
137 : /** **/
138 : /********************************************************************/
139 : GEN
140 283453 : FpC_add(GEN x, GEN y, GEN p)
141 3450765 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
142 :
143 : GEN
144 0 : FpV_add(GEN x, GEN y, GEN p)
145 0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
146 :
147 : GEN
148 0 : FpM_add(GEN x, GEN y, GEN p)
149 0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
150 :
151 : GEN
152 1159163 : Flv_add(GEN x, GEN y, ulong p)
153 : {
154 1159163 : if (p==2)
155 2478925 : pari_APPLY_ulong(x[i]^y[i])
156 : else
157 13903675 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
158 : }
159 :
160 : void
161 509440 : Flv_add_inplace(GEN x, GEN y, ulong p)
162 : {
163 509440 : long i, l = lg(x);
164 509440 : if (p==2)
165 1569662 : for (i = 1; i < l; i++) x[i] ^= y[i];
166 : else
167 164878 : for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
168 509440 : }
169 :
170 : ulong
171 3871 : Flv_sum(GEN x, ulong p)
172 : {
173 3871 : long i, l = lg(x);
174 3871 : ulong s = 0;
175 3871 : if (p==2)
176 17689 : for (i = 1; i < l; i++) s ^= x[i];
177 : else
178 0 : for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
179 3871 : return s;
180 : }
181 :
182 : GEN
183 827139 : FpC_sub(GEN x, GEN y, GEN p)
184 11174101 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
185 :
186 : GEN
187 0 : FpV_sub(GEN x, GEN y, GEN p)
188 0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
189 :
190 : GEN
191 0 : FpM_sub(GEN x, GEN y, GEN p)
192 0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
193 :
194 : GEN
195 220682232 : Flv_sub(GEN x, GEN y, ulong p)
196 990181113 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
197 :
198 : void
199 0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
200 : {
201 0 : long i, l = lg(x);
202 0 : for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
203 0 : }
204 :
205 : GEN
206 51293 : Flm_Fl_add(GEN x, ulong y, ulong p)
207 : {
208 51293 : long l = lg(x), i, j;
209 51293 : GEN z = cgetg(l,t_MAT);
210 :
211 51293 : if (l==1) return z;
212 51293 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
213 241346 : for (i=1; i<l; i++)
214 : {
215 190054 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
216 190053 : gel(z,i) = zi;
217 1329706 : for (j=1; j<l; j++) zi[j] = xi[j];
218 190053 : zi[i] = Fl_add(zi[i], y, p);
219 : }
220 51292 : return z;
221 : }
222 : GEN
223 3619 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
224 :
225 : GEN
226 23751 : Flm_add(GEN x, GEN y, ulong p)
227 363189 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
228 :
229 : GEN
230 25574441 : Flm_sub(GEN x, GEN y, ulong p)
231 246236885 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
232 :
233 : /********************************************************************/
234 : /** **/
235 : /** MULTIPLICATION **/
236 : /** **/
237 : /********************************************************************/
238 : GEN
239 964652 : FpC_Fp_mul(GEN x, GEN y, GEN p)
240 11607618 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
241 :
242 : GEN
243 1227048 : Flv_Fl_mul(GEN x, ulong y, ulong p)
244 16903726 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
245 :
246 : GEN
247 4706 : Flv_Fl_div(GEN x, ulong y, ulong p)
248 : {
249 4706 : return Flv_Fl_mul(x, Fl_inv(y, p), p);
250 : }
251 : void
252 0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
253 : {
254 0 : Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
255 0 : }
256 : GEN
257 763925 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
258 763925 : long i, j, h, l = lg(X);
259 763925 : GEN A = cgetg(l, t_MAT);
260 763925 : if (l == 1) return A;
261 763925 : h = lgcols(X);
262 4046149 : for (j=1; j<l; j++)
263 : {
264 3282350 : GEN a = cgetg(h, t_COL), x = gel(X, j);
265 23552338 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
266 3282224 : gel(A,j) = a;
267 : }
268 763799 : return A;
269 : }
270 :
271 : /* x *= y */
272 : void
273 7881363 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
274 : {
275 : long i;
276 7881363 : if (HIGHWORD(y | p))
277 8587572 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
278 : else
279 32186270 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
280 7881360 : }
281 : void
282 0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
283 0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
284 :
285 : /* set y *= x */
286 : void
287 0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
288 : {
289 0 : long i, j, m, l = lg(y);
290 0 : if (l == 1) return;
291 0 : m = lgcols(y);
292 0 : if (HIGHWORD(x | p))
293 0 : for(j=1; j<l; j++)
294 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
295 : else
296 0 : for(j=1; j<l; j++)
297 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 17171790 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
303 : {
304 17171790 : long i, j, m, l = lg(y);
305 17171790 : GEN z = cgetg(l, t_MAT);
306 17170308 : if (l == 1) return z;
307 17170308 : m = lgcols(y);
308 156097493 : for(j=1; j<l; j++) {
309 138841321 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 372116223 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
311 : }
312 17256172 : return z;
313 : }
314 :
315 : /* return x * y */
316 : static GEN
317 8456091 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
318 : {
319 8456091 : long i, j, m, l = lg(y);
320 8456091 : GEN z = cgetg(l, t_MAT);
321 8456083 : if (l == 1) return z;
322 8456083 : m = lgcols(y);
323 48001491 : for(j=1; j<l; j++) {
324 39545412 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
325 127424422 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
326 : }
327 8456079 : return z;
328 : }
329 :
330 : /* return x * y */
331 : GEN
332 25556795 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
333 : {
334 25556795 : if (HIGHWORD(p))
335 17171665 : return Flm_Fl_mul_pre_i(y, x, p, pi);
336 : else
337 8385130 : return Flm_Fl_mul_OK(y, x, p);
338 : }
339 :
340 : /* return x * y */
341 : GEN
342 71006 : Flm_Fl_mul(GEN y, ulong x, ulong p)
343 : {
344 71006 : if (HIGHWORD(p))
345 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
346 : else
347 70932 : return Flm_Fl_mul_OK(y, x, p);
348 : }
349 :
350 : GEN
351 4069689 : Flv_neg(GEN x, ulong p)
352 31168713 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
353 :
354 : void
355 6579 : Flv_neg_inplace(GEN v, ulong p)
356 : {
357 : long i;
358 192965 : for (i = 1; i < lg(v); ++i)
359 186386 : v[i] = Fl_neg(v[i], p);
360 6579 : }
361 :
362 : GEN
363 322716 : Flm_neg(GEN x, ulong p)
364 4361700 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
365 :
366 : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
367 : static GEN
368 18557391 : FpMrow_FpC_mul_i(GEN x, GEN y, long lx, long i, GEN p)
369 : {
370 18557391 : pari_sp av = avma;
371 18557391 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
372 : long k;
373 223190837 : for (k = 2; k < lx; k++)
374 : {
375 204634610 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
376 204632388 : if (signe(t)) c = addii(c, t);
377 : }
378 18556227 : return gc_INT(av, modii(c, p));
379 : }
380 :
381 : static long
382 97688539 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
383 : {
384 : long k;
385 97688539 : long c = coeff(x,i,1) * y[1];
386 1498281551 : for (k = 2; k < lx; k++)
387 1400593012 : c += coeff(x,i,k) * y[k];
388 97688539 : return c;
389 : }
390 :
391 : GEN
392 6450537 : zm_zc_mul(GEN x, GEN y)
393 : {
394 6450537 : long lx = lg(x), l, i;
395 : GEN z;
396 6450537 : if (lx == 1) return cgetg(1, t_VECSMALL);
397 6450537 : l = lg(gel(x,1));
398 6450537 : z = cgetg(l,t_VECSMALL);
399 104139076 : for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
400 6450537 : return z;
401 : }
402 :
403 : GEN
404 2114 : zm_mul(GEN x, GEN y)
405 : {
406 2114 : long i,j,lx=lg(x), ly=lg(y);
407 : GEN z;
408 2114 : if (ly==1) return cgetg(1,t_MAT);
409 2114 : z = cgetg(ly,t_MAT);
410 2114 : if (lx==1)
411 : {
412 0 : for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
413 0 : return z;
414 : }
415 30849 : for (j=1; j<ly; j++)
416 28735 : gel(z,j) = zm_zc_mul(x, gel(y,j));
417 2114 : return z;
418 : }
419 :
420 : static ulong
421 650184270 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
422 : {
423 650184270 : ulong c = ucoeff(x,i,1) * uel(y,1);
424 : long k;
425 7895016357 : for (k = 2; k < lx; k++) {
426 7244832087 : c += ucoeff(x,i,k) * uel(y,k);
427 7244832087 : if (c & HIGHBIT) c %= p;
428 : }
429 650184270 : return c % p;
430 : }
431 :
432 : static ulong
433 537540927 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
434 : {
435 : ulong l0, l1, v1, h0, h1;
436 537540927 : long k = 1;
437 : LOCAL_OVERFLOW;
438 : LOCAL_HIREMAINDER;
439 537540927 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
440 8856935655 : while (++k < lx) {
441 8319394728 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
442 8319394728 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
443 : }
444 537540927 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
445 55258587 : else return remlll_pre(v1, h1, l1, p, pi);
446 : }
447 :
448 : static GEN
449 310747 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
450 : {
451 : long i,j;
452 310747 : GEN z = NULL;
453 :
454 2742128 : for (j=1; j<lx; j++)
455 : {
456 2431382 : if (!y[j]) continue;
457 720987 : if (!z) z = Flv_copy(gel(x,j));
458 14689886 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
459 : }
460 310746 : if (!z) z = zero_zv(l-1);
461 310747 : return z;
462 : }
463 :
464 : static GEN
465 1752784 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
466 : {
467 1752784 : GEN z = cgetg(l,t_COL);
468 : long i;
469 16617491 : for (i = 1; i < l; i++) gel(z,i) = FpMrow_FpC_mul_i(x, y, lx, i, p);
470 1752691 : return z;
471 : }
472 :
473 : static void
474 76556618 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
475 : {
476 : long i;
477 726762670 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
478 76578474 : }
479 : static GEN
480 73336397 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
481 : {
482 73336397 : GEN z = cgetg(l,t_VECSMALL);
483 73335056 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
484 73337399 : return z;
485 : }
486 :
487 : static void
488 95925051 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
489 : {
490 : long i;
491 635516980 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
492 96248131 : }
493 : static GEN
494 91561604 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
495 : {
496 91561604 : GEN z = cgetg(l,t_VECSMALL);
497 91378369 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
498 91613589 : return z;
499 : }
500 :
501 : GEN
502 2004160 : FpM_mul(GEN x, GEN y, GEN p)
503 : {
504 2004160 : long lx=lg(x), ly=lg(y);
505 : GEN z;
506 2004160 : pari_sp av = avma;
507 2004160 : if (ly==1) return cgetg(1,t_MAT);
508 2004160 : if (lx==1) return zeromat(0, ly-1);
509 2004160 : if (lgefint(p) == 3)
510 : {
511 1832482 : ulong pp = uel(p,2);
512 1832482 : if (pp == 2)
513 : {
514 720751 : x = ZM_to_F2m(x);
515 720806 : y = ZM_to_F2m(y);
516 720813 : z = F2m_to_ZM(F2m_mul(x,y));
517 : }
518 : else
519 : {
520 1111731 : x = ZM_to_Flm(x, pp);
521 1111752 : y = ZM_to_Flm(y, pp);
522 1111748 : z = Flm_to_ZM(Flm_mul(x,y, pp));
523 : }
524 : } else
525 171678 : z = FpM_red(ZM_mul(x, y), p);
526 2004236 : return gc_upto(av, z);
527 : }
528 :
529 : static GEN
530 28225482 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
531 : {
532 : long j;
533 28225482 : GEN z = cgetg(ly, t_MAT);
534 28224332 : if (SMALL_ULONG(p))
535 93186188 : for (j=1; j<ly; j++)
536 72997571 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
537 : else
538 99394613 : for (j=1; j<ly; j++)
539 91356686 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
540 28226544 : return z;
541 : }
542 :
543 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
544 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
545 : static void
546 66254 : add_slices_ip(long m, long n,
547 : GEN A, long ma, long da, long na, long ea,
548 : GEN B, long mb, long db, long nb, long eb,
549 : GEN M, long dy, long dx, ulong p)
550 : {
551 66254 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
552 : GEN C;
553 :
554 2780930 : for (j = 1; j <= min_e; j++) {
555 2714674 : C = gel(M, j + dx) + dy;
556 123894325 : for (i = 1; i <= min_d; i++)
557 121179651 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
558 121179679 : ucoeff(B, mb + i, nb + j), p);
559 2779206 : for (; i <= da; i++)
560 64560 : uel(C, i) = ucoeff(A, ma + i, na + j);
561 2714646 : for (; i <= db; i++)
562 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
563 2714646 : for (; i <= m; i++)
564 0 : uel(C, i) = 0;
565 : }
566 69869 : for (; j <= ea; j++) {
567 3613 : C = gel(M, j + dx) + dy;
568 195159 : for (i = 1; i <= da; i++)
569 191546 : uel(C, i) = ucoeff(A, ma + i, na + j);
570 3613 : for (; i <= m; i++)
571 0 : uel(C, i) = 0;
572 : }
573 66256 : for (; j <= eb; j++) {
574 0 : C = gel(M, j + dx) + dy;
575 0 : for (i = 1; i <= db; i++)
576 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
577 0 : for (; i <= m; i++)
578 0 : uel(C, i) = 0;
579 : }
580 66256 : for (; j <= n; j++) {
581 0 : C = gel(M, j + dx) + dy;
582 0 : for (i = 1; i <= m; i++)
583 0 : uel(C, i) = 0;
584 : }
585 66256 : }
586 :
587 : static GEN
588 33128 : add_slices(long m, long n,
589 : GEN A, long ma, long da, long na, long ea,
590 : GEN B, long mb, long db, long nb, long eb, ulong p)
591 : {
592 : GEN M;
593 : long j;
594 33128 : M = cgetg(n + 1, t_MAT);
595 1415027 : for (j = 1; j <= n; j++)
596 1381900 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
597 33127 : add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
598 33128 : return M;
599 : }
600 :
601 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
602 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
603 : static GEN
604 57974 : subtract_slices(long m, long n,
605 : GEN A, long ma, long da, long na, long ea,
606 : GEN B, long mb, long db, long nb, long eb, ulong p)
607 : {
608 57974 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
609 57974 : GEN M = cgetg(n + 1, t_MAT), C;
610 :
611 2526676 : for (j = 1; j <= min_e; j++) {
612 2468702 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
613 116108861 : for (i = 1; i <= min_d; i++)
614 113640164 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
615 113639949 : ucoeff(B, mb + i, nb + j), p);
616 2569360 : for (; i <= da; i++)
617 100448 : uel(C, i) = ucoeff(A, ma + i, na + j);
618 2543952 : for (; i <= db; i++)
619 75040 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
620 2468912 : for (; i <= m; i++)
621 0 : uel(C, i) = 0;
622 : }
623 57974 : for (; j <= ea; j++) {
624 0 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
625 0 : for (i = 1; i <= da; i++)
626 0 : uel(C, i) = ucoeff(A, ma + i, na + j);
627 0 : for (; i <= m; i++)
628 0 : uel(C, i) = 0;
629 : }
630 59393 : for (; j <= eb; j++) {
631 1419 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
632 82481 : for (i = 1; i <= db; i++)
633 81062 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
634 1419 : for (; i <= m; i++)
635 0 : uel(C, i) = 0;
636 : }
637 59393 : for (; j <= n; j++)
638 1419 : gel(M, j) = zero_Flv(m);
639 57974 : return M;
640 : }
641 :
642 : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
643 :
644 : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
645 : static GEN
646 8282 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
647 : {
648 : pari_sp av;
649 8282 : long m1 = (m + 1)/2, m2 = m/2,
650 8282 : n1 = (n + 1)/2, n2 = n/2,
651 8282 : p1 = (p + 1)/2, p2 = p/2;
652 : GEN A11, A12, A22, B11, B21, B22,
653 : S1, S2, S3, S4, T1, T2, T3, T4,
654 : M1, M2, M3, M4, M5, M6, M7,
655 : V1, V2, V3, C;
656 : long j;
657 8282 : C = cgetg(p + 1, t_MAT);
658 676685 : for (j = 1; j <= p; j++)
659 668403 : gel(C, j) = cgetg(m + 1, t_VECSMALL);
660 8282 : av = avma;
661 8282 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
662 8282 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
663 8282 : M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
664 8282 : if (gc_needed(av, 1))
665 0 : (void)gc_all(av, 2, &T2, &M2); /* destroy S1 */
666 8282 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
667 8282 : if (gc_needed(av, 1))
668 0 : (void)gc_all(av, 2, &M2, &T3); /* destroy T2 */
669 8282 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
670 8282 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
671 8282 : M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
672 8282 : if (gc_needed(av, 1))
673 0 : (void)gc_all(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
674 8282 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
675 8282 : if (gc_needed(av, 1))
676 0 : (void)gc_all(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
677 8282 : A11 = matslice(A, 1, m1, 1, n1);
678 8282 : B11 = matslice(B, 1, n1, 1, p1);
679 8282 : M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
680 8282 : if (gc_needed(av, 1))
681 0 : (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
682 8282 : A12 = matslice(A, 1, m1, n1 + 1, n);
683 8282 : B21 = matslice(B, n1 + 1, n, 1, p1);
684 8282 : M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
685 8282 : if (gc_needed(av, 1))
686 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
687 8282 : add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
688 8282 : if (gc_needed(av, 1))
689 0 : (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy M4 */
690 8282 : M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
691 8282 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
692 8282 : if (gc_needed(av, 1))
693 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &M1, &M5, &S4); /* destroy S3 */
694 8282 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
695 8282 : if (gc_needed(av, 1))
696 0 : (void)gc_all(av, 6, &M2, &M3, &M1, &M5, &S4, &T4); /* destroy T3 */
697 8282 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
698 8282 : if (gc_needed(av, 1))
699 0 : (void)gc_all(av, 5, &M2, &M3, &S4, &T4, &V1); /* destroy M1, M5 */
700 8282 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
701 8282 : M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
702 8282 : if (gc_needed(av, 1))
703 0 : (void)gc_all(av, 5, &M2, &M3, &T4, &V1, &M6); /* destroy S4, B22 */
704 8282 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
705 8282 : M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
706 8282 : if (gc_needed(av, 1))
707 0 : (void)gc_all(av, 5, &M2, &M3, &V1, &M6, &M7); /* destroy A22, T4 */
708 8282 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
709 8282 : add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
710 8282 : if (gc_needed(av, 1))
711 0 : (void)gc_all(av, 4, &M2, &M3, &V1, &M7); /* destroy V3, M6 */
712 8282 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
713 8282 : if (gc_needed(av, 1))
714 0 : (void)gc_all(av, 3, &M3, &M7, &V2); /* destroy V1, M2 */
715 8282 : add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
716 8282 : add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
717 8282 : set_avma(av); return C;
718 : }
719 :
720 : /* Strassen-Winograd used for dim >= ZM_sw_bound */
721 : static GEN
722 28232366 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
723 : {
724 28232366 : ulong e = expu(p);
725 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
726 23699382 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
727 : #else
728 4533444 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
729 : #endif
730 28232826 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
731 28224544 : return Flm_mul_classical(x, y, l, lx, ly, p, pi);
732 : else
733 8282 : return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
734 : }
735 :
736 : GEN
737 13663898 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
738 : {
739 13663898 : long lx=lg(x), ly=lg(y);
740 13663898 : if (ly==1) return cgetg(1,t_MAT);
741 13663898 : if (lx==1) return zero_Flm(0, ly-1);
742 13663898 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
743 : }
744 :
745 : GEN
746 14430356 : Flm_mul(GEN x, GEN y, ulong p)
747 : {
748 14430356 : long lx=lg(x), ly=lg(y);
749 14430356 : if (ly==1) return cgetg(1,t_MAT);
750 14430069 : if (lx==1) return zero_Flm(0, ly-1);
751 14430069 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
752 : }
753 :
754 : GEN
755 79186 : Flm_sqr(GEN x, ulong p)
756 : {
757 79186 : long lx = lg(x);
758 79186 : if (lx==1) return cgetg(1,t_MAT);
759 79186 : return Flm_mul_i(x, x, lx, lx, lx, p, get_Fl_red(p));
760 : }
761 :
762 : struct _Flm
763 : {
764 : ulong p;
765 : long n;
766 : };
767 :
768 : static GEN
769 13861 : _Flm_mul(void *E , GEN x, GEN y)
770 13861 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
771 : static GEN
772 79186 : _Flm_sqr(void *E, GEN x)
773 79186 : { return Flm_sqr(x,((struct _Flm*)E)->p); }
774 : static GEN
775 826 : _Flm_one(void *E)
776 826 : { return matid_Flm(((struct _Flm*)E)->n); }
777 : GEN
778 45309 : Flm_powu(GEN x, ulong n, ulong p)
779 : {
780 : struct _Flm d;
781 45309 : if (!n) return matid(lg(x)-1);
782 45309 : d.p = p;
783 45309 : return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
784 : }
785 : GEN
786 826 : Flm_powers(GEN x, ulong n, ulong p)
787 : {
788 : struct _Flm d;
789 826 : d.p = p;
790 826 : d.n = lg(x)-1;
791 826 : return gen_powers(x, n, 1, (void*)&d,
792 : &_Flm_sqr, &_Flm_mul, &_Flm_one);
793 : }
794 :
795 : static GEN
796 0 : _FpM_mul(void *p , GEN x, GEN y)
797 0 : { return FpM_mul(x,y,(GEN)p); }
798 : static GEN
799 0 : _FpM_sqr(void *p, GEN x)
800 0 : { return FpM_mul(x,x,(GEN)p); }
801 : GEN
802 0 : FpM_powu(GEN x, ulong n, GEN p)
803 : {
804 0 : if (!n) return matid(lg(x)-1);
805 0 : if (lgefint(p) == 3)
806 : {
807 0 : pari_sp av = avma;
808 0 : ulong pp = uel(p,2);
809 : GEN z;
810 0 : if (pp == 2)
811 0 : z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
812 : else
813 0 : z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
814 0 : return gc_upto(av, z);
815 : }
816 0 : return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
817 : }
818 :
819 : /*Multiple a column vector by a line vector to make a matrix*/
820 : GEN
821 14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
822 : {
823 14 : long i,j, lx=lg(x), ly=lg(y);
824 : GEN z;
825 14 : if (ly==1) return cgetg(1,t_MAT);
826 14 : z = cgetg(ly, t_MAT);
827 49 : for (j=1; j < ly; j++)
828 : {
829 35 : GEN zj = cgetg(lx,t_VECSMALL);
830 112 : for (i=1; i<lx; i++)
831 77 : uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
832 35 : gel(z,j) = zj;
833 : }
834 14 : return z;
835 : }
836 :
837 : /*Multiple a column vector by a line vector to make a matrix*/
838 : GEN
839 5777 : FpC_FpV_mul(GEN x, GEN y, GEN p)
840 : {
841 5777 : long i,j, lx=lg(x), ly=lg(y);
842 : GEN z;
843 5777 : if (ly==1) return cgetg(1,t_MAT);
844 5777 : z = cgetg(ly,t_MAT);
845 63885 : for (j=1; j < ly; j++)
846 : {
847 58108 : GEN zj = cgetg(lx,t_COL);
848 1371116 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
849 58108 : gel(z, j) = zj;
850 : }
851 5777 : return z;
852 : }
853 :
854 : /* Multiply a line vector by a column and return a scalar (t_INT) */
855 : GEN
856 8491009 : FpV_dotproduct(GEN x, GEN y, GEN p)
857 : {
858 8491009 : long i, lx = lg(x);
859 : pari_sp av;
860 : GEN c;
861 8491009 : if (lx == 1) return gen_0;
862 8488496 : av = avma; c = mulii(gel(x,1),gel(y,1));
863 30285855 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
864 8488433 : return gc_INT(av, modii(c,p));
865 : }
866 : GEN
867 0 : FpV_dotsquare(GEN x, GEN p)
868 : {
869 0 : long i, lx = lg(x);
870 : pari_sp av;
871 : GEN c;
872 0 : if (lx == 1) return gen_0;
873 0 : av = avma; c = sqri(gel(x,1));
874 0 : for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
875 0 : return gc_INT(av, modii(c,p));
876 : }
877 :
878 : INLINE ulong
879 9383092 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
880 : {
881 9383092 : ulong c = uel(x,0) * uel(y,0);
882 : long k;
883 66488802 : for (k = 1; k < lx; k++) {
884 57105710 : c += uel(x,k) * uel(y,k);
885 57105710 : if (c & HIGHBIT) c %= p;
886 : }
887 9383092 : return c % p;
888 : }
889 :
890 : INLINE ulong
891 748114 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
892 : {
893 : ulong l0, l1, v1, h0, h1;
894 748114 : long i = 0;
895 : LOCAL_OVERFLOW;
896 : LOCAL_HIREMAINDER;
897 748114 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
898 15480608 : while (++i < lx) {
899 14732494 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
900 14732494 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
901 : }
902 748114 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
903 464890 : else return remlll_pre(v1, h1, l1, p, pi);
904 : }
905 :
906 : ulong
907 531366 : Flv_dotproduct(GEN x, GEN y, ulong p)
908 : {
909 531366 : long lx = lg(x)-1;
910 531366 : if (lx == 0) return 0;
911 531366 : if (SMALL_ULONG(p))
912 531366 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
913 : else
914 0 : return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
915 : }
916 :
917 : ulong
918 66256 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
919 : {
920 66256 : long lx = lg(x)-1;
921 66256 : if (lx == 0) return 0;
922 66256 : if (SMALL_ULONG(p))
923 60044 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
924 : else
925 6212 : return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
926 : }
927 :
928 : ulong
929 9601223 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
930 : {
931 9601223 : long lx = minss(lgpol(x), lgpol(y));
932 9601219 : if (lx == 0) return 0;
933 9533584 : if (pi)
934 741902 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
935 : else
936 8791682 : return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
937 : }
938 : ulong
939 0 : Flx_dotproduct(GEN x, GEN y, ulong p)
940 0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
941 :
942 : GEN
943 1752784 : FpM_FpC_mul(GEN x, GEN y, GEN p)
944 : {
945 1752784 : long lx = lg(x);
946 1752784 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
947 : }
948 : GEN
949 887025 : Flm_Flc_mul(GEN x, GEN y, ulong p)
950 : {
951 887025 : long l, lx = lg(x);
952 887025 : if (lx==1) return cgetg(1,t_VECSMALL);
953 887025 : l = lgcols(x);
954 887025 : if (p==2)
955 310747 : return Flm_Flc_mul_i_2(x, y, lx, l);
956 576278 : else if (SMALL_ULONG(p))
957 338938 : return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
958 : else
959 237340 : return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
960 : }
961 :
962 : /* allow pi = 0 */
963 : GEN
964 6791 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
965 : {
966 6791 : long l, lx = lg(x);
967 : GEN z;
968 6791 : if (lx==1) return cgetg(1,t_VECSMALL);
969 6791 : l = lgcols(x);
970 6791 : z = cgetg(l, t_VECSMALL);
971 6791 : if (SMALL_ULONG(p))
972 4407 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
973 : else
974 2384 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
975 6791 : return z;
976 : }
977 :
978 : /* allow pi = 0 */
979 : GEN
980 7730042 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
981 : {
982 7730042 : long l, lx = lg(x);
983 : GEN z;
984 7730042 : if (lx==1) return pol0_Flx(sv);
985 7730042 : l = lgcols(x);
986 7728189 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
987 7726579 : if (SMALL_ULONG(p))
988 3216184 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
989 : else
990 4510395 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
991 7741836 : return Flx_renormalize(z, l + 1);
992 : }
993 :
994 : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
995 : GEN
996 1154696 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
997 : {
998 1154696 : long i, l, lx = lg(x);
999 : GEN z;
1000 1154696 : if (lx==1) return pol_0(v);
1001 1154696 : l = lgcols(x);
1002 1154696 : z = new_chunk(l+1);
1003 1509855 : for (i=l-1; i; i--)
1004 : {
1005 1494990 : pari_sp av = avma;
1006 1494990 : GEN t = FpMrow_FpC_mul_i(x,y,lx,i,p);
1007 1494987 : if (signe(t))
1008 : {
1009 1139830 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
1010 1139830 : gel(z,i+1) = t;
1011 1139830 : break;
1012 : }
1013 355157 : set_avma(av);
1014 : }
1015 1154695 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1016 1139830 : z[0] = evaltyp(t_POL) | _evallg(i+2);
1017 1139830 : z[1] = evalsigne(1) | evalvarn(v);
1018 3337432 : for (; i; i--) gel(z,i+1) = FpMrow_FpC_mul_i(x,y,lx,i,p);
1019 1139834 : return z;
1020 : }
1021 :
1022 : /********************************************************************/
1023 : /** **/
1024 : /** TRANSPOSITION **/
1025 : /** **/
1026 : /********************************************************************/
1027 :
1028 : /* == zm_transpose */
1029 : GEN
1030 730726 : Flm_transpose(GEN x)
1031 : {
1032 730726 : long i, dx, lx = lg(x);
1033 : GEN y;
1034 730726 : if (lx == 1) return cgetg(1,t_MAT);
1035 730593 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1036 8830557 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1037 730592 : return y;
1038 : }
1039 :
1040 : /********************************************************************/
1041 : /** **/
1042 : /** SCALAR MATRICES **/
1043 : /** **/
1044 : /********************************************************************/
1045 :
1046 : GEN
1047 2009 : gen_matid(long n, void *E, const struct bb_field *S)
1048 : {
1049 2009 : GEN y = cgetg(n+1,t_MAT), _0, _1;
1050 : long i;
1051 2009 : if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
1052 2009 : _0 = S->s(E,0);
1053 2009 : _1 = S->s(E,1);
1054 8883 : for (i=1; i<=n; i++)
1055 : {
1056 6874 : GEN z = const_col(n, _0); gel(z,i) = _1;
1057 6874 : gel(y, i) = z;
1058 : }
1059 2009 : return y;
1060 : }
1061 :
1062 : GEN
1063 35 : matid_F2xqM(long n, GEN T)
1064 : {
1065 : void *E;
1066 35 : const struct bb_field *S = get_F2xq_field(&E, T);
1067 35 : return gen_matid(n, E, S);
1068 : }
1069 : GEN
1070 7 : matid_FlxqM(long n, GEN T, ulong p)
1071 : {
1072 : void *E;
1073 7 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1074 7 : return gen_matid(n, E, S);
1075 : }
1076 :
1077 : GEN
1078 1435516 : matid_Flm(long n)
1079 : {
1080 1435516 : GEN y = cgetg(n+1,t_MAT);
1081 : long i;
1082 1435511 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1083 9132861 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1084 1435490 : return y;
1085 : }
1086 :
1087 : GEN
1088 42 : scalar_Flm(long s, long n)
1089 : {
1090 : long i;
1091 42 : GEN y = cgetg(n+1,t_MAT);
1092 560 : for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
1093 42 : return y;
1094 : }
1095 :
1096 : /********************************************************************/
1097 : /** **/
1098 : /** CONVERSIONS **/
1099 : /** **/
1100 : /********************************************************************/
1101 : GEN
1102 54925059 : ZV_to_Flv(GEN x, ulong p)
1103 560007963 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1104 :
1105 : GEN
1106 13588905 : ZM_to_Flm(GEN x, ulong p)
1107 67287552 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1108 :
1109 : GEN
1110 14722 : ZVV_to_FlvV(GEN x, ulong p)
1111 112431 : { pari_APPLY_type(t_VEC, ZV_to_Flv(gel(x,i), p)) }
1112 :
1113 : GEN
1114 4780 : ZMV_to_FlmV(GEN x, ulong m)
1115 37006 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1116 :
1117 : /* TO INTMOD */
1118 : static GEN
1119 18792807 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1120 : static GEN
1121 708315 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1122 :
1123 : GEN
1124 4003 : Fp_to_mod(GEN z, GEN p)
1125 : {
1126 4003 : retmkintmod(modii(z, p), icopy(p));
1127 : }
1128 :
1129 : GEN
1130 966014 : FpX_to_mod_raw(GEN z, GEN p)
1131 : {
1132 966014 : long i, l = lg(z);
1133 : GEN x;
1134 966014 : if (l == 2)
1135 : {
1136 510440 : x = cgetg(3,t_POL); x[1] = z[1];
1137 510440 : gel(x,2) = mkintmod(gen_0,p); return x;
1138 : }
1139 455574 : x = cgetg(l,t_POL); x[1] = z[1];
1140 3824884 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1141 455574 : return normalizepol_lg(x,l);
1142 : }
1143 :
1144 : /* z in Z[X], return z * Mod(1,p), normalized*/
1145 : GEN
1146 1232742 : FpX_to_mod(GEN z, GEN p)
1147 : {
1148 1232742 : long i, l = lg(z);
1149 : GEN x;
1150 1232742 : if (l == 2)
1151 : {
1152 1197 : x = cgetg(3,t_POL); x[1] = z[1];
1153 1197 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1154 : }
1155 1231545 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1156 16551661 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1157 1231545 : return normalizepol_lg(x,l);
1158 : }
1159 :
1160 : GEN
1161 84021 : FpXC_to_mod(GEN z, GEN p)
1162 : {
1163 84021 : long i,l = lg(z);
1164 84021 : GEN x = cgetg(l, t_COL);
1165 84021 : p = icopy(p);
1166 474159 : for (i=1; i<l; i++)
1167 390138 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1168 84021 : return x;
1169 : }
1170 :
1171 : static GEN
1172 0 : FpXC_to_mod_raw(GEN x, GEN p)
1173 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1174 :
1175 : GEN
1176 0 : FpXM_to_mod(GEN z, GEN p)
1177 : {
1178 0 : long i,l = lg(z);
1179 0 : GEN x = cgetg(l, t_MAT);
1180 0 : p = icopy(p);
1181 0 : for (i=1; i<l; i++)
1182 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1183 0 : return x;
1184 : }
1185 :
1186 : /* z in Z^n, return z * Mod(1,p), normalized*/
1187 : GEN
1188 36180 : FpV_to_mod(GEN z, GEN p)
1189 : {
1190 36180 : long i,l = lg(z);
1191 36180 : GEN x = cgetg(l, t_VEC);
1192 36180 : if (l == 1) return x;
1193 36180 : p = icopy(p);
1194 72955 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1195 36180 : return x;
1196 : }
1197 : /* z in Z^n, return z * Mod(1,p), normalized*/
1198 : GEN
1199 105 : FpC_to_mod(GEN z, GEN p)
1200 : {
1201 105 : long i, l = lg(z);
1202 105 : GEN x = cgetg(l, t_COL);
1203 105 : if (l == 1) return x;
1204 91 : p = icopy(p);
1205 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1206 91 : return x;
1207 : }
1208 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1209 : GEN
1210 226 : FpM_to_mod(GEN z, GEN p)
1211 : {
1212 226 : long i, j, m, l = lg(z);
1213 226 : GEN x = cgetg(l,t_MAT), y, zi;
1214 226 : if (l == 1) return x;
1215 205 : m = lgcols(z);
1216 205 : p = icopy(p);
1217 2456 : for (i=1; i<l; i++)
1218 : {
1219 2251 : gel(x,i) = cgetg(m,t_COL);
1220 2251 : y = gel(x,i); zi= gel(z,i);
1221 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1222 : }
1223 205 : return x;
1224 : }
1225 : GEN
1226 28 : Flc_to_mod(GEN z, ulong pp)
1227 : {
1228 28 : long i, l = lg(z);
1229 28 : GEN p, x = cgetg(l, t_COL);
1230 28 : if (l == 1) return x;
1231 28 : p = utoipos(pp);
1232 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1233 28 : return x;
1234 : }
1235 : GEN
1236 60324 : Flm_to_mod(GEN z, ulong pp)
1237 : {
1238 60324 : long i, j, m, l = lg(z);
1239 60324 : GEN p, x = cgetg(l,t_MAT), y, zi;
1240 60324 : if (l == 1) return x;
1241 60296 : m = lgcols(z);
1242 60296 : p = utoipos(pp);
1243 209518 : for (i=1; i<l; i++)
1244 : {
1245 149222 : gel(x,i) = cgetg(m,t_COL);
1246 149222 : y = gel(x,i); zi= gel(z,i);
1247 857453 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1248 : }
1249 60296 : return x;
1250 : }
1251 :
1252 : GEN
1253 574 : FpVV_to_mod(GEN z, GEN p)
1254 : {
1255 574 : long i, j, m, l = lg(z);
1256 574 : GEN x = cgetg(l,t_VEC), y, zi;
1257 574 : if (l == 1) return x;
1258 574 : m = lgcols(z);
1259 574 : p = icopy(p);
1260 1246 : for (i=1; i<l; i++)
1261 : {
1262 672 : gel(x,i) = cgetg(m,t_VEC);
1263 672 : y = gel(x,i); zi= gel(z,i);
1264 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1265 : }
1266 574 : return x;
1267 : }
1268 :
1269 : /* z in Z^n, return z * Mod(1,p), normalized*/
1270 : GEN
1271 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1272 : {
1273 7 : long i,l = lg(z);
1274 7 : GEN x = cgetg(l, t_COL);
1275 7 : if (l == 1) return x;
1276 7 : p = icopy(p);
1277 7 : T = FpX_to_mod_raw(T, p);
1278 21 : for (i=1; i<l; i++)
1279 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1280 7 : return x;
1281 : }
1282 :
1283 : static GEN
1284 551383 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1285 : {
1286 551383 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1287 551383 : return mkpolmod(z, T);
1288 : }
1289 :
1290 : /* z in Z^n, return z * Mod(1,p), normalized*/
1291 : GEN
1292 14 : FqC_to_mod(GEN z, GEN T, GEN p)
1293 : {
1294 : GEN x;
1295 14 : long i,l = lg(z);
1296 14 : if (!T) return FpC_to_mod(z, p);
1297 14 : x = cgetg(l, t_COL);
1298 14 : if (l == 1) return x;
1299 14 : p = icopy(p);
1300 14 : T = FpX_to_mod_raw(T, p);
1301 252 : for (i=1; i<l; i++)
1302 238 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1303 14 : return x;
1304 : }
1305 :
1306 : /* z in Z^n, return z * Mod(1,p), normalized*/
1307 : static GEN
1308 106631 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1309 657776 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1310 :
1311 : /* z in Z^n, return z * Mod(1,p), normalized*/
1312 : GEN
1313 26040 : FqM_to_mod(GEN z, GEN T, GEN p)
1314 : {
1315 : GEN x;
1316 26040 : long i,l = lg(z);
1317 26040 : if (!T) return FpM_to_mod(z, p);
1318 26040 : x = cgetg(l, t_MAT);
1319 26040 : if (l == 1) return x;
1320 26026 : p = icopy(p);
1321 26026 : T = FpX_to_mod_raw(T, p);
1322 132657 : for (i=1; i<l; i++)
1323 106631 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1324 26026 : return x;
1325 : }
1326 :
1327 : /********************************************************************/
1328 : /* */
1329 : /* Blackbox linear algebra */
1330 : /* */
1331 : /********************************************************************/
1332 :
1333 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1334 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1335 : * (e_j) is the canonical basis.
1336 : * A sparse matrix (zMs) is a t_VEC of zCs */
1337 :
1338 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1339 : * integer representing an element of Fp. This is important since p can be
1340 : * large and p+E[i] would not fit in a C long. */
1341 :
1342 : /* RgCs and RgMs are similar, except that the type of the component is
1343 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1344 : * of E. */
1345 :
1346 : /* Most functions take an argument nbrow which is the number of lines of the
1347 : * column/matrix, which cannot be derived from the data. */
1348 :
1349 : GEN
1350 0 : zCs_to_ZC(GEN R, long nbrow)
1351 : {
1352 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1353 0 : long j, l = lg(C);
1354 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1355 0 : return c;
1356 : }
1357 :
1358 : GEN
1359 0 : zMs_to_ZM(GEN x, long nbrow)
1360 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1361 :
1362 : GEN
1363 0 : zMs_compress(GEN M, long nbrow)
1364 : {
1365 0 : pari_sp av = avma;
1366 : GEN p, q, W;
1367 0 : long c, i, l = lg(M), n;
1368 0 : GEN v = zero_F2v(nbrow);
1369 0 : for (i = 1; i < l; i++)
1370 : {
1371 0 : GEN Vi = gmael(M,i,1);
1372 0 : long li = lg(Vi), j;
1373 0 : for (j = 1; j < li; j++)
1374 0 : F2v_set(v, uel(Vi,j));
1375 : }
1376 0 : n = F2v_hamming(v);
1377 0 : p = cgetg(n+1, t_VECSMALL);
1378 0 : q = zero_Flv(nbrow);
1379 0 : for (c = 1, i = 1; i <= nbrow; i++)
1380 0 : if (F2v_coeff(v,i))
1381 : {
1382 0 : q[i] = c;
1383 0 : p[c++] = i;
1384 : }
1385 0 : W = cgetg(l, t_MAT);
1386 0 : for (i = 1; i < l; i++)
1387 : {
1388 0 : GEN Vi = gmael(M,i,1), Wi;
1389 0 : long li = lg(Vi), j;
1390 0 : Wi = cgetg(li, t_VECSMALL);
1391 0 : for (j = 1; j < li; j++)
1392 0 : uel(Wi, j) = uel(q,uel(Vi,j));
1393 0 : gel(W, i) = mkmat2(Wi, gmael(M,i,2));
1394 : }
1395 0 : return gc_GEN(av, mkvec2(W, p));
1396 : }
1397 :
1398 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1399 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1400 : GEN
1401 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1402 : {
1403 147 : pari_sp ltop = avma;
1404 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1405 147 : if (ZV_equal0(B)) return zerocol(n);
1406 147 : while (++col <= n)
1407 : {
1408 147 : pari_sp btop = avma, av;
1409 : long i, lQ;
1410 147 : GEN V, Q, M, W = B;
1411 147 : GEN b = cgetg(m+2, t_POL);
1412 147 : b[1] = evalsigne(1)|evalvarn(0);
1413 147 : gel(b, 2) = gel(W, col);
1414 46219 : for (i = 3; i<m+2; i++)
1415 46072 : gel(b, i) = cgeti(lgefint(p));
1416 147 : av = avma;
1417 46219 : for (i = 3; i<m+2; i++)
1418 : {
1419 46072 : W = f(E, W);
1420 46072 : affii(gel(W, col),gel(b, i));
1421 46072 : if (gc_needed(av,1))
1422 : {
1423 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1424 72 : W = gc_upto(av, W);
1425 : }
1426 : }
1427 147 : b = FpX_renormalize(b, m+2);
1428 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1429 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1430 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1431 147 : W = B; lQ =lg(Q);
1432 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1433 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1434 147 : av = avma;
1435 22742 : for (i = lQ-3; i > 1; i--)
1436 : {
1437 22595 : W = f(E, W);
1438 22595 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1439 22595 : if (gc_needed(av,1))
1440 : {
1441 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1442 110 : (void)gc_all(av, 2, &V, &W);
1443 : }
1444 : }
1445 147 : V = FpC_red(V, p);
1446 147 : W = FpC_sub(f(E,V), B, p);
1447 147 : if (ZV_equal0(W)) return gc_GEN(ltop, V);
1448 0 : av = avma;
1449 0 : for (i = 1; i <= n; ++i)
1450 : {
1451 0 : V = W;
1452 0 : W = f(E, V);
1453 0 : if (ZV_equal0(W))
1454 0 : return gc_GEN(ltop, shallowtrans(V));
1455 0 : (void)gc_all(av, 2, &V, &W);
1456 : }
1457 0 : set_avma(btop);
1458 : }
1459 0 : return NULL;
1460 : }
1461 :
1462 : GEN
1463 0 : zMs_ZC_mul(GEN M, GEN B)
1464 : {
1465 : long i, j;
1466 0 : long n = lg(B)-1;
1467 0 : GEN V = zerocol(n);
1468 0 : for (i = 1; i <= n; ++i)
1469 0 : if (signe(gel(B, i)))
1470 : {
1471 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1472 0 : long l = lg(C);
1473 0 : for (j = 1; j < l; ++j)
1474 : {
1475 0 : long k = C[j];
1476 0 : switch(E[j])
1477 : {
1478 0 : case 1:
1479 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1480 0 : break;
1481 0 : case -1:
1482 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1483 0 : break;
1484 0 : default:
1485 0 : gel(V, k) = gel(V,k)==gen_0 ? mulis(gel(B, i), E[j]) : addii(gel(V, k), mulis(gel(B, i), E[j]));
1486 0 : break;
1487 : }
1488 : }
1489 : }
1490 0 : return V;
1491 : }
1492 :
1493 : GEN
1494 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1495 :
1496 : GEN
1497 68961 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1498 : {
1499 68961 : long i, j, lM = lg(M);
1500 68961 : GEN V = cgetg(lM,t_VEC);
1501 28331907 : for (i = 1; i < lM; ++i)
1502 : {
1503 28262946 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1504 28262946 : pari_sp av = avma;
1505 28262946 : long lC = lg(C);
1506 28262946 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1507 28262946 : z = mulis(gel(B, C[1]), E[1]);
1508 228726530 : for (j = 2; j < lC; ++j)
1509 : {
1510 200463584 : GEN b = gel(B, C[j]);
1511 200463584 : switch(E[j])
1512 : {
1513 131566519 : case 1: z = addii(z, b); break;
1514 57522273 : case -1: z = subii(z, b); break;
1515 11374792 : default: z = addii(z, mulis(b, E[j])); break;
1516 : }
1517 : }
1518 28262946 : gel(V,i) = gc_INT(av, p? Fp_red(z, p): z);
1519 : }
1520 68961 : return V;
1521 : }
1522 : GEN
1523 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1524 :
1525 : GEN
1526 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1527 : {
1528 0 : pari_sp av = avma, av2;
1529 0 : GEN xi, xb, pi = gen_1, P;
1530 : long i;
1531 0 : if (!C) {
1532 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1533 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1534 : }
1535 0 : P = utoipos(p);
1536 0 : av2 = avma;
1537 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1538 0 : xb = Flm_to_ZM(xi);
1539 0 : for (i = 2; i <= e; i++)
1540 : {
1541 0 : pi = muliu(pi, p); /* = p^(i-1) */
1542 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1543 0 : if (gc_needed(av,2))
1544 : {
1545 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1546 0 : (void)gc_all(av2,3, &pi,&b,&xb);
1547 : }
1548 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1549 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1550 : }
1551 0 : return gc_upto(av, xb);
1552 : }
1553 :
1554 : struct wrapper_modp_s {
1555 : GEN (*f)(void*, GEN);
1556 : void *E;
1557 : GEN p;
1558 : };
1559 :
1560 : /* compute f(x) mod p */
1561 : static GEN
1562 0 : wrap_relcomb_modp(void *E, GEN x)
1563 : {
1564 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1565 0 : return FpC_red(W->f(W->E, x), W->p);
1566 : }
1567 : static GEN
1568 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1569 :
1570 : static GEN
1571 68814 : wrap_relker(void*E, GEN x)
1572 : {
1573 68814 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1574 68814 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1575 : }
1576 :
1577 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1578 : GEN
1579 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1580 : {
1581 : struct wrapper_modp_s W;
1582 0 : pari_sp av = avma;
1583 0 : GEN xb, xi, pi = gen_1;
1584 : long i;
1585 0 : W.E = E;
1586 0 : W.f = f;
1587 0 : W.p = p;
1588 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1589 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1590 0 : xb = xi;
1591 0 : for (i = 2; i <= e; i++)
1592 : {
1593 0 : pi = mulii(pi, p); /* = p^(i-1) */
1594 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1595 0 : if (gc_needed(av,2))
1596 : {
1597 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1598 0 : (void)gc_all(av,3, &pi,&B,&xb);
1599 : }
1600 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1601 0 : if (!xi) return NULL;
1602 0 : if (typ(xi) == t_VEC) return gc_upto(av, xi);
1603 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1604 : }
1605 0 : return gc_upto(av, xb);
1606 : }
1607 :
1608 : static GEN
1609 23036 : vecprow(GEN A, GEN prow)
1610 : {
1611 23036 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1612 : }
1613 :
1614 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1615 : * or the index of a column which is linearly dependent from the others as a
1616 : * t_VECSMALL with a single component. */
1617 : GEN
1618 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1619 : {
1620 0 : pari_sp av = avma;
1621 : GEN pcol, prow;
1622 0 : long nbi=lg(M)-1, lR;
1623 : long i, n;
1624 : GEN Mp, Ap, Rp;
1625 : pari_timer ti;
1626 0 : if (DEBUGLEVEL) timer_start(&ti);
1627 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1628 0 : if (!pcol) return gc_NULL(av);
1629 0 : if (DEBUGLEVEL)
1630 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1631 0 : n = lg(pcol)-1;
1632 0 : Mp = cgetg(n+1, t_MAT);
1633 0 : for(i=1; i<=n; i++)
1634 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1635 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1636 0 : if (DEBUGLEVEL) timer_start(&ti);
1637 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1638 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1639 0 : if (!Rp) return gc_NULL(av);
1640 0 : lR = lg(Rp)-1;
1641 0 : if (typ(Rp) == t_COL)
1642 : {
1643 0 : GEN R = zerocol(nbi+1);
1644 0 : for (i=1; i<=lR; i++)
1645 0 : gel(R,pcol[i]) = gel(Rp,i);
1646 0 : return gc_GEN(av, R);
1647 : }
1648 0 : for (i = 1; i <= lR; ++i)
1649 0 : if (signe(gel(Rp, i)))
1650 0 : return gc_leaf(av, mkvecsmall(pcol[i]));
1651 : return NULL; /* LCOV_EXCL_LINE */
1652 : }
1653 :
1654 : GEN
1655 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1656 : {
1657 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1658 : }
1659 :
1660 : GEN
1661 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1662 : {
1663 : GEN res;
1664 0 : pari_CATCH(e_INV)
1665 : {
1666 0 : GEN E = pari_err_last();
1667 0 : GEN x = err_get_compo(E,2);
1668 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1669 0 : if (DEBUGLEVEL)
1670 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1671 0 : res = NULL;
1672 : } pari_TRY {
1673 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1674 0 : } pari_ENDCATCH
1675 0 : return res;
1676 : }
1677 :
1678 : static GEN
1679 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1680 : {
1681 147 : long i, j, oldf = 0, f = 0;
1682 147 : long lrow = lg(prow), lM = lg(M);
1683 147 : GEN W = const_vecsmall(lM-1,1);
1684 147 : GEN R = cgetg(lrow, t_VEC);
1685 228354 : for (i=1; i<lrow; i++)
1686 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1687 : do
1688 : {
1689 442 : oldf = f;
1690 374995 : for(i=1; i<lM; i++)
1691 374553 : if (W[i])
1692 : {
1693 131729 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1694 131729 : long c=0, cj=0, lF = lg(F);
1695 1093639 : for(j=1; j<lF; j++)
1696 961910 : if (!gel(R,F[j])) { c++; cj=j; }
1697 131729 : if (c>=2) continue;
1698 112311 : if (c==1)
1699 : {
1700 32899 : pari_sp av = avma;
1701 32899 : GEN s = gen_0;
1702 272954 : for(j=1; j<lF; j++)
1703 240055 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1704 : /* s /= -E[cj] mod p */
1705 32899 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1706 32899 : gel(R,F[cj]) = gc_INT(av, s);
1707 32899 : f++;
1708 : }
1709 112311 : W[i]=0;
1710 : }
1711 442 : } while(oldf!=f);
1712 228354 : for (i=1; i<lrow; i++)
1713 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1714 147 : return R;
1715 : }
1716 :
1717 : /* Return a linear form Y such that YM is essentially 0 */
1718 : GEN
1719 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1720 : {
1721 147 : pari_sp av = avma, av2;
1722 : GEN Mp, pcol, prow;
1723 : long i, n;
1724 : pari_timer ti;
1725 : struct wrapper_modp_s W;
1726 147 : if (DEBUGLEVEL) timer_start(&ti);
1727 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1728 147 : if (!pcol) return gc_NULL(av);
1729 147 : if (DEBUGLEVEL)
1730 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1731 147 : n = lg(pcol)-1;
1732 147 : Mp = cgetg(n+1, t_MAT);
1733 23183 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1734 147 : W.E = (void*)Mp;
1735 147 : W.p = p;
1736 147 : av2 = avma;
1737 0 : for(;; set_avma(av2))
1738 0 : {
1739 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1740 147 : if (DEBUGLEVEL) timer_start(&ti);
1741 147 : pari_CATCH(e_INV)
1742 : {
1743 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1744 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1745 0 : if (DEBUGLEVEL)
1746 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1747 0 : Rp = NULL;
1748 : } pari_TRY {
1749 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1750 147 : } pari_ENDCATCH
1751 147 : if (!Rp) continue;
1752 147 : if (typ(Rp)==t_COL)
1753 : {
1754 147 : Rp = FpC_sub(Rp,B,p);
1755 147 : if (ZV_equal0(Rp)) continue;
1756 : }
1757 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1758 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1759 147 : return gc_GEN(av, R);
1760 : }
1761 : }
1762 :
1763 : GEN
1764 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1765 : {
1766 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1767 : }
|
