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 25306746 : FpC_red(GEN x, GEN p)
28 218081437 : { 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 539 : FpV_red(GEN x, GEN p)
33 27742 : { pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
34 :
35 : GEN
36 6457441 : FpC_center(GEN x, GEN p, GEN pov2)
37 41660616 : { 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 4488319 : FpM_red(GEN x, GEN p)
46 19518091 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2586908 : FpM_center(GEN x, GEN p, GEN pov2)
50 7929326 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
51 :
52 : /* p != 3; assume entries in [0,p[ and ps2 = p>>1. */
53 : static void
54 77 : _FpC_center_inplace(GEN z, GEN p, GEN ps2)
55 : {
56 77 : long i, l = lg(z);
57 357 : for (i = 1; i < l; i++)
58 : { /* HACK: assume p != 3, which ensures u = gen_[0-2] is never written to */
59 280 : GEN u = gel(z,i);
60 280 : if (abscmpii(u, ps2) > 0) subiiz(u, p, u);
61 : }
62 77 : }
63 : static void
64 7 : _F3C_center_inplace(GEN z)
65 : {
66 7 : long i, l = lg(z);
67 35 : for (i = 1; i < l; i++) /* z[i] = 0, 1 : no-op */
68 28 : if (equaliu(gel(z,i), 2)) gel(z,i) = gen_m1;
69 7 : }
70 : void
71 0 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
72 : {
73 0 : if (equaliu(p, 3))
74 0 : _F3C_center_inplace(z);
75 : else
76 0 : _FpC_center_inplace(z, p, ps2);
77 0 : }
78 :
79 : void
80 42 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
81 : {
82 42 : long i, l = lg(z);
83 42 : if (equaliu(p, 3))
84 14 : for (i = 1; i < l; i++) _F3C_center_inplace(gel(z,i));
85 : else
86 112 : for (i = 1; i < l; i++) _FpC_center_inplace(gel(z,i), p, pov2);
87 42 : }
88 : GEN
89 11347 : Flm_center(GEN x, ulong p, ulong ps2)
90 175875 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
91 :
92 : GEN
93 119 : random_FpV(long d, GEN p)
94 : {
95 : long i;
96 119 : GEN y = cgetg(d+1,t_VEC);
97 25047 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
98 119 : return y;
99 : }
100 :
101 : GEN
102 413 : random_FpC(long d, GEN p)
103 : {
104 : long i;
105 413 : GEN y = cgetg(d+1,t_COL);
106 3409 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
107 413 : return y;
108 : }
109 :
110 : GEN
111 41229 : random_Flv(long n, ulong p)
112 : {
113 41229 : GEN y = cgetg(n+1, t_VECSMALL);
114 : long i;
115 260972 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
116 41230 : return y;
117 : }
118 :
119 : /********************************************************************/
120 : /** **/
121 : /** ADD, SUB **/
122 : /** **/
123 : /********************************************************************/
124 : GEN
125 242687 : FpC_add(GEN x, GEN y, GEN p)
126 3154231 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
127 :
128 : GEN
129 0 : FpV_add(GEN x, GEN y, GEN p)
130 0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
131 :
132 : GEN
133 0 : FpM_add(GEN x, GEN y, GEN p)
134 0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
135 :
136 : GEN
137 773116 : Flv_add(GEN x, GEN y, ulong p)
138 : {
139 773116 : if (p==2)
140 1349952 : pari_APPLY_ulong(x[i]^y[i])
141 : else
142 8233010 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
143 : }
144 :
145 : void
146 392929 : Flv_add_inplace(GEN x, GEN y, ulong p)
147 : {
148 392929 : long i, l = lg(x);
149 392929 : if (p==2)
150 1233910 : for (i = 1; i < l; i++) x[i] ^= y[i];
151 : else
152 164878 : for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
153 392929 : }
154 :
155 : ulong
156 3934 : Flv_sum(GEN x, ulong p)
157 : {
158 3934 : long i, l = lg(x);
159 3934 : ulong s = 0;
160 3934 : if (p==2)
161 17738 : for (i = 1; i < l; i++) s ^= x[i];
162 : else
163 0 : for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
164 3934 : return s;
165 : }
166 :
167 : GEN
168 493826 : FpC_sub(GEN x, GEN y, GEN p)
169 6679155 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
170 :
171 : GEN
172 0 : FpV_sub(GEN x, GEN y, GEN p)
173 0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
174 :
175 : GEN
176 0 : FpM_sub(GEN x, GEN y, GEN p)
177 0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
178 :
179 : GEN
180 210987909 : Flv_sub(GEN x, GEN y, ulong p)
181 957759377 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
182 :
183 : void
184 0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
185 : {
186 0 : long i, l = lg(x);
187 0 : for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
188 0 : }
189 :
190 : GEN
191 49653 : Flm_Fl_add(GEN x, ulong y, ulong p)
192 : {
193 49653 : long l = lg(x), i, j;
194 49653 : GEN z = cgetg(l,t_MAT);
195 :
196 49653 : if (l==1) return z;
197 49653 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
198 236092 : for (i=1; i<l; i++)
199 : {
200 186436 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
201 186437 : gel(z,i) = zi;
202 1306808 : for (j=1; j<l; j++) zi[j] = xi[j];
203 186437 : zi[i] = Fl_add(zi[i], y, p);
204 : }
205 49656 : return z;
206 : }
207 : GEN
208 3689 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
209 :
210 : GEN
211 9224 : Flm_add(GEN x, GEN y, ulong p)
212 158026 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
213 :
214 : GEN
215 25030890 : Flm_sub(GEN x, GEN y, ulong p)
216 236000450 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
217 :
218 : /********************************************************************/
219 : /** **/
220 : /** MULTIPLICATION **/
221 : /** **/
222 : /********************************************************************/
223 : GEN
224 1117848 : FpC_Fp_mul(GEN x, GEN y, GEN p)
225 12181709 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
226 :
227 : GEN
228 855599 : Flv_Fl_mul(GEN x, ulong y, ulong p)
229 11180917 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
230 :
231 : GEN
232 907 : Flv_Fl_div(GEN x, ulong y, ulong p)
233 : {
234 907 : return Flv_Fl_mul(x, Fl_inv(y, p), p);
235 : }
236 : void
237 0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
238 : {
239 0 : Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
240 0 : }
241 : GEN
242 754923 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
243 754923 : long i, j, h, l = lg(X);
244 754923 : GEN A = cgetg(l, t_MAT);
245 754923 : if (l == 1) return A;
246 754923 : h = lgcols(X);
247 3960113 : for (j=1; j<l; j++)
248 : {
249 3205365 : GEN a = cgetg(h, t_COL), x = gel(X, j);
250 23157885 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
251 3205190 : gel(A,j) = a;
252 : }
253 754748 : return A;
254 : }
255 :
256 : /* x *= y */
257 : void
258 7502833 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
259 : {
260 : long i;
261 7502833 : if (HIGHWORD(y | p))
262 8526788 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
263 : else
264 31726377 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
265 7502834 : }
266 : void
267 0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
268 0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
269 :
270 : /* set y *= x */
271 : void
272 0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
273 : {
274 0 : long i, j, m, l = lg(y);
275 0 : if (l == 1) return;
276 0 : m = lgcols(y);
277 0 : if (HIGHWORD(x | p))
278 0 : for(j=1; j<l; j++)
279 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
280 : else
281 0 : for(j=1; j<l; j++)
282 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
283 : }
284 :
285 : /* return x * y */
286 : static GEN
287 16911634 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
288 : {
289 16911634 : long i, j, m, l = lg(y);
290 16911634 : GEN z = cgetg(l, t_MAT);
291 16909256 : if (l == 1) return z;
292 16909256 : m = lgcols(y);
293 150621579 : for(j=1; j<l; j++) {
294 133623329 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
295 358651167 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
296 : }
297 16998250 : return z;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 8214717 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
303 : {
304 8214717 : long i, j, m, l = lg(y);
305 8214717 : GEN z = cgetg(l, t_MAT);
306 8214714 : if (l == 1) return z;
307 8214714 : m = lgcols(y);
308 44408764 : for(j=1; j<l; j++) {
309 36194061 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 116327343 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
311 : }
312 8214703 : return z;
313 : }
314 :
315 : /* return x * y */
316 : GEN
317 25090052 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
318 : {
319 25090052 : if (HIGHWORD(p))
320 16911437 : return Flm_Fl_mul_pre_i(y, x, p, pi);
321 : else
322 8178615 : return Flm_Fl_mul_OK(y, x, p);
323 : }
324 :
325 : /* return x * y */
326 : GEN
327 36107 : Flm_Fl_mul(GEN y, ulong x, ulong p)
328 : {
329 36107 : if (HIGHWORD(p))
330 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
331 : else
332 36033 : return Flm_Fl_mul_OK(y, x, p);
333 : }
334 :
335 : GEN
336 3983419 : Flv_neg(GEN x, ulong p)
337 30321279 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
338 :
339 : void
340 6188 : Flv_neg_inplace(GEN v, ulong p)
341 : {
342 : long i;
343 186959 : for (i = 1; i < lg(v); ++i)
344 180771 : v[i] = Fl_neg(v[i], p);
345 6188 : }
346 :
347 : GEN
348 314137 : Flm_neg(GEN x, ulong p)
349 4266975 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
350 :
351 : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
352 : static GEN
353 18669294 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
354 : {
355 18669294 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
356 : long k;
357 220689935 : for (k = 2; k < lx; k++)
358 : {
359 202022270 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
360 202017035 : if (signe(t)) c = addii(c, t);
361 : }
362 18667665 : return c;
363 : }
364 :
365 : static long
366 97688656 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
367 : {
368 : long k;
369 97688656 : long c = coeff(x,i,1) * y[1];
370 1498282053 : for (k = 2; k < lx; k++)
371 1400593397 : c += coeff(x,i,k) * y[k];
372 97688656 : return c;
373 : }
374 :
375 : GEN
376 6450502 : zm_zc_mul(GEN x, GEN y)
377 : {
378 6450502 : long lx = lg(x), l, i;
379 : GEN z;
380 6450502 : if (lx == 1) return cgetg(1, t_VECSMALL);
381 6450502 : l = lg(gel(x,1));
382 6450502 : z = cgetg(l,t_VECSMALL);
383 104139158 : for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
384 6450502 : return z;
385 : }
386 :
387 : GEN
388 2114 : zm_mul(GEN x, GEN y)
389 : {
390 2114 : long i,j,lx=lg(x), ly=lg(y);
391 : GEN z;
392 2114 : if (ly==1) return cgetg(1,t_MAT);
393 2114 : z = cgetg(ly,t_MAT);
394 2114 : if (lx==1)
395 : {
396 0 : for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
397 0 : return z;
398 : }
399 30849 : for (j=1; j<ly; j++)
400 28735 : gel(z,j) = zm_zc_mul(x, gel(y,j));
401 2114 : return z;
402 : }
403 :
404 : static ulong
405 634079489 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
406 : {
407 634079489 : ulong c = ucoeff(x,i,1) * uel(y,1);
408 : long k;
409 7665535329 : for (k = 2; k < lx; k++) {
410 7031455840 : c += ucoeff(x,i,k) * uel(y,k);
411 7031455840 : if (c & HIGHBIT) c %= p;
412 : }
413 634079489 : return c % p;
414 : }
415 :
416 : static ulong
417 537152034 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
418 : {
419 : ulong l0, l1, v1, h0, h1;
420 537152034 : long k = 1;
421 : LOCAL_OVERFLOW;
422 : LOCAL_HIREMAINDER;
423 537152034 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
424 9903292481 : while (++k < lx) {
425 9366140447 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
426 9366140447 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
427 : }
428 537152034 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
429 61380709 : else return remlll_pre(v1, h1, l1, p, pi);
430 : }
431 :
432 : static GEN
433 179834 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
434 : {
435 : long i,j;
436 179834 : GEN z = NULL;
437 :
438 1610932 : for (j=1; j<lx; j++)
439 : {
440 1431098 : if (!y[j]) continue;
441 443889 : if (!z) z = Flv_copy(gel(x,j));
442 12370868 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
443 : }
444 179834 : if (!z) z = zero_zv(l-1);
445 179834 : return z;
446 : }
447 :
448 : static GEN
449 1659228 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
450 : {
451 1659228 : GEN z = cgetg(l,t_COL);
452 : long i;
453 16028506 : for (i = 1; i < l; i++)
454 : {
455 14369317 : pari_sp av = avma;
456 14369317 : GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
457 14369072 : gel(z,i) = gerepileuptoint(av, modii(c,p));
458 : }
459 1659189 : return z;
460 : }
461 :
462 : static void
463 73414572 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
464 : {
465 : long i;
466 707531756 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
467 73429954 : }
468 : static GEN
469 70104469 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
470 : {
471 70104469 : GEN z = cgetg(l,t_VECSMALL);
472 70103648 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
473 70105427 : return z;
474 : }
475 :
476 : static void
477 93099415 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
478 : {
479 : long i;
480 632932234 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
481 93695189 : }
482 : static GEN
483 88995116 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
484 : {
485 88995116 : GEN z = cgetg(l,t_VECSMALL);
486 88863671 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
487 89047681 : return z;
488 : }
489 :
490 : GEN
491 1353907 : FpM_mul(GEN x, GEN y, GEN p)
492 : {
493 1353907 : long lx=lg(x), ly=lg(y);
494 : GEN z;
495 1353907 : pari_sp av = avma;
496 1353907 : if (ly==1) return cgetg(1,t_MAT);
497 1353907 : if (lx==1) return zeromat(0, ly-1);
498 1353907 : if (lgefint(p) == 3)
499 : {
500 1308242 : ulong pp = uel(p,2);
501 1308242 : if (pp == 2)
502 : {
503 674688 : x = ZM_to_F2m(x);
504 674706 : y = ZM_to_F2m(y);
505 674717 : z = F2m_to_ZM(F2m_mul(x,y));
506 : }
507 : else
508 : {
509 633554 : x = ZM_to_Flm(x, pp);
510 633575 : y = ZM_to_Flm(y, pp);
511 633578 : z = Flm_to_ZM(Flm_mul(x,y, pp));
512 : }
513 : } else
514 45665 : z = FpM_red(ZM_mul(x, y), p);
515 1353968 : return gerepileupto(av, z);
516 : }
517 :
518 : static GEN
519 27464649 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
520 : {
521 : long j;
522 27464649 : GEN z = cgetg(ly, t_MAT);
523 27463024 : if (SMALL_ULONG(p))
524 89420120 : for (j=1; j<ly; j++)
525 69822077 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
526 : else
527 96659042 : for (j=1; j<ly; j++)
528 88792366 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
529 27464719 : return z;
530 : }
531 :
532 : /*
533 : Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
534 : as an (m x n)-matrix, padding the input with zeroes as necessary.
535 : */
536 : static void
537 66663 : add_slices_ip(long m, long n,
538 : GEN A, long ma, long da, long na, long ea,
539 : GEN B, long mb, long db, long nb, long eb,
540 : GEN M, long dy, long dx, ulong p)
541 : {
542 66663 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
543 : GEN C;
544 :
545 2806979 : for (j = 1; j <= min_e; j++) {
546 2740316 : C = gel(M, j + dx) + dy;
547 125566102 : for (i = 1; i <= min_d; i++)
548 122825786 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
549 122825842 : ucoeff(B, mb + i, nb + j), p);
550 2807880 : for (; i <= da; i++)
551 67620 : uel(C, i) = ucoeff(A, ma + i, na + j);
552 2740260 : for (; i <= db; i++)
553 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
554 2740260 : for (; i <= m; i++)
555 0 : uel(C, i) = 0;
556 : }
557 70366 : for (; j <= ea; j++) {
558 3703 : C = gel(M, j + dx) + dy;
559 201369 : for (i = 1; i <= da; i++)
560 197666 : uel(C, i) = ucoeff(A, ma + i, na + j);
561 3703 : for (; i <= m; i++)
562 0 : uel(C, i) = 0;
563 : }
564 66663 : for (; j <= eb; j++) {
565 0 : C = gel(M, j + dx) + dy;
566 0 : for (i = 1; i <= db; i++)
567 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
568 0 : for (; i <= m; i++)
569 0 : uel(C, i) = 0;
570 : }
571 66663 : for (; j <= n; j++) {
572 0 : C = gel(M, j + dx) + dy;
573 0 : for (i = 1; i <= m; i++)
574 0 : uel(C, i) = 0;
575 : }
576 66663 : }
577 :
578 : static GEN
579 33332 : add_slices(long m, long n,
580 : GEN A, long ma, long da, long na, long ea,
581 : GEN B, long mb, long db, long nb, long eb, ulong p)
582 : {
583 : GEN M;
584 : long j;
585 33332 : M = cgetg(n + 1, t_MAT);
586 1428190 : for (j = 1; j <= n; j++)
587 1394860 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
588 33330 : add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
589 33331 : return M;
590 : }
591 :
592 : /*
593 : Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
594 : as an (m x n)-matrix, padding the input with zeroes as necessary.
595 : */
596 : static GEN
597 58330 : subtract_slices(long m, long n,
598 : GEN A, long ma, long da, long na, long ea,
599 : GEN B, long mb, long db, long nb, long eb, ulong p)
600 : {
601 58330 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
602 58331 : GEN M = cgetg(n + 1, t_MAT), C;
603 :
604 2549618 : for (j = 1; j <= min_e; j++) {
605 2491287 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
606 117613247 : for (i = 1; i <= min_d; i++)
607 115121984 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
608 115123186 : ucoeff(B, mb + i, nb + j), p);
609 2593569 : for (; i <= da; i++)
610 103508 : uel(C, i) = ucoeff(A, ma + i, na + j);
611 2568210 : for (; i <= db; i++)
612 78149 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
613 2490061 : for (; i <= m; i++)
614 0 : uel(C, i) = 0;
615 : }
616 58331 : for (; j <= ea; j++) {
617 0 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
618 0 : for (i = 1; i <= da; i++)
619 0 : uel(C, i) = ucoeff(A, ma + i, na + j);
620 0 : for (; i <= m; i++)
621 0 : uel(C, i) = 0;
622 : }
623 59795 : for (; j <= eb; j++) {
624 1464 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
625 85631 : for (i = 1; i <= db; i++)
626 84167 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
627 1464 : for (; i <= m; i++)
628 0 : uel(C, i) = 0;
629 : }
630 59795 : for (; j <= n; j++)
631 1464 : gel(M, j) = zero_Flv(m);
632 58331 : return M;
633 : }
634 :
635 : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
636 :
637 : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
638 : static GEN
639 8333 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
640 : {
641 : pari_sp av;
642 8333 : long m1 = (m + 1)/2, m2 = m/2,
643 8333 : n1 = (n + 1)/2, n2 = n/2,
644 8333 : p1 = (p + 1)/2, p2 = p/2;
645 : GEN A11, A12, A22, B11, B21, B22,
646 : S1, S2, S3, S4, T1, T2, T3, T4,
647 : M1, M2, M3, M4, M5, M6, M7,
648 : V1, V2, V3, C;
649 : long j;
650 8333 : C = cgetg(p + 1, t_MAT);
651 683152 : for (j = 1; j <= p; j++)
652 674820 : gel(C, j) = cgetg(m + 1, t_VECSMALL);
653 8332 : av = avma;
654 8332 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
655 8333 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
656 8333 : M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
657 8333 : if (gc_needed(av, 1))
658 0 : gerepileall(av, 2, &T2, &M2); /* destroy S1 */
659 8333 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
660 8333 : if (gc_needed(av, 1))
661 0 : gerepileall(av, 2, &M2, &T3); /* destroy T2 */
662 8333 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
663 8333 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
664 8333 : M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
665 8333 : if (gc_needed(av, 1))
666 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
667 8333 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
668 8333 : if (gc_needed(av, 1))
669 0 : gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
670 8333 : A11 = matslice(A, 1, m1, 1, n1);
671 8333 : B11 = matslice(B, 1, n1, 1, p1);
672 8333 : M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
673 8333 : if (gc_needed(av, 1))
674 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
675 8333 : A12 = matslice(A, 1, m1, n1 + 1, n);
676 8333 : B21 = matslice(B, n1 + 1, n, 1, p1);
677 8333 : M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
678 8333 : if (gc_needed(av, 1))
679 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
680 8333 : add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
681 8333 : if (gc_needed(av, 1))
682 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy M4 */
683 8333 : M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
684 8333 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
685 8333 : if (gc_needed(av, 1))
686 0 : gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4); /* destroy S3 */
687 8333 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
688 8333 : if (gc_needed(av, 1))
689 0 : gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4); /* destroy T3 */
690 8333 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
691 8333 : if (gc_needed(av, 1))
692 0 : gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1); /* destroy M1, M5 */
693 8333 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
694 8333 : M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
695 8333 : if (gc_needed(av, 1))
696 0 : gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6); /* destroy S4, B22 */
697 8333 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
698 8333 : M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
699 8333 : if (gc_needed(av, 1))
700 0 : gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7); /* destroy A22, T4 */
701 8333 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
702 8333 : add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
703 8333 : if (gc_needed(av, 1))
704 0 : gerepileall(av, 4, &M2, &M3, &V1, &M7); /* destroy V3, M6 */
705 8333 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
706 8332 : if (gc_needed(av, 1))
707 0 : gerepileall(av, 3, &M3, &M7, &V2); /* destroy V1, M2 */
708 8332 : add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
709 8333 : add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
710 8333 : set_avma(av); return C;
711 : }
712 :
713 : /* Strassen-Winograd used for dim >= ZM_sw_bound */
714 : static GEN
715 27471754 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
716 : {
717 27471754 : ulong e = expu(p);
718 : #ifdef LONG_IS_64BIT
719 23071554 : long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
720 : #else
721 4400350 : long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
722 : #endif
723 27471904 : if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
724 27463571 : return Flm_mul_classical(x, y, l, lx, ly, p, pi);
725 : else
726 8333 : return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
727 : }
728 :
729 : GEN
730 13379183 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
731 : {
732 13379183 : long lx=lg(x), ly=lg(y);
733 13379183 : if (ly==1) return cgetg(1,t_MAT);
734 13379183 : if (lx==1) return zero_Flm(0, ly-1);
735 13379183 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
736 : }
737 :
738 : GEN
739 14032020 : Flm_mul(GEN x, GEN y, ulong p)
740 : {
741 14032020 : long lx=lg(x), ly=lg(y);
742 14032020 : if (ly==1) return cgetg(1,t_MAT);
743 14031733 : if (lx==1) return zero_Flm(0, ly-1);
744 14031733 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
745 : }
746 :
747 : struct _Flm
748 : {
749 : ulong p;
750 : long n;
751 : };
752 :
753 : static GEN
754 8470 : _Flm_mul(void *E , GEN x, GEN y)
755 8470 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
756 : static GEN
757 43946 : _Flm_sqr(void *E, GEN x)
758 43946 : { return Flm_mul(x,x,((struct _Flm*)E)->p); }
759 : static GEN
760 861 : _Flm_one(void *E)
761 861 : { return matid_Flm(((struct _Flm*)E)->n); }
762 : GEN
763 24902 : Flm_powu(GEN x, ulong n, ulong p)
764 : {
765 : struct _Flm d;
766 24902 : if (!n) return matid(lg(x)-1);
767 24902 : d.p = p;
768 24902 : return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
769 : }
770 : GEN
771 861 : Flm_powers(GEN x, ulong n, ulong p)
772 : {
773 : struct _Flm d;
774 861 : d.p = p;
775 861 : d.n = lg(x)-1;
776 861 : return gen_powers(x, n, 1, (void*)&d,
777 : &_Flm_sqr, &_Flm_mul, &_Flm_one);
778 : }
779 :
780 : static GEN
781 0 : _FpM_mul(void *p , GEN x, GEN y)
782 0 : { return FpM_mul(x,y,(GEN)p); }
783 : static GEN
784 0 : _FpM_sqr(void *p, GEN x)
785 0 : { return FpM_mul(x,x,(GEN)p); }
786 : GEN
787 0 : FpM_powu(GEN x, ulong n, GEN p)
788 : {
789 0 : if (!n) return matid(lg(x)-1);
790 0 : if (lgefint(p) == 3)
791 : {
792 0 : pari_sp av = avma;
793 0 : ulong pp = uel(p,2);
794 : GEN z;
795 0 : if (pp == 2)
796 0 : z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
797 : else
798 0 : z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
799 0 : return gerepileupto(av, z);
800 : }
801 0 : return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
802 : }
803 :
804 : /*Multiple a column vector by a line vector to make a matrix*/
805 : GEN
806 14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
807 : {
808 14 : long i,j, lx=lg(x), ly=lg(y);
809 : GEN z;
810 14 : if (ly==1) return cgetg(1,t_MAT);
811 14 : z = cgetg(ly, t_MAT);
812 49 : for (j=1; j < ly; j++)
813 : {
814 35 : GEN zj = cgetg(lx,t_VECSMALL);
815 112 : for (i=1; i<lx; i++)
816 77 : uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
817 35 : gel(z,j) = zj;
818 : }
819 14 : return z;
820 : }
821 :
822 : /*Multiple a column vector by a line vector to make a matrix*/
823 : GEN
824 5670 : FpC_FpV_mul(GEN x, GEN y, GEN p)
825 : {
826 5670 : long i,j, lx=lg(x), ly=lg(y);
827 : GEN z;
828 5670 : if (ly==1) return cgetg(1,t_MAT);
829 5670 : z = cgetg(ly,t_MAT);
830 63861 : for (j=1; j < ly; j++)
831 : {
832 58191 : GEN zj = cgetg(lx,t_COL);
833 1379123 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
834 58191 : gel(z, j) = zj;
835 : }
836 5670 : return z;
837 : }
838 :
839 : /* Multiply a line vector by a column and return a scalar (t_INT) */
840 : GEN
841 4791974 : FpV_dotproduct(GEN x, GEN y, GEN p)
842 : {
843 4791974 : long i, lx = lg(x);
844 : pari_sp av;
845 : GEN c;
846 4791974 : if (lx == 1) return gen_0;
847 4791456 : av = avma; c = mulii(gel(x,1),gel(y,1));
848 16450081 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
849 4791363 : return gerepileuptoint(av, modii(c,p));
850 : }
851 : GEN
852 0 : FpV_dotsquare(GEN x, GEN p)
853 : {
854 0 : long i, lx = lg(x);
855 : pari_sp av;
856 : GEN c;
857 0 : if (lx == 1) return gen_0;
858 0 : av = avma; c = sqri(gel(x,1));
859 0 : for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
860 0 : return gerepileuptoint(av, modii(c,p));
861 : }
862 :
863 : INLINE ulong
864 8984906 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
865 : {
866 8984906 : ulong c = uel(x,0) * uel(y,0);
867 : long k;
868 64209493 : for (k = 1; k < lx; k++) {
869 55224587 : c += uel(x,k) * uel(y,k);
870 55224587 : if (c & HIGHBIT) c %= p;
871 : }
872 8984906 : return c % p;
873 : }
874 :
875 : INLINE ulong
876 739816 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
877 : {
878 : ulong l0, l1, v1, h0, h1;
879 739816 : long i = 0;
880 : LOCAL_OVERFLOW;
881 : LOCAL_HIREMAINDER;
882 739816 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
883 15409934 : while (++i < lx) {
884 14670118 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
885 14670118 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
886 : }
887 739816 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
888 464436 : else return remlll_pre(v1, h1, l1, p, pi);
889 : }
890 :
891 : ulong
892 356660 : Flv_dotproduct(GEN x, GEN y, ulong p)
893 : {
894 356660 : long lx = lg(x)-1;
895 356660 : if (lx == 0) return 0;
896 356660 : if (SMALL_ULONG(p))
897 356660 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
898 : else
899 0 : return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
900 : }
901 :
902 : ulong
903 58818 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
904 : {
905 58818 : long lx = lg(x)-1;
906 58818 : if (lx == 0) return 0;
907 58818 : if (SMALL_ULONG(p))
908 53260 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
909 : else
910 5558 : return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
911 : }
912 :
913 : ulong
914 9376130 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
915 : {
916 9376130 : long lx = minss(lgpol(x), lgpol(y));
917 9376238 : if (lx == 0) return 0;
918 9309249 : if (pi)
919 734258 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
920 : else
921 8574991 : return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
922 : }
923 : ulong
924 0 : Flx_dotproduct(GEN x, GEN y, ulong p)
925 0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
926 :
927 : GEN
928 1659228 : FpM_FpC_mul(GEN x, GEN y, GEN p)
929 : {
930 1659228 : long lx = lg(x);
931 1659228 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
932 : }
933 : GEN
934 699280 : Flm_Flc_mul(GEN x, GEN y, ulong p)
935 : {
936 699280 : long l, lx = lg(x);
937 699280 : if (lx==1) return cgetg(1,t_VECSMALL);
938 699280 : l = lgcols(x);
939 699279 : if (p==2)
940 179834 : return Flm_Flc_mul_i_2(x, y, lx, l);
941 519445 : else if (SMALL_ULONG(p))
942 282532 : return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
943 : else
944 236913 : return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
945 : }
946 :
947 : /* allow pi = 0 */
948 : GEN
949 6693 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
950 : {
951 6693 : long l, lx = lg(x);
952 : GEN z;
953 6693 : if (lx==1) return cgetg(1,t_VECSMALL);
954 6693 : l = lgcols(x);
955 6693 : z = cgetg(l, t_VECSMALL);
956 6693 : if (SMALL_ULONG(p))
957 4407 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
958 : else
959 2286 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
960 6693 : return z;
961 : }
962 :
963 : /* allow pi = 0 */
964 : GEN
965 7533691 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
966 : {
967 7533691 : long l, lx = lg(x);
968 : GEN z;
969 7533691 : if (lx==1) return pol0_Flx(sv);
970 7533691 : l = lgcols(x);
971 7530825 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
972 7533520 : if (SMALL_ULONG(p))
973 3306518 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
974 : else
975 4227002 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
976 7547280 : return Flx_renormalize(z, l + 1);
977 : }
978 :
979 : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
980 : GEN
981 1405600 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
982 : {
983 1405600 : long i, l, lx = lg(x);
984 : GEN z;
985 1405600 : if (lx==1) return pol_0(v);
986 1405600 : l = lgcols(x);
987 1405600 : z = new_chunk(l+1);
988 2137123 : for (i=l-1; i; i--)
989 : {
990 1954182 : pari_sp av = avma;
991 1954182 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
992 1954180 : p1 = modii(p1, p);
993 1954180 : if (signe(p1))
994 : {
995 1222657 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
996 1222657 : gel(z,i+1) = gerepileuptoint(av, p1);
997 1222659 : break;
998 : }
999 731523 : set_avma(av);
1000 : }
1001 1405600 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1002 1222659 : z[0] = evaltyp(t_POL) | evallg(i+2);
1003 1222659 : z[1] = evalsigne(1) | evalvarn(v);
1004 3568451 : for (; i; i--)
1005 : {
1006 2345793 : pari_sp av = avma;
1007 2345793 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1008 2345790 : gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
1009 : }
1010 1222658 : return z;
1011 : }
1012 :
1013 : /********************************************************************/
1014 : /** **/
1015 : /** TRANSPOSITION **/
1016 : /** **/
1017 : /********************************************************************/
1018 :
1019 : /* == zm_transpose */
1020 : GEN
1021 713637 : Flm_transpose(GEN x)
1022 : {
1023 713637 : long i, dx, lx = lg(x);
1024 : GEN y;
1025 713637 : if (lx == 1) return cgetg(1,t_MAT);
1026 713504 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1027 8472065 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1028 713505 : return y;
1029 : }
1030 :
1031 : /********************************************************************/
1032 : /** **/
1033 : /** SCALAR MATRICES **/
1034 : /** **/
1035 : /********************************************************************/
1036 :
1037 : GEN
1038 4022 : gen_matid(long n, void *E, const struct bb_field *S)
1039 : {
1040 4022 : GEN y = cgetg(n+1,t_MAT), _0, _1;
1041 : long i;
1042 4022 : if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
1043 4022 : _0 = S->s(E,0);
1044 4022 : _1 = S->s(E,1);
1045 17903 : for (i=1; i<=n; i++)
1046 : {
1047 13881 : GEN z = const_col(n, _0); gel(z,i) = _1;
1048 13881 : gel(y, i) = z;
1049 : }
1050 4022 : return y;
1051 : }
1052 :
1053 : GEN
1054 35 : matid_F2xqM(long n, GEN T)
1055 : {
1056 : void *E;
1057 35 : const struct bb_field *S = get_F2xq_field(&E, T);
1058 35 : return gen_matid(n, E, S);
1059 : }
1060 : GEN
1061 56 : matid_FlxqM(long n, GEN T, ulong p)
1062 : {
1063 : void *E;
1064 56 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1065 56 : return gen_matid(n, E, S);
1066 : }
1067 :
1068 : GEN
1069 1406355 : matid_Flm(long n)
1070 : {
1071 1406355 : GEN y = cgetg(n+1,t_MAT);
1072 : long i;
1073 1406348 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1074 8817738 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1075 1406358 : return y;
1076 : }
1077 :
1078 : GEN
1079 42 : scalar_Flm(long s, long n)
1080 : {
1081 : long i;
1082 42 : GEN y = cgetg(n+1,t_MAT);
1083 560 : for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
1084 42 : return y;
1085 : }
1086 :
1087 : /********************************************************************/
1088 : /** **/
1089 : /** CONVERSIONS **/
1090 : /** **/
1091 : /********************************************************************/
1092 : GEN
1093 48289985 : ZV_to_Flv(GEN x, ulong p)
1094 509364903 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1095 :
1096 : GEN
1097 11850511 : ZM_to_Flm(GEN x, ulong p)
1098 59226339 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1099 :
1100 : GEN
1101 2418 : ZMV_to_FlmV(GEN x, ulong m)
1102 18920 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1103 :
1104 : /* TO INTMOD */
1105 : static GEN
1106 20250476 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1107 : static GEN
1108 383844 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1109 :
1110 : GEN
1111 3495 : Fp_to_mod(GEN z, GEN p)
1112 : {
1113 3495 : retmkintmod(modii(z, p), icopy(p));
1114 : }
1115 :
1116 : GEN
1117 997241 : FpX_to_mod_raw(GEN z, GEN p)
1118 : {
1119 997241 : long i,l = lg(z);
1120 : GEN x;
1121 997241 : if (l == 2)
1122 : {
1123 522424 : x = cgetg(3,t_POL); x[1] = z[1];
1124 522424 : gel(x,2) = mkintmod(gen_0,p); return x;
1125 : }
1126 474817 : x = cgetg(l,t_POL);
1127 3918481 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1128 474817 : x[1] = z[1]; return normalizepol_lg(x,l);
1129 : }
1130 :
1131 : /* z in Z[X], return z * Mod(1,p), normalized*/
1132 : GEN
1133 1668749 : FpX_to_mod(GEN z, GEN p)
1134 : {
1135 1668749 : long i,l = lg(z);
1136 : GEN x;
1137 1668749 : if (l == 2)
1138 : {
1139 3157 : x = cgetg(3,t_POL); x[1] = z[1];
1140 3157 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1141 : }
1142 1665592 : x = cgetg(l,t_POL);
1143 1665592 : if (l >2) p = icopy(p);
1144 18369436 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1145 1665592 : x[1] = z[1]; return normalizepol_lg(x,l);
1146 : }
1147 :
1148 : GEN
1149 84021 : FpXC_to_mod(GEN z, GEN p)
1150 : {
1151 84021 : long i,l = lg(z);
1152 84021 : GEN x = cgetg(l, t_COL);
1153 84021 : p = icopy(p);
1154 474166 : for (i=1; i<l; i++)
1155 390145 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1156 84021 : return x;
1157 : }
1158 :
1159 : static GEN
1160 0 : FpXC_to_mod_raw(GEN x, GEN p)
1161 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1162 :
1163 : GEN
1164 0 : FpXM_to_mod(GEN z, GEN p)
1165 : {
1166 0 : long i,l = lg(z);
1167 0 : GEN x = cgetg(l, t_MAT);
1168 0 : p = icopy(p);
1169 0 : for (i=1; i<l; i++)
1170 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1171 0 : return x;
1172 : }
1173 :
1174 : /* z in Z^n, return z * Mod(1,p), normalized*/
1175 : GEN
1176 35963 : FpV_to_mod(GEN z, GEN p)
1177 : {
1178 35963 : long i,l = lg(z);
1179 35963 : GEN x = cgetg(l, t_VEC);
1180 35963 : if (l == 1) return x;
1181 35963 : p = icopy(p);
1182 72325 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1183 35963 : return x;
1184 : }
1185 : /* z in Z^n, return z * Mod(1,p), normalized*/
1186 : GEN
1187 105 : FpC_to_mod(GEN z, GEN p)
1188 : {
1189 105 : long i, l = lg(z);
1190 105 : GEN x = cgetg(l, t_COL);
1191 105 : if (l == 1) return x;
1192 91 : p = icopy(p);
1193 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1194 91 : return x;
1195 : }
1196 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1197 : GEN
1198 226 : FpM_to_mod(GEN z, GEN p)
1199 : {
1200 226 : long i, j, m, l = lg(z);
1201 226 : GEN x = cgetg(l,t_MAT), y, zi;
1202 226 : if (l == 1) return x;
1203 205 : m = lgcols(z);
1204 205 : p = icopy(p);
1205 2456 : for (i=1; i<l; i++)
1206 : {
1207 2251 : gel(x,i) = cgetg(m,t_COL);
1208 2251 : y = gel(x,i); zi= gel(z,i);
1209 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1210 : }
1211 205 : return x;
1212 : }
1213 : GEN
1214 28 : Flc_to_mod(GEN z, ulong pp)
1215 : {
1216 28 : long i, l = lg(z);
1217 28 : GEN p, x = cgetg(l, t_COL);
1218 28 : if (l == 1) return x;
1219 28 : p = utoipos(pp);
1220 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1221 28 : return x;
1222 : }
1223 : GEN
1224 56292 : Flm_to_mod(GEN z, ulong pp)
1225 : {
1226 56292 : long i, j, m, l = lg(z);
1227 56292 : GEN p, x = cgetg(l,t_MAT), y, zi;
1228 56292 : if (l == 1) return x;
1229 56264 : m = lgcols(z);
1230 56264 : p = utoipos(pp);
1231 179348 : for (i=1; i<l; i++)
1232 : {
1233 123084 : gel(x,i) = cgetg(m,t_COL);
1234 123084 : y = gel(x,i); zi= gel(z,i);
1235 506844 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1236 : }
1237 56264 : return x;
1238 : }
1239 :
1240 : GEN
1241 574 : FpVV_to_mod(GEN z, GEN p)
1242 : {
1243 574 : long i, j, m, l = lg(z);
1244 574 : GEN x = cgetg(l,t_VEC), y, zi;
1245 574 : if (l == 1) return x;
1246 574 : m = lgcols(z);
1247 574 : p = icopy(p);
1248 1246 : for (i=1; i<l; i++)
1249 : {
1250 672 : gel(x,i) = cgetg(m,t_VEC);
1251 672 : y = gel(x,i); zi= gel(z,i);
1252 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1253 : }
1254 574 : return x;
1255 : }
1256 :
1257 : /* z in Z^n, return z * Mod(1,p), normalized*/
1258 : GEN
1259 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1260 : {
1261 7 : long i,l = lg(z);
1262 7 : GEN x = cgetg(l, t_COL);
1263 7 : if (l == 1) return x;
1264 7 : p = icopy(p);
1265 7 : T = FpX_to_mod_raw(T, p);
1266 21 : for (i=1; i<l; i++)
1267 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1268 7 : return x;
1269 : }
1270 :
1271 : static GEN
1272 582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1273 : {
1274 582533 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1275 582533 : return mkpolmod(z, T);
1276 : }
1277 :
1278 : /* z in Z^n, return z * Mod(1,p), normalized*/
1279 : GEN
1280 28 : FqC_to_mod(GEN z, GEN T, GEN p)
1281 : {
1282 : GEN x;
1283 28 : long i,l = lg(z);
1284 28 : if (!T) return FpC_to_mod(z, p);
1285 28 : x = cgetg(l, t_COL);
1286 28 : if (l == 1) return x;
1287 28 : p = icopy(p);
1288 28 : T = FpX_to_mod_raw(T, p);
1289 504 : for (i=1; i<l; i++)
1290 476 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1291 28 : return x;
1292 : }
1293 :
1294 : /* z in Z^n, return z * Mod(1,p), normalized*/
1295 : static GEN
1296 107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1297 690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1298 :
1299 : /* z in Z^n, return z * Mod(1,p), normalized*/
1300 : GEN
1301 26145 : FqM_to_mod(GEN z, GEN T, GEN p)
1302 : {
1303 : GEN x;
1304 26145 : long i,l = lg(z);
1305 26145 : if (!T) return FpM_to_mod(z, p);
1306 26145 : x = cgetg(l, t_MAT);
1307 26145 : if (l == 1) return x;
1308 26117 : p = icopy(p);
1309 26117 : T = FpX_to_mod_raw(T, p);
1310 134092 : for (i=1; i<l; i++)
1311 107975 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1312 26117 : return x;
1313 : }
1314 :
1315 : /********************************************************************/
1316 : /* */
1317 : /* Blackbox linear algebra */
1318 : /* */
1319 : /********************************************************************/
1320 :
1321 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1322 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1323 : * (e_j) is the canonical basis.
1324 : * A sparse matrix (zMs) is a t_VEC of zCs */
1325 :
1326 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1327 : * integer representing an element of Fp. This is important since p can be
1328 : * large and p+E[i] would not fit in a C long. */
1329 :
1330 : /* RgCs and RgMs are similar, except that the type of the component is
1331 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1332 : * of E. */
1333 :
1334 : /* Most functions take an argument nbrow which is the number of lines of the
1335 : * column/matrix, which cannot be derived from the data. */
1336 :
1337 : GEN
1338 0 : zCs_to_ZC(GEN R, long nbrow)
1339 : {
1340 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1341 0 : long j, l = lg(C);
1342 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1343 0 : return c;
1344 : }
1345 :
1346 : GEN
1347 0 : zMs_to_ZM(GEN x, long nbrow)
1348 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1349 :
1350 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1351 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1352 : GEN
1353 119 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1354 : {
1355 119 : pari_sp ltop = avma;
1356 119 : long col = 0, n = lg(B)-1, m = 2*n+1;
1357 119 : if (ZV_equal0(B)) return zerocol(n);
1358 119 : while (++col <= n)
1359 : {
1360 119 : pari_sp btop = avma, av;
1361 : long i, lQ;
1362 119 : GEN V, Q, M, W = B;
1363 119 : GEN b = cgetg(m+2, t_POL);
1364 119 : b[1] = evalsigne(1)|evalvarn(0);
1365 119 : gel(b, 2) = gel(W, col);
1366 49975 : for (i = 3; i<m+2; i++)
1367 49856 : gel(b, i) = cgeti(lgefint(p));
1368 119 : av = avma;
1369 49975 : for (i = 3; i<m+2; i++)
1370 : {
1371 49856 : W = f(E, W);
1372 49856 : affii(gel(W, col),gel(b, i));
1373 49856 : if (gc_needed(av,1))
1374 : {
1375 1084 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1376 1084 : W = gerepileupto(av, W);
1377 : }
1378 : }
1379 119 : b = FpX_renormalize(b, m+2);
1380 119 : if (lgpol(b)==0) {set_avma(btop); continue; }
1381 119 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1382 119 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1383 119 : W = B; lQ =lg(Q);
1384 119 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1385 119 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1386 119 : av = avma;
1387 24599 : for (i = lQ-3; i > 1; i--)
1388 : {
1389 24480 : W = f(E, W);
1390 24480 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1391 24480 : if (gc_needed(av,1))
1392 : {
1393 864 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1394 864 : gerepileall(av, 2, &V, &W);
1395 : }
1396 : }
1397 119 : V = FpC_red(V, p);
1398 119 : W = FpC_sub(f(E,V), B, p);
1399 119 : if (ZV_equal0(W)) return gerepilecopy(ltop, V);
1400 0 : av = avma;
1401 0 : for (i = 1; i <= n; ++i)
1402 : {
1403 0 : V = W;
1404 0 : W = f(E, V);
1405 0 : if (ZV_equal0(W))
1406 0 : return gerepilecopy(ltop, shallowtrans(V));
1407 0 : gerepileall(av, 2, &V, &W);
1408 : }
1409 0 : set_avma(btop);
1410 : }
1411 0 : return NULL;
1412 : }
1413 :
1414 : GEN
1415 0 : zMs_ZC_mul(GEN M, GEN B)
1416 : {
1417 : long i, j;
1418 0 : long n = lg(B)-1;
1419 0 : GEN V = zerocol(n);
1420 0 : for (i = 1; i <= n; ++i)
1421 0 : if (signe(gel(B, i)))
1422 : {
1423 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1424 0 : long l = lg(C);
1425 0 : for (j = 1; j < l; ++j)
1426 : {
1427 0 : long k = C[j];
1428 0 : switch(E[j])
1429 : {
1430 0 : case 1:
1431 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1432 0 : break;
1433 0 : case -1:
1434 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1435 0 : break;
1436 0 : default:
1437 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]));
1438 0 : break;
1439 : }
1440 : }
1441 : }
1442 0 : return V;
1443 : }
1444 :
1445 : GEN
1446 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1447 :
1448 : GEN
1449 74574 : ZV_zMs_mul(GEN B, GEN M)
1450 : {
1451 : long i, j;
1452 74574 : long m = lg(M)-1;
1453 74574 : GEN V = cgetg(m+1,t_VEC);
1454 38973800 : for (i = 1; i <= m; ++i)
1455 : {
1456 38899226 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1457 38899226 : long l = lg(C);
1458 : GEN z;
1459 38899226 : if (l == 1)
1460 : {
1461 0 : gel(V,i) = gen_0;
1462 0 : continue;
1463 : }
1464 38899226 : z = mulis(gel(B, C[1]), E[1]);
1465 327097845 : for (j = 2; j < l; ++j)
1466 : {
1467 288198619 : long k = C[j];
1468 288198619 : switch(E[j])
1469 : {
1470 194218066 : case 1:
1471 194218066 : z = addii(z, gel(B,k));
1472 194218066 : break;
1473 78587594 : case -1:
1474 78587594 : z = subii(z, gel(B,k));
1475 78587594 : break;
1476 15392959 : default:
1477 15392959 : z = addii(z, mulis(gel(B,k), E[j]));
1478 15392959 : break;
1479 : }
1480 : }
1481 38899226 : gel(V,i) = z;
1482 : }
1483 74574 : return V;
1484 : }
1485 :
1486 : GEN
1487 119 : FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
1488 :
1489 : GEN
1490 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1491 : {
1492 0 : pari_sp av = avma, av2;
1493 0 : GEN xi, xb, pi = gen_1, P;
1494 : long i;
1495 0 : if (!C) {
1496 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1497 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1498 : }
1499 0 : P = utoipos(p);
1500 0 : av2 = avma;
1501 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1502 0 : xb = Flm_to_ZM(xi);
1503 0 : for (i = 2; i <= e; i++)
1504 : {
1505 0 : pi = muliu(pi, p); /* = p^(i-1) */
1506 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1507 0 : if (gc_needed(av,2))
1508 : {
1509 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1510 0 : gerepileall(av2,3, &pi,&b,&xb);
1511 : }
1512 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1513 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1514 : }
1515 0 : return gerepileupto(av, xb);
1516 : }
1517 :
1518 : struct wrapper_modp_s {
1519 : GEN (*f)(void*, GEN);
1520 : void *E;
1521 : GEN p;
1522 : };
1523 :
1524 : /* compute f(x) mod p */
1525 : static GEN
1526 74455 : wrap_relcomb_modp(void *E, GEN x)
1527 : {
1528 74455 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1529 74455 : return FpC_red(W->f(W->E, x), W->p);
1530 : }
1531 : static GEN
1532 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1533 :
1534 : static GEN
1535 74455 : wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
1536 :
1537 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1538 : GEN
1539 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1540 : {
1541 : struct wrapper_modp_s W;
1542 0 : pari_sp av = avma;
1543 0 : GEN xb, xi, pi = gen_1;
1544 : long i;
1545 0 : W.E = E;
1546 0 : W.f = f;
1547 0 : W.p = p;
1548 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1549 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1550 0 : xb = xi;
1551 0 : for (i = 2; i <= e; i++)
1552 : {
1553 0 : pi = mulii(pi, p); /* = p^(i-1) */
1554 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1555 0 : if (gc_needed(av,2))
1556 : {
1557 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1558 0 : gerepileall(av,3, &pi,&B,&xb);
1559 : }
1560 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1561 0 : if (!xi) return NULL;
1562 0 : if (typ(xi) == t_VEC) return gerepileupto(av, xi);
1563 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1564 : }
1565 0 : return gerepileupto(av, xb);
1566 : }
1567 :
1568 : static GEN
1569 24928 : vecprow(GEN A, GEN prow)
1570 : {
1571 24928 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1572 : }
1573 :
1574 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1575 : * or the index of a column which is linearly dependent from the others as a
1576 : * t_VECSMALL with a single component. */
1577 : GEN
1578 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1579 : {
1580 0 : pari_sp av = avma;
1581 : GEN pcol, prow;
1582 0 : long nbi=lg(M)-1, lR;
1583 : long i, n;
1584 : GEN Mp, Ap, Rp;
1585 : pari_timer ti;
1586 0 : if (DEBUGLEVEL) timer_start(&ti);
1587 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1588 0 : if (!pcol) return gc_NULL(av);
1589 0 : if (DEBUGLEVEL)
1590 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1591 0 : n = lg(pcol)-1;
1592 0 : Mp = cgetg(n+1, t_MAT);
1593 0 : for(i=1; i<=n; i++)
1594 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1595 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1596 0 : if (DEBUGLEVEL) timer_start(&ti);
1597 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1598 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1599 0 : if (!Rp) return gc_NULL(av);
1600 0 : lR = lg(Rp)-1;
1601 0 : if (typ(Rp) == t_COL)
1602 : {
1603 0 : GEN R = zerocol(nbi+1);
1604 0 : for (i=1; i<=lR; i++)
1605 0 : gel(R,pcol[i]) = gel(Rp,i);
1606 0 : return gerepilecopy(av, R);
1607 : }
1608 0 : for (i = 1; i <= lR; ++i)
1609 0 : if (signe(gel(Rp, i)))
1610 0 : return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
1611 : return NULL; /* LCOV_EXCL_LINE */
1612 : }
1613 :
1614 : GEN
1615 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1616 : {
1617 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1618 : }
1619 :
1620 : GEN
1621 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1622 : {
1623 : GEN res;
1624 0 : pari_CATCH(e_INV)
1625 : {
1626 0 : GEN E = pari_err_last();
1627 0 : GEN x = err_get_compo(E,2);
1628 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1629 0 : if (DEBUGLEVEL)
1630 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1631 0 : res = NULL;
1632 : } pari_TRY {
1633 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1634 0 : } pari_ENDCATCH
1635 0 : return res;
1636 : }
1637 :
1638 : static GEN
1639 119 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1640 : {
1641 119 : long i, j, oldf = 0, f = 0;
1642 119 : long lrow = lg(prow), lM = lg(M);
1643 119 : GEN W = const_vecsmall(lM-1,1);
1644 119 : GEN R = cgetg(lrow, t_VEC);
1645 235298 : for (i=1; i<lrow; i++)
1646 235179 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1647 : do
1648 : {
1649 393 : oldf = f;
1650 436630 : for(i=1; i<lM; i++)
1651 436237 : if (W[i])
1652 : {
1653 147868 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1654 147868 : long c=0, cj=0, lF = lg(F);
1655 1294915 : for(j=1; j<lF; j++)
1656 1147047 : if (!gel(R,F[j])) { c++; cj=j; }
1657 147868 : if (c>=2) continue;
1658 127155 : if (c==1)
1659 : {
1660 41413 : pari_sp av = avma;
1661 41413 : GEN s = gen_0;
1662 359128 : for(j=1; j<lF; j++)
1663 317715 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1664 : /* s /= -E[cj] mod p */
1665 41413 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1666 41413 : gel(R,F[cj]) = gerepileuptoint(av, s);
1667 41413 : f++;
1668 : }
1669 127155 : W[i]=0;
1670 : }
1671 393 : } while(oldf!=f);
1672 235298 : for (i=1; i<lrow; i++)
1673 235179 : if (!gel(R,i)) gel(R,i) = gen_0;
1674 119 : return R;
1675 : }
1676 :
1677 : /* Return a linear form Y such that YM is essentially 0 */
1678 : GEN
1679 119 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1680 : {
1681 119 : pari_sp av = avma, av2;
1682 : GEN pcol, prow;
1683 : long i, n;
1684 : GEN Mp, B, MB, R, Rp;
1685 : pari_timer ti;
1686 : struct wrapper_modp_s W;
1687 119 : if (DEBUGLEVEL) timer_start(&ti);
1688 119 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1689 119 : if (!pcol) return gc_NULL(av);
1690 119 : if (DEBUGLEVEL)
1691 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1692 119 : n = lg(pcol)-1;
1693 119 : Mp = cgetg(n+1, t_MAT);
1694 25047 : for(i=1; i<=n; i++)
1695 24928 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1696 119 : W.E = (void*) Mp;
1697 119 : W.f = wrap_relker;
1698 119 : W.p = p;
1699 119 : av2 = avma;
1700 : for(;;)
1701 : {
1702 119 : set_avma(av2);
1703 119 : B = random_FpV(n, p);
1704 119 : MB = FpV_FpMs_mul(B, Mp, p);
1705 119 : if (DEBUGLEVEL) timer_start(&ti);
1706 119 : pari_CATCH(e_INV)
1707 : {
1708 0 : GEN E = pari_err_last();
1709 0 : GEN x = err_get_compo(E,2);
1710 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1711 0 : if (DEBUGLEVEL)
1712 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1713 0 : Rp = NULL;
1714 : } pari_TRY {
1715 119 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
1716 119 : } pari_ENDCATCH
1717 119 : if (!Rp) continue;
1718 119 : if (typ(Rp)==t_COL)
1719 : {
1720 119 : Rp = FpC_sub(Rp,B,p);
1721 119 : if (ZV_equal0(Rp)) continue;
1722 : }
1723 119 : R = FpMs_structelim_back(M, Rp, prow, p);
1724 119 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1725 119 : return gerepilecopy(av, R);
1726 : }
1727 : }
1728 :
1729 : GEN
1730 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1731 : {
1732 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1733 : }
|
