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 27267548 : FpC_red(GEN x, GEN p)
28 187595286 : { 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 6802184 : FpC_center(GEN x, GEN p, GEN pov2)
37 42369001 : { 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 5052878 : FpM_red(GEN x, GEN p)
46 21777491 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2740765 : FpM_center(GEN x, GEN p, GEN pov2)
50 8226908 : { 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 147 : random_FpV(long d, GEN p)
94 : {
95 : long i;
96 147 : GEN y = cgetg(d+1,t_VEC);
97 23186 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
98 147 : 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 41289 : random_Flv(long n, ulong p)
112 : {
113 41289 : GEN y = cgetg(n+1, t_VECSMALL);
114 : long i;
115 261217 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
116 41289 : return y;
117 : }
118 :
119 : /********************************************************************/
120 : /** **/
121 : /** ADD, SUB **/
122 : /** **/
123 : /********************************************************************/
124 : GEN
125 242716 : FpC_add(GEN x, GEN y, GEN p)
126 3154376 : { 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 775759 : Flv_add(GEN x, GEN y, ulong p)
138 : {
139 775759 : if (p==2)
140 1350438 : pari_APPLY_ulong(x[i]^y[i])
141 : else
142 8278219 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
143 : }
144 :
145 : void
146 436243 : Flv_add_inplace(GEN x, GEN y, ulong p)
147 : {
148 436243 : long i, l = lg(x);
149 436243 : if (p==2)
150 1348322 : 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 436243 : }
154 :
155 : ulong
156 3878 : Flv_sum(GEN x, ulong p)
157 : {
158 3878 : long i, l = lg(x);
159 3878 : ulong s = 0;
160 3878 : if (p==2)
161 17710 : 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 3878 : return s;
165 : }
166 :
167 : GEN
168 493882 : FpC_sub(GEN x, GEN y, GEN p)
169 6675433 : { 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 210519254 : Flv_sub(GEN x, GEN y, ulong p)
181 953984056 : { 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 49341 : Flm_Fl_add(GEN x, ulong y, ulong p)
192 : {
193 49341 : long l = lg(x), i, j;
194 49341 : GEN z = cgetg(l,t_MAT);
195 :
196 49341 : if (l==1) return z;
197 49341 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
198 235268 : for (i=1; i<l; i++)
199 : {
200 185927 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
201 185927 : gel(z,i) = zi;
202 1306076 : for (j=1; j<l; j++) zi[j] = xi[j];
203 185927 : zi[i] = Fl_add(zi[i], y, p);
204 : }
205 49341 : 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 9228 : Flm_add(GEN x, GEN y, ulong p)
212 158062 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
213 :
214 : GEN
215 24899858 : Flm_sub(GEN x, GEN y, ulong p)
216 235403614 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
217 :
218 : /********************************************************************/
219 : /** **/
220 : /** MULTIPLICATION **/
221 : /** **/
222 : /********************************************************************/
223 : GEN
224 882819 : FpC_Fp_mul(GEN x, GEN y, GEN p)
225 11093890 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
226 :
227 : GEN
228 858974 : Flv_Fl_mul(GEN x, ulong y, ulong p)
229 11247974 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
230 :
231 : GEN
232 4697 : Flv_Fl_div(GEN x, ulong y, ulong p)
233 : {
234 4697 : 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 754609 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
243 754609 : long i, j, h, l = lg(X);
244 754609 : GEN A = cgetg(l, t_MAT);
245 754608 : if (l == 1) return A;
246 754608 : h = lgcols(X);
247 3987659 : for (j=1; j<l; j++)
248 : {
249 3233109 : GEN a = cgetg(h, t_COL), x = gel(X, j);
250 23226049 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
251 3233051 : gel(A,j) = a;
252 : }
253 754550 : return A;
254 : }
255 :
256 : /* x *= y */
257 : void
258 7724327 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
259 : {
260 : long i;
261 7724327 : if (HIGHWORD(y | p))
262 8555822 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
263 : else
264 31957781 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
265 7724325 : }
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 16760287 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
288 : {
289 16760287 : long i, j, m, l = lg(y);
290 16760287 : GEN z = cgetg(l, t_MAT);
291 16758716 : if (l == 1) return z;
292 16758716 : m = lgcols(y);
293 150038453 : for(j=1; j<l; j++) {
294 133166976 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
295 357560400 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
296 : }
297 16871477 : return z;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 8138890 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
303 : {
304 8138890 : long i, j, m, l = lg(y);
305 8138890 : GEN z = cgetg(l, t_MAT);
306 8138884 : if (l == 1) return z;
307 8138884 : m = lgcols(y);
308 44281551 : for(j=1; j<l; j++) {
309 36142670 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 115289126 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
311 : }
312 8138881 : return z;
313 : }
314 :
315 : /* return x * y */
316 : GEN
317 24863158 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
318 : {
319 24863158 : if (HIGHWORD(p))
320 16760349 : return Flm_Fl_mul_pre_i(y, x, p, pi);
321 : else
322 8102809 : return Flm_Fl_mul_OK(y, x, p);
323 : }
324 :
325 : /* return x * y */
326 : GEN
327 36111 : Flm_Fl_mul(GEN y, ulong x, ulong p)
328 : {
329 36111 : if (HIGHWORD(p))
330 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
331 : else
332 36037 : return Flm_Fl_mul_OK(y, x, p);
333 : }
334 :
335 : GEN
336 4002282 : Flv_neg(GEN x, ulong p)
337 30389489 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
338 :
339 : void
340 6202 : Flv_neg_inplace(GEN v, ulong p)
341 : {
342 : long i;
343 187994 : for (i = 1; i < lg(v); ++i)
344 181792 : v[i] = Fl_neg(v[i], p);
345 6202 : }
346 :
347 : GEN
348 314725 : Flm_neg(GEN x, ulong p)
349 4286470 : { 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 18697724 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
354 : {
355 18697724 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
356 : long k;
357 220737623 : for (k = 2; k < lx; k++)
358 : {
359 202041306 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
360 202037257 : if (signe(t)) c = addii(c, t);
361 : }
362 18696317 : 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 631913917 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
406 : {
407 631913917 : ulong c = ucoeff(x,i,1) * uel(y,1);
408 : long k;
409 7656149047 : for (k = 2; k < lx; k++) {
410 7024235130 : c += ucoeff(x,i,k) * uel(y,k);
411 7024235130 : if (c & HIGHBIT) c %= p;
412 : }
413 631913917 : return c % p;
414 : }
415 :
416 : static ulong
417 525833592 : 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 525833592 : long k = 1;
421 : LOCAL_OVERFLOW;
422 : LOCAL_HIREMAINDER;
423 525833592 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
424 8900361891 : while (++k < lx) {
425 8374528299 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
426 8374528299 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
427 : }
428 525833592 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
429 55477956 : else return remlll_pre(v1, h1, l1, p, pi);
430 : }
431 :
432 : static GEN
433 179884 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
434 : {
435 : long i,j;
436 179884 : GEN z = NULL;
437 :
438 1611132 : for (j=1; j<lx; j++)
439 : {
440 1431248 : if (!y[j]) continue;
441 443125 : if (!z) z = Flv_copy(gel(x,j));
442 12319412 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
443 : }
444 179884 : if (!z) z = zero_zv(l-1);
445 179884 : return z;
446 : }
447 :
448 : static GEN
449 1657401 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
450 : {
451 1657401 : GEN z = cgetg(l,t_COL);
452 : long i;
453 16016099 : for (i = 1; i < l; i++)
454 : {
455 14358742 : pari_sp av = avma;
456 14358742 : GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
457 14358503 : gel(z,i) = gerepileuptoint(av, modii(c,p));
458 : }
459 1657357 : return z;
460 : }
461 :
462 : static void
463 73278065 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
464 : {
465 : long i;
466 705194569 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
467 73282596 : }
468 : static GEN
469 69966934 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
470 : {
471 69966934 : GEN z = cgetg(l,t_VECSMALL);
472 69965866 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
473 69967245 : return z;
474 : }
475 :
476 : static void
477 92719841 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
478 : {
479 : long i;
480 620447397 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
481 92843526 : }
482 : static GEN
483 88578177 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
484 : {
485 88578177 : GEN z = cgetg(l,t_VECSMALL);
486 88433983 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
487 88631976 : return z;
488 : }
489 :
490 : GEN
491 1350694 : FpM_mul(GEN x, GEN y, GEN p)
492 : {
493 1350694 : long lx=lg(x), ly=lg(y);
494 : GEN z;
495 1350694 : pari_sp av = avma;
496 1350694 : if (ly==1) return cgetg(1,t_MAT);
497 1350694 : if (lx==1) return zeromat(0, ly-1);
498 1350694 : if (lgefint(p) == 3)
499 : {
500 1305012 : ulong pp = uel(p,2);
501 1305012 : if (pp == 2)
502 : {
503 671969 : x = ZM_to_F2m(x);
504 672025 : y = ZM_to_F2m(y);
505 672029 : z = F2m_to_ZM(F2m_mul(x,y));
506 : }
507 : else
508 : {
509 633043 : x = ZM_to_Flm(x, pp);
510 633045 : y = ZM_to_Flm(y, pp);
511 633047 : z = Flm_to_ZM(Flm_mul(x,y, pp));
512 : }
513 : } else
514 45682 : z = FpM_red(ZM_mul(x, y), p);
515 1350755 : return gerepileupto(av, z);
516 : }
517 :
518 : static GEN
519 27419606 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
520 : {
521 : long j;
522 27419606 : GEN z = cgetg(ly, t_MAT);
523 27419226 : if (SMALL_ULONG(p))
524 89290183 : for (j=1; j<ly; j++)
525 69681945 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
526 : else
527 96186099 : for (j=1; j<ly; j++)
528 88373811 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
529 27420526 : 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 66664 : 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 66664 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
543 : GEN C;
544 :
545 2807310 : for (j = 1; j <= min_e; j++) {
546 2740646 : C = gel(M, j + dx) + dy;
547 125641603 : for (i = 1; i <= min_d; i++)
548 122900957 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
549 122901706 : ucoeff(B, mb + i, nb + j), p);
550 2807517 : for (; i <= da; i++)
551 67620 : uel(C, i) = ucoeff(A, ma + i, na + j);
552 2739897 : for (; i <= db; i++)
553 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
554 2739897 : for (; i <= m; i++)
555 0 : uel(C, i) = 0;
556 : }
557 70367 : 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 66664 : 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 66664 : 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 66664 : }
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 1428206 : for (j = 1; j <= n; j++)
587 1394873 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
588 33333 : add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
589 33332 : 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 58330 : GEN M = cgetg(n + 1, t_MAT), C;
603 :
604 2549591 : for (j = 1; j <= min_e; j++) {
605 2491260 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
606 117631675 : for (i = 1; i <= min_d; i++)
607 115140427 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
608 115140289 : ucoeff(B, mb + i, nb + j), p);
609 2594894 : for (; i <= da; i++)
610 103508 : uel(C, i) = ucoeff(A, ma + i, na + j);
611 2569535 : for (; i <= db; i++)
612 78149 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
613 2491386 : 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 683181 : for (j = 1; j <= p; j++)
652 674849 : 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 8333 : if (gc_needed(av, 1))
707 0 : gerepileall(av, 3, &M3, &M7, &V2); /* destroy V1, M2 */
708 8333 : 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 27427243 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
716 : {
717 27427243 : ulong e = expu(p);
718 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
719 23039390 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
720 : #else
721 4388015 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
722 : #endif
723 27427405 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
724 27419072 : 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 13309015 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
731 : {
732 13309015 : long lx=lg(x), ly=lg(y);
733 13309015 : if (ly==1) return cgetg(1,t_MAT);
734 13309015 : if (lx==1) return zero_Flm(0, ly-1);
735 13309015 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
736 : }
737 :
738 : GEN
739 14059286 : Flm_mul(GEN x, GEN y, ulong p)
740 : {
741 14059286 : long lx=lg(x), ly=lg(y);
742 14059286 : if (ly==1) return cgetg(1,t_MAT);
743 14058999 : if (lx==1) return zero_Flm(0, ly-1);
744 14058999 : 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 5607 : FpC_FpV_mul(GEN x, GEN y, GEN p)
825 : {
826 5607 : long i,j, lx=lg(x), ly=lg(y);
827 : GEN z;
828 5607 : if (ly==1) return cgetg(1,t_MAT);
829 5607 : z = cgetg(ly,t_MAT);
830 63546 : for (j=1; j < ly; j++)
831 : {
832 57939 : GEN zj = cgetg(lx,t_COL);
833 1377866 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
834 57939 : gel(z, j) = zj;
835 : }
836 5607 : return z;
837 : }
838 :
839 : /* Multiply a line vector by a column and return a scalar (t_INT) */
840 : GEN
841 8640523 : FpV_dotproduct(GEN x, GEN y, GEN p)
842 : {
843 8640523 : long i, lx = lg(x);
844 : pari_sp av;
845 : GEN c;
846 8640523 : if (lx == 1) return gen_0;
847 8638017 : av = avma; c = mulii(gel(x,1),gel(y,1));
848 27991857 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
849 8638011 : 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 8998716 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
865 : {
866 8998716 : ulong c = uel(x,0) * uel(y,0);
867 : long k;
868 64307171 : for (k = 1; k < lx; k++) {
869 55308455 : c += uel(x,k) * uel(y,k);
870 55308455 : if (c & HIGHBIT) c %= p;
871 : }
872 8998716 : return c % p;
873 : }
874 :
875 : INLINE ulong
876 740068 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
877 : {
878 : ulong l0, l1, v1, h0, h1;
879 740068 : long i = 0;
880 : LOCAL_OVERFLOW;
881 : LOCAL_HIREMAINDER;
882 740068 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
883 15410510 : while (++i < lx) {
884 14670442 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
885 14670442 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
886 : }
887 740068 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
888 464435 : 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 9391887 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
915 : {
916 9391887 : long lx = minss(lgpol(x), lgpol(y));
917 9391874 : if (lx == 0) return 0;
918 9323300 : if (pi)
919 734510 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
920 : else
921 8588790 : 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 1657402 : FpM_FpC_mul(GEN x, GEN y, GEN p)
929 : {
930 1657402 : long lx = lg(x);
931 1657402 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
932 : }
933 : GEN
934 702440 : Flm_Flc_mul(GEN x, GEN y, ulong p)
935 : {
936 702440 : long l, lx = lg(x);
937 702440 : if (lx==1) return cgetg(1,t_VECSMALL);
938 702440 : l = lgcols(x);
939 702440 : if (p==2)
940 179884 : return Flm_Flc_mul_i_2(x, y, lx, l);
941 522556 : else if (SMALL_ULONG(p))
942 285098 : return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
943 : else
944 237458 : 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 7541685 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
966 : {
967 7541685 : long l, lx = lg(x);
968 : GEN z;
969 7541685 : if (lx==1) return pol0_Flx(sv);
970 7541685 : l = lgcols(x);
971 7540954 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
972 7538649 : if (SMALL_ULONG(p))
973 3307024 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
974 : else
975 4231625 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
976 7549346 : 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 1418753 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
982 : {
983 1418753 : long i, l, lx = lg(x);
984 : GEN z;
985 1418753 : if (lx==1) return pol_0(v);
986 1418753 : l = lgcols(x);
987 1418753 : z = new_chunk(l+1);
988 2160130 : for (i=l-1; i; i--)
989 : {
990 1976734 : pari_sp av = avma;
991 1976734 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
992 1976732 : p1 = modii(p1, p);
993 1976733 : if (signe(p1))
994 : {
995 1235356 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
996 1235356 : gel(z,i+1) = gerepileuptoint(av, p1);
997 1235356 : break;
998 : }
999 741377 : set_avma(av);
1000 : }
1001 1418752 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1002 1235356 : z[0] = evaltyp(t_POL) | _evallg(i+2);
1003 1235356 : z[1] = evalsigne(1) | evalvarn(v);
1004 3597575 : for (; i; i--)
1005 : {
1006 2362219 : pari_sp av = avma;
1007 2362219 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1008 2362217 : gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
1009 : }
1010 1235356 : return z;
1011 : }
1012 :
1013 : /********************************************************************/
1014 : /** **/
1015 : /** TRANSPOSITION **/
1016 : /** **/
1017 : /********************************************************************/
1018 :
1019 : /* == zm_transpose */
1020 : GEN
1021 714662 : Flm_transpose(GEN x)
1022 : {
1023 714662 : long i, dx, lx = lg(x);
1024 : GEN y;
1025 714662 : if (lx == 1) return cgetg(1,t_MAT);
1026 714529 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1027 8502018 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1028 714529 : 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 1385224 : matid_Flm(long n)
1070 : {
1071 1385224 : GEN y = cgetg(n+1,t_MAT);
1072 : long i;
1073 1385223 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1074 8746466 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1075 1385212 : 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 49307373 : ZV_to_Flv(GEN x, ulong p)
1094 511887709 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1095 :
1096 : GEN
1097 12293761 : ZM_to_Flm(GEN x, ulong p)
1098 60601793 : { 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 20251526 : 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 3443 : Fp_to_mod(GEN z, GEN p)
1112 : {
1113 3443 : retmkintmod(modii(z, p), icopy(p));
1114 : }
1115 :
1116 : GEN
1117 997248 : FpX_to_mod_raw(GEN z, GEN p)
1118 : {
1119 997248 : long i, l = lg(z);
1120 : GEN x;
1121 997248 : 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 474824 : x = cgetg(l,t_POL); x[1] = z[1];
1127 3918509 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1128 474824 : return normalizepol_lg(x,l);
1129 : }
1130 :
1131 : /* z in Z[X], return z * Mod(1,p), normalized*/
1132 : GEN
1133 1669085 : FpX_to_mod(GEN z, GEN p)
1134 : {
1135 1669085 : long i, l = lg(z);
1136 : GEN x;
1137 1669085 : if (l == 2)
1138 : {
1139 3171 : x = cgetg(3,t_POL); x[1] = z[1];
1140 3171 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1141 : }
1142 1665914 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1143 18370815 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1144 1665914 : return normalizepol_lg(x,l);
1145 : }
1146 :
1147 : GEN
1148 84021 : FpXC_to_mod(GEN z, GEN p)
1149 : {
1150 84021 : long i,l = lg(z);
1151 84021 : GEN x = cgetg(l, t_COL);
1152 84021 : p = icopy(p);
1153 474166 : for (i=1; i<l; i++)
1154 390145 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1155 84021 : return x;
1156 : }
1157 :
1158 : static GEN
1159 0 : FpXC_to_mod_raw(GEN x, GEN p)
1160 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1161 :
1162 : GEN
1163 0 : FpXM_to_mod(GEN z, GEN p)
1164 : {
1165 0 : long i,l = lg(z);
1166 0 : GEN x = cgetg(l, t_MAT);
1167 0 : p = icopy(p);
1168 0 : for (i=1; i<l; i++)
1169 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1170 0 : return x;
1171 : }
1172 :
1173 : /* z in Z^n, return z * Mod(1,p), normalized*/
1174 : GEN
1175 35935 : FpV_to_mod(GEN z, GEN p)
1176 : {
1177 35935 : long i,l = lg(z);
1178 35935 : GEN x = cgetg(l, t_VEC);
1179 35935 : if (l == 1) return x;
1180 35935 : p = icopy(p);
1181 72269 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1182 35935 : return x;
1183 : }
1184 : /* z in Z^n, return z * Mod(1,p), normalized*/
1185 : GEN
1186 105 : FpC_to_mod(GEN z, GEN p)
1187 : {
1188 105 : long i, l = lg(z);
1189 105 : GEN x = cgetg(l, t_COL);
1190 105 : if (l == 1) return x;
1191 91 : p = icopy(p);
1192 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1193 91 : return x;
1194 : }
1195 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1196 : GEN
1197 226 : FpM_to_mod(GEN z, GEN p)
1198 : {
1199 226 : long i, j, m, l = lg(z);
1200 226 : GEN x = cgetg(l,t_MAT), y, zi;
1201 226 : if (l == 1) return x;
1202 205 : m = lgcols(z);
1203 205 : p = icopy(p);
1204 2456 : for (i=1; i<l; i++)
1205 : {
1206 2251 : gel(x,i) = cgetg(m,t_COL);
1207 2251 : y = gel(x,i); zi= gel(z,i);
1208 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1209 : }
1210 205 : return x;
1211 : }
1212 : GEN
1213 28 : Flc_to_mod(GEN z, ulong pp)
1214 : {
1215 28 : long i, l = lg(z);
1216 28 : GEN p, x = cgetg(l, t_COL);
1217 28 : if (l == 1) return x;
1218 28 : p = utoipos(pp);
1219 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1220 28 : return x;
1221 : }
1222 : GEN
1223 56292 : Flm_to_mod(GEN z, ulong pp)
1224 : {
1225 56292 : long i, j, m, l = lg(z);
1226 56292 : GEN p, x = cgetg(l,t_MAT), y, zi;
1227 56292 : if (l == 1) return x;
1228 56264 : m = lgcols(z);
1229 56264 : p = utoipos(pp);
1230 179348 : for (i=1; i<l; i++)
1231 : {
1232 123084 : gel(x,i) = cgetg(m,t_COL);
1233 123084 : y = gel(x,i); zi= gel(z,i);
1234 506844 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1235 : }
1236 56264 : return x;
1237 : }
1238 :
1239 : GEN
1240 574 : FpVV_to_mod(GEN z, GEN p)
1241 : {
1242 574 : long i, j, m, l = lg(z);
1243 574 : GEN x = cgetg(l,t_VEC), y, zi;
1244 574 : if (l == 1) return x;
1245 574 : m = lgcols(z);
1246 574 : p = icopy(p);
1247 1246 : for (i=1; i<l; i++)
1248 : {
1249 672 : gel(x,i) = cgetg(m,t_VEC);
1250 672 : y = gel(x,i); zi= gel(z,i);
1251 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1252 : }
1253 574 : return x;
1254 : }
1255 :
1256 : /* z in Z^n, return z * Mod(1,p), normalized*/
1257 : GEN
1258 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1259 : {
1260 7 : long i,l = lg(z);
1261 7 : GEN x = cgetg(l, t_COL);
1262 7 : if (l == 1) return x;
1263 7 : p = icopy(p);
1264 7 : T = FpX_to_mod_raw(T, p);
1265 21 : for (i=1; i<l; i++)
1266 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1267 7 : return x;
1268 : }
1269 :
1270 : static GEN
1271 582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1272 : {
1273 582533 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1274 582533 : return mkpolmod(z, T);
1275 : }
1276 :
1277 : /* z in Z^n, return z * Mod(1,p), normalized*/
1278 : GEN
1279 28 : FqC_to_mod(GEN z, GEN T, GEN p)
1280 : {
1281 : GEN x;
1282 28 : long i,l = lg(z);
1283 28 : if (!T) return FpC_to_mod(z, p);
1284 28 : x = cgetg(l, t_COL);
1285 28 : if (l == 1) return x;
1286 28 : p = icopy(p);
1287 28 : T = FpX_to_mod_raw(T, p);
1288 504 : for (i=1; i<l; i++)
1289 476 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1290 28 : return x;
1291 : }
1292 :
1293 : /* z in Z^n, return z * Mod(1,p), normalized*/
1294 : static GEN
1295 107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1296 690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1297 :
1298 : /* z in Z^n, return z * Mod(1,p), normalized*/
1299 : GEN
1300 26145 : FqM_to_mod(GEN z, GEN T, GEN p)
1301 : {
1302 : GEN x;
1303 26145 : long i,l = lg(z);
1304 26145 : if (!T) return FpM_to_mod(z, p);
1305 26145 : x = cgetg(l, t_MAT);
1306 26145 : if (l == 1) return x;
1307 26117 : p = icopy(p);
1308 26117 : T = FpX_to_mod_raw(T, p);
1309 134092 : for (i=1; i<l; i++)
1310 107975 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1311 26117 : return x;
1312 : }
1313 :
1314 : /********************************************************************/
1315 : /* */
1316 : /* Blackbox linear algebra */
1317 : /* */
1318 : /********************************************************************/
1319 :
1320 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1321 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1322 : * (e_j) is the canonical basis.
1323 : * A sparse matrix (zMs) is a t_VEC of zCs */
1324 :
1325 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1326 : * integer representing an element of Fp. This is important since p can be
1327 : * large and p+E[i] would not fit in a C long. */
1328 :
1329 : /* RgCs and RgMs are similar, except that the type of the component is
1330 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1331 : * of E. */
1332 :
1333 : /* Most functions take an argument nbrow which is the number of lines of the
1334 : * column/matrix, which cannot be derived from the data. */
1335 :
1336 : GEN
1337 0 : zCs_to_ZC(GEN R, long nbrow)
1338 : {
1339 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1340 0 : long j, l = lg(C);
1341 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1342 0 : return c;
1343 : }
1344 :
1345 : GEN
1346 0 : zMs_to_ZM(GEN x, long nbrow)
1347 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1348 :
1349 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1350 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1351 : GEN
1352 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1353 : {
1354 147 : pari_sp ltop = avma;
1355 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1356 147 : if (ZV_equal0(B)) return zerocol(n);
1357 147 : while (++col <= n)
1358 : {
1359 147 : pari_sp btop = avma, av;
1360 : long i, lQ;
1361 147 : GEN V, Q, M, W = B;
1362 147 : GEN b = cgetg(m+2, t_POL);
1363 147 : b[1] = evalsigne(1)|evalvarn(0);
1364 147 : gel(b, 2) = gel(W, col);
1365 46225 : for (i = 3; i<m+2; i++)
1366 46078 : gel(b, i) = cgeti(lgefint(p));
1367 147 : av = avma;
1368 46225 : for (i = 3; i<m+2; i++)
1369 : {
1370 46078 : W = f(E, W);
1371 46078 : affii(gel(W, col),gel(b, i));
1372 46078 : if (gc_needed(av,1))
1373 : {
1374 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1375 72 : W = gerepileupto(av, W);
1376 : }
1377 : }
1378 147 : b = FpX_renormalize(b, m+2);
1379 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1380 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1381 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1382 147 : W = B; lQ =lg(Q);
1383 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1384 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1385 147 : av = avma;
1386 22745 : for (i = lQ-3; i > 1; i--)
1387 : {
1388 22598 : W = f(E, W);
1389 22598 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1390 22598 : if (gc_needed(av,1))
1391 : {
1392 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1393 110 : gerepileall(av, 2, &V, &W);
1394 : }
1395 : }
1396 147 : V = FpC_red(V, p);
1397 147 : W = FpC_sub(f(E,V), B, p);
1398 147 : if (ZV_equal0(W)) return gerepilecopy(ltop, V);
1399 0 : av = avma;
1400 0 : for (i = 1; i <= n; ++i)
1401 : {
1402 0 : V = W;
1403 0 : W = f(E, V);
1404 0 : if (ZV_equal0(W))
1405 0 : return gerepilecopy(ltop, shallowtrans(V));
1406 0 : gerepileall(av, 2, &V, &W);
1407 : }
1408 0 : set_avma(btop);
1409 : }
1410 0 : return NULL;
1411 : }
1412 :
1413 : GEN
1414 0 : zMs_ZC_mul(GEN M, GEN B)
1415 : {
1416 : long i, j;
1417 0 : long n = lg(B)-1;
1418 0 : GEN V = zerocol(n);
1419 0 : for (i = 1; i <= n; ++i)
1420 0 : if (signe(gel(B, i)))
1421 : {
1422 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1423 0 : long l = lg(C);
1424 0 : for (j = 1; j < l; ++j)
1425 : {
1426 0 : long k = C[j];
1427 0 : switch(E[j])
1428 : {
1429 0 : case 1:
1430 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1431 0 : break;
1432 0 : case -1:
1433 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1434 0 : break;
1435 0 : default:
1436 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]));
1437 0 : break;
1438 : }
1439 : }
1440 : }
1441 0 : return V;
1442 : }
1443 :
1444 : GEN
1445 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1446 :
1447 : GEN
1448 68970 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1449 : {
1450 68970 : long i, j, lM = lg(M);
1451 68970 : GEN V = cgetg(lM,t_VEC);
1452 28332699 : for (i = 1; i < lM; ++i)
1453 : {
1454 28263729 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1455 28263729 : pari_sp av = avma;
1456 28263729 : long lC = lg(C);
1457 28263729 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1458 28263729 : z = mulis(gel(B, C[1]), E[1]);
1459 228730226 : for (j = 2; j < lC; ++j)
1460 : {
1461 200466497 : GEN b = gel(B, C[j]);
1462 200466497 : switch(E[j])
1463 : {
1464 131566519 : case 1: z = addii(z, b); break;
1465 57524496 : case -1: z = subii(z, b); break;
1466 11375482 : default: z = addii(z, mulis(b, E[j])); break;
1467 : }
1468 : }
1469 28263729 : gel(V,i) = gerepileuptoint(av, p? Fp_red(z, p): z);
1470 : }
1471 68970 : return V;
1472 : }
1473 : GEN
1474 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1475 :
1476 : GEN
1477 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1478 : {
1479 0 : pari_sp av = avma, av2;
1480 0 : GEN xi, xb, pi = gen_1, P;
1481 : long i;
1482 0 : if (!C) {
1483 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1484 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1485 : }
1486 0 : P = utoipos(p);
1487 0 : av2 = avma;
1488 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1489 0 : xb = Flm_to_ZM(xi);
1490 0 : for (i = 2; i <= e; i++)
1491 : {
1492 0 : pi = muliu(pi, p); /* = p^(i-1) */
1493 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1494 0 : if (gc_needed(av,2))
1495 : {
1496 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1497 0 : gerepileall(av2,3, &pi,&b,&xb);
1498 : }
1499 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1500 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1501 : }
1502 0 : return gerepileupto(av, xb);
1503 : }
1504 :
1505 : struct wrapper_modp_s {
1506 : GEN (*f)(void*, GEN);
1507 : void *E;
1508 : GEN p;
1509 : };
1510 :
1511 : /* compute f(x) mod p */
1512 : static GEN
1513 0 : wrap_relcomb_modp(void *E, GEN x)
1514 : {
1515 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1516 0 : return FpC_red(W->f(W->E, x), W->p);
1517 : }
1518 : static GEN
1519 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1520 :
1521 : static GEN
1522 68823 : wrap_relker(void*E, GEN x)
1523 : {
1524 68823 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1525 68823 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1526 : }
1527 :
1528 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1529 : GEN
1530 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1531 : {
1532 : struct wrapper_modp_s W;
1533 0 : pari_sp av = avma;
1534 0 : GEN xb, xi, pi = gen_1;
1535 : long i;
1536 0 : W.E = E;
1537 0 : W.f = f;
1538 0 : W.p = p;
1539 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1540 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1541 0 : xb = xi;
1542 0 : for (i = 2; i <= e; i++)
1543 : {
1544 0 : pi = mulii(pi, p); /* = p^(i-1) */
1545 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1546 0 : if (gc_needed(av,2))
1547 : {
1548 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1549 0 : gerepileall(av,3, &pi,&B,&xb);
1550 : }
1551 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1552 0 : if (!xi) return NULL;
1553 0 : if (typ(xi) == t_VEC) return gerepileupto(av, xi);
1554 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1555 : }
1556 0 : return gerepileupto(av, xb);
1557 : }
1558 :
1559 : static GEN
1560 23039 : vecprow(GEN A, GEN prow)
1561 : {
1562 23039 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1563 : }
1564 :
1565 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1566 : * or the index of a column which is linearly dependent from the others as a
1567 : * t_VECSMALL with a single component. */
1568 : GEN
1569 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1570 : {
1571 0 : pari_sp av = avma;
1572 : GEN pcol, prow;
1573 0 : long nbi=lg(M)-1, lR;
1574 : long i, n;
1575 : GEN Mp, Ap, Rp;
1576 : pari_timer ti;
1577 0 : if (DEBUGLEVEL) timer_start(&ti);
1578 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1579 0 : if (!pcol) return gc_NULL(av);
1580 0 : if (DEBUGLEVEL)
1581 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1582 0 : n = lg(pcol)-1;
1583 0 : Mp = cgetg(n+1, t_MAT);
1584 0 : for(i=1; i<=n; i++)
1585 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1586 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1587 0 : if (DEBUGLEVEL) timer_start(&ti);
1588 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1589 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1590 0 : if (!Rp) return gc_NULL(av);
1591 0 : lR = lg(Rp)-1;
1592 0 : if (typ(Rp) == t_COL)
1593 : {
1594 0 : GEN R = zerocol(nbi+1);
1595 0 : for (i=1; i<=lR; i++)
1596 0 : gel(R,pcol[i]) = gel(Rp,i);
1597 0 : return gerepilecopy(av, R);
1598 : }
1599 0 : for (i = 1; i <= lR; ++i)
1600 0 : if (signe(gel(Rp, i)))
1601 0 : return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
1602 : return NULL; /* LCOV_EXCL_LINE */
1603 : }
1604 :
1605 : GEN
1606 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1607 : {
1608 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1609 : }
1610 :
1611 : GEN
1612 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1613 : {
1614 : GEN res;
1615 0 : pari_CATCH(e_INV)
1616 : {
1617 0 : GEN E = pari_err_last();
1618 0 : GEN x = err_get_compo(E,2);
1619 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1620 0 : if (DEBUGLEVEL)
1621 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1622 0 : res = NULL;
1623 : } pari_TRY {
1624 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1625 0 : } pari_ENDCATCH
1626 0 : return res;
1627 : }
1628 :
1629 : static GEN
1630 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1631 : {
1632 147 : long i, j, oldf = 0, f = 0;
1633 147 : long lrow = lg(prow), lM = lg(M);
1634 147 : GEN W = const_vecsmall(lM-1,1);
1635 147 : GEN R = cgetg(lrow, t_VEC);
1636 228354 : for (i=1; i<lrow; i++)
1637 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1638 : do
1639 : {
1640 442 : oldf = f;
1641 374995 : for(i=1; i<lM; i++)
1642 374553 : if (W[i])
1643 : {
1644 131759 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1645 131759 : long c=0, cj=0, lF = lg(F);
1646 1093902 : for(j=1; j<lF; j++)
1647 962143 : if (!gel(R,F[j])) { c++; cj=j; }
1648 131759 : if (c>=2) continue;
1649 112311 : if (c==1)
1650 : {
1651 32896 : pari_sp av = avma;
1652 32896 : GEN s = gen_0;
1653 272945 : for(j=1; j<lF; j++)
1654 240049 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1655 : /* s /= -E[cj] mod p */
1656 32896 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1657 32896 : gel(R,F[cj]) = gerepileuptoint(av, s);
1658 32896 : f++;
1659 : }
1660 112311 : W[i]=0;
1661 : }
1662 442 : } while(oldf!=f);
1663 228354 : for (i=1; i<lrow; i++)
1664 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1665 147 : return R;
1666 : }
1667 :
1668 : /* Return a linear form Y such that YM is essentially 0 */
1669 : GEN
1670 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1671 : {
1672 147 : pari_sp av = avma, av2;
1673 : GEN Mp, pcol, prow;
1674 : long i, n;
1675 : pari_timer ti;
1676 : struct wrapper_modp_s W;
1677 147 : if (DEBUGLEVEL) timer_start(&ti);
1678 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1679 147 : if (!pcol) return gc_NULL(av);
1680 147 : if (DEBUGLEVEL)
1681 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1682 147 : n = lg(pcol)-1;
1683 147 : Mp = cgetg(n+1, t_MAT);
1684 23186 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1685 147 : W.E = (void*)Mp;
1686 147 : W.p = p;
1687 147 : av2 = avma;
1688 0 : for(;; set_avma(av2))
1689 0 : {
1690 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1691 147 : if (DEBUGLEVEL) timer_start(&ti);
1692 147 : pari_CATCH(e_INV)
1693 : {
1694 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1695 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1696 0 : if (DEBUGLEVEL)
1697 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1698 0 : Rp = NULL;
1699 : } pari_TRY {
1700 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1701 147 : } pari_ENDCATCH
1702 147 : if (!Rp) continue;
1703 147 : if (typ(Rp)==t_COL)
1704 : {
1705 147 : Rp = FpC_sub(Rp,B,p);
1706 147 : if (ZV_equal0(Rp)) continue;
1707 : }
1708 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1709 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1710 147 : return gerepilecopy(av, R);
1711 : }
1712 : }
1713 :
1714 : GEN
1715 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1716 : {
1717 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1718 : }
|
