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 29832102 : FpC_red(GEN x, GEN p)
28 261375836 : { 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 6865034 : FpC_center(GEN x, GEN p, GEN pov2)
37 43335105 : { 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 5194309 : FpM_red(GEN x, GEN p)
46 24399128 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2728252 : FpM_center(GEN x, GEN p, GEN pov2)
50 8273169 : { pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
51 :
52 : /* HACK: assume p > 3, which ensures u = gen_[0-2] is never written to */
53 : static void
54 168 : Fp_center_inplace(GEN u, GEN p, GEN ps2)
55 : {
56 168 : if (abscmpii(u, ps2) > 0)
57 : {
58 28 : pari_sp av = avma;
59 28 : affii(subii(u, p), u);
60 28 : set_avma(av);
61 : }
62 168 : }
63 :
64 : /* p > 3; assume entries in [0,p[ and ps2 = p>>1. */
65 : static void
66 42 : _FpC_center_inplace(GEN z, GEN p, GEN ps2)
67 : {
68 42 : long i, l = lg(z);
69 210 : for (i = 1; i < l; i++) Fp_center_inplace(gel(z,i), p, ps2);
70 42 : }
71 : static void
72 7 : _F3C_center_inplace(GEN z)
73 : {
74 7 : long i, l = lg(z);
75 35 : for (i = 1; i < l; i++) /* z[i] = 0, 1 : no-op */
76 28 : if (equaliu(gel(z,i), 2)) gel(z,i) = gen_m1;
77 7 : }
78 : void
79 0 : FpC_center_inplace(GEN z, GEN p, GEN ps2)
80 : {
81 0 : switch (itou_or_0(p))
82 : {
83 0 : case 2: break;
84 0 : case 3:_F3C_center_inplace(z); break;
85 0 : default: _FpC_center_inplace(z, p, ps2);
86 : }
87 0 : }
88 :
89 : void
90 42 : FpM_center_inplace(GEN z, GEN p, GEN pov2)
91 : {
92 42 : long i, l = lg(z);
93 42 : switch (itou_or_0(p))
94 : {
95 21 : case 2: break;
96 7 : case 3:
97 14 : for (i = 1; i < l; i++) _F3C_center_inplace(gel(z,i));
98 7 : break;
99 14 : default:
100 56 : for (i = 1; i < l; i++) _FpC_center_inplace(gel(z,i), p, pov2);
101 : }
102 42 : }
103 : GEN
104 11347 : Flm_center(GEN x, ulong p, ulong ps2)
105 175875 : { pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
106 :
107 : GEN
108 147 : random_FpV(long d, GEN p)
109 : {
110 : long i;
111 147 : GEN y = cgetg(d+1,t_VEC);
112 23183 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
113 147 : return y;
114 : }
115 :
116 : GEN
117 924 : random_FpC(long d, GEN p)
118 : {
119 : long i;
120 924 : GEN y = cgetg(d+1,t_COL);
121 6188 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
122 924 : return y;
123 : }
124 :
125 : GEN
126 41412 : random_Flv(long n, ulong p)
127 : {
128 41412 : GEN y = cgetg(n+1, t_VECSMALL);
129 : long i;
130 262142 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
131 41412 : return y;
132 : }
133 :
134 : /********************************************************************/
135 : /** **/
136 : /** ADD, SUB **/
137 : /** **/
138 : /********************************************************************/
139 : GEN
140 283601 : FpC_add(GEN x, GEN y, GEN p)
141 3551393 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
142 :
143 : GEN
144 0 : FpV_add(GEN x, GEN y, GEN p)
145 0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
146 :
147 : GEN
148 0 : FpM_add(GEN x, GEN y, GEN p)
149 0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
150 :
151 : GEN
152 1812476 : Flv_add(GEN x, GEN y, ulong p)
153 : {
154 1812476 : if (p==2)
155 7096824 : pari_APPLY_ulong(x[i]^y[i])
156 : else
157 41998925 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
158 : }
159 :
160 : void
161 446313 : Flv_add_inplace(GEN x, GEN y, ulong p)
162 : {
163 446313 : long i, l = lg(x);
164 446313 : if (p==2)
165 1387785 : for (i = 1; i < l; i++) x[i] ^= y[i];
166 : else
167 164878 : for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
168 446313 : }
169 :
170 : ulong
171 3871 : Flv_sum(GEN x, ulong p)
172 : {
173 3871 : long i, l = lg(x);
174 3871 : ulong s = 0;
175 3871 : if (p==2)
176 17689 : for (i = 1; i < l; i++) s ^= x[i];
177 : else
178 0 : for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
179 3871 : return s;
180 : }
181 :
182 : GEN
183 965634 : FpC_sub(GEN x, GEN y, GEN p)
184 18959595 : { pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
185 :
186 : GEN
187 0 : FpV_sub(GEN x, GEN y, GEN p)
188 0 : { pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
189 :
190 : GEN
191 0 : FpM_sub(GEN x, GEN y, GEN p)
192 0 : { pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
193 :
194 : GEN
195 231499901 : Flv_sub(GEN x, GEN y, ulong p)
196 1037115583 : { pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
197 :
198 : void
199 0 : Flv_sub_inplace(GEN x, GEN y, ulong p)
200 : {
201 0 : long i, l = lg(x);
202 0 : for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
203 0 : }
204 :
205 : GEN
206 51208 : Flm_Fl_add(GEN x, ulong y, ulong p)
207 : {
208 51208 : long l = lg(x), i, j;
209 51208 : GEN z = cgetg(l,t_MAT);
210 :
211 51208 : if (l==1) return z;
212 51208 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
213 241347 : for (i=1; i<l; i++)
214 : {
215 190139 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
216 190139 : gel(z,i) = zi;
217 1331328 : for (j=1; j<l; j++) zi[j] = xi[j];
218 190139 : zi[i] = Fl_add(zi[i], y, p);
219 : }
220 51208 : return z;
221 : }
222 : GEN
223 3724 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
224 :
225 : GEN
226 29322 : Flm_add(GEN x, GEN y, ulong p)
227 695951 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
228 :
229 : GEN
230 25774417 : Flm_sub(GEN x, GEN y, ulong p)
231 257268568 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
232 :
233 : /********************************************************************/
234 : /** **/
235 : /** MULTIPLICATION **/
236 : /** **/
237 : /********************************************************************/
238 : GEN
239 960456 : FpC_Fp_mul(GEN x, GEN y, GEN p)
240 11704192 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
241 :
242 : GEN
243 1753334 : Flv_Fl_mul(GEN x, ulong y, ulong p)
244 41052004 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
245 :
246 : GEN
247 4692 : Flv_Fl_div(GEN x, ulong y, ulong p)
248 : {
249 4692 : return Flv_Fl_mul(x, Fl_inv(y, p), p);
250 : }
251 : void
252 0 : Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
253 : {
254 0 : Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
255 0 : }
256 : GEN
257 765487 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
258 765487 : long i, j, h, l = lg(X);
259 765487 : GEN A = cgetg(l, t_MAT);
260 765487 : if (l == 1) return A;
261 765487 : h = lgcols(X);
262 4076735 : for (j=1; j<l; j++)
263 : {
264 3311412 : GEN a = cgetg(h, t_COL), x = gel(X, j);
265 24386427 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
266 3311248 : gel(A,j) = a;
267 : }
268 765323 : return A;
269 : }
270 :
271 : /* x *= y */
272 : void
273 7776228 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
274 : {
275 : long i;
276 7776228 : if (HIGHWORD(y | p))
277 8598712 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
278 : else
279 32083762 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
280 7776229 : }
281 : void
282 0 : Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
283 0 : { Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
284 :
285 : /* set y *= x */
286 : void
287 0 : Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
288 : {
289 0 : long i, j, m, l = lg(y);
290 0 : if (l == 1) return;
291 0 : m = lgcols(y);
292 0 : if (HIGHWORD(x | p))
293 0 : for(j=1; j<l; j++)
294 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
295 : else
296 0 : for(j=1; j<l; j++)
297 0 : for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 17159681 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
303 : {
304 17159681 : long i, j, m, l = lg(y);
305 17159681 : GEN z = cgetg(l, t_MAT);
306 17158819 : if (l == 1) return z;
307 17158819 : m = lgcols(y);
308 153875699 : for(j=1; j<l; j++) {
309 136654657 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 367516655 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
311 : }
312 17221042 : return z;
313 : }
314 :
315 : /* return x * y */
316 : static GEN
317 8641886 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
318 : {
319 8641886 : long i, j, m, l = lg(y);
320 8641886 : GEN z = cgetg(l, t_MAT);
321 8641878 : if (l == 1) return z;
322 8641878 : m = lgcols(y);
323 53933946 : for(j=1; j<l; j++) {
324 45292074 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
325 163459509 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
326 : }
327 8641872 : return z;
328 : }
329 :
330 : /* return x * y */
331 : GEN
332 25718849 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
333 : {
334 25718849 : if (HIGHWORD(p))
335 17159575 : return Flm_Fl_mul_pre_i(y, x, p, pi);
336 : else
337 8559274 : return Flm_Fl_mul_OK(y, x, p);
338 : }
339 :
340 : /* return x * y */
341 : GEN
342 82649 : Flm_Fl_mul(GEN y, ulong x, ulong p)
343 : {
344 82649 : if (HIGHWORD(p))
345 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
346 : else
347 82575 : return Flm_Fl_mul_OK(y, x, p);
348 : }
349 :
350 : GEN
351 4102230 : Flv_neg(GEN x, ulong p)
352 31299717 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
353 :
354 : void
355 6419 : Flv_neg_inplace(GEN v, ulong p)
356 : {
357 : long i;
358 184519 : for (i = 1; i < lg(v); ++i)
359 178100 : v[i] = Fl_neg(v[i], p);
360 6419 : }
361 :
362 : GEN
363 324925 : Flm_neg(GEN x, ulong p)
364 4396410 : { pari_APPLY_same(Flv_neg(gel(x,i), p)) }
365 :
366 : /* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
367 : static GEN
368 18918687 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
369 : {
370 18918687 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
371 : long k;
372 230075399 : for (k = 2; k < lx; k++)
373 : {
374 211158038 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
375 211156014 : if (signe(t)) c = addii(c, t);
376 : }
377 18917361 : return c;
378 : }
379 :
380 : static long
381 97688616 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
382 : {
383 : long k;
384 97688616 : long c = coeff(x,i,1) * y[1];
385 1498281691 : for (k = 2; k < lx; k++)
386 1400593075 : c += coeff(x,i,k) * y[k];
387 97688616 : return c;
388 : }
389 :
390 : GEN
391 6450565 : zm_zc_mul(GEN x, GEN y)
392 : {
393 6450565 : long lx = lg(x), l, i;
394 : GEN z;
395 6450565 : if (lx == 1) return cgetg(1, t_VECSMALL);
396 6450565 : l = lg(gel(x,1));
397 6450565 : z = cgetg(l,t_VECSMALL);
398 104139181 : for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
399 6450565 : return z;
400 : }
401 :
402 : GEN
403 2114 : zm_mul(GEN x, GEN y)
404 : {
405 2114 : long i,j,lx=lg(x), ly=lg(y);
406 : GEN z;
407 2114 : if (ly==1) return cgetg(1,t_MAT);
408 2114 : z = cgetg(ly,t_MAT);
409 2114 : if (lx==1)
410 : {
411 0 : for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
412 0 : return z;
413 : }
414 30849 : for (j=1; j<ly; j++)
415 28735 : gel(z,j) = zm_zc_mul(x, gel(y,j));
416 2114 : return z;
417 : }
418 :
419 : static ulong
420 767884757 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
421 : {
422 767884757 : ulong c = ucoeff(x,i,1) * uel(y,1);
423 : long k;
424 12009370574 : for (k = 2; k < lx; k++) {
425 11241485817 : c += ucoeff(x,i,k) * uel(y,k);
426 11241485817 : if (c & HIGHBIT) c %= p;
427 : }
428 767884757 : return c % p;
429 : }
430 :
431 : static ulong
432 533629576 : Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
433 : {
434 : ulong l0, l1, v1, h0, h1;
435 533629576 : long k = 1;
436 : LOCAL_OVERFLOW;
437 : LOCAL_HIREMAINDER;
438 533629576 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
439 8862862008 : while (++k < lx) {
440 8329232432 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
441 8329232432 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
442 : }
443 533629576 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
444 55070678 : else return remlll_pre(v1, h1, l1, p, pi);
445 : }
446 :
447 : static GEN
448 329856 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
449 : {
450 : long i,j;
451 329856 : GEN z = NULL;
452 :
453 3725922 : for (j=1; j<lx; j++)
454 : {
455 3396066 : if (!y[j]) continue;
456 1106401 : if (!z) z = Flv_copy(gel(x,j));
457 29798028 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
458 : }
459 329856 : if (!z) z = zero_zv(l-1);
460 329856 : return z;
461 : }
462 :
463 : static GEN
464 1754822 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
465 : {
466 1754822 : GEN z = cgetg(l,t_COL);
467 : long i;
468 16781938 : for (i = 1; i < l; i++)
469 : {
470 15027166 : pari_sp av = avma;
471 15027166 : GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
472 15026946 : gel(z,i) = gc_INT(av, modii(c,p));
473 : }
474 1754772 : return z;
475 : }
476 :
477 : static void
478 86251338 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
479 : {
480 : long i;
481 854152152 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
482 86265251 : }
483 : static GEN
484 82960323 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
485 : {
486 82960323 : GEN z = cgetg(l,t_VECSMALL);
487 82959446 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
488 82960426 : return z;
489 : }
490 :
491 : static void
492 95037839 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
493 : {
494 : long i;
495 631495337 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
496 95088003 : }
497 : static GEN
498 90851508 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
499 : {
500 90851508 : GEN z = cgetg(l,t_VECSMALL);
501 90768410 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
502 90897010 : return z;
503 : }
504 :
505 : GEN
506 1634247 : FpM_mul(GEN x, GEN y, GEN p)
507 : {
508 1634247 : long lx=lg(x), ly=lg(y);
509 : GEN z;
510 1634247 : pari_sp av = avma;
511 1634247 : if (ly==1) return cgetg(1,t_MAT);
512 1634247 : if (lx==1) return zeromat(0, ly-1);
513 1634247 : if (lgefint(p) == 3)
514 : {
515 1588031 : ulong pp = uel(p,2);
516 1588031 : if (pp == 2)
517 : {
518 711599 : x = ZM_to_F2m(x);
519 711646 : y = ZM_to_F2m(y);
520 711654 : z = F2m_to_ZM(F2m_mul(x,y));
521 : }
522 : else
523 : {
524 876432 : x = ZM_to_Flm(x, pp);
525 876439 : y = ZM_to_Flm(y, pp);
526 876441 : z = Flm_to_ZM(Flm_mul(x,y, pp));
527 : }
528 : } else
529 46216 : z = FpM_red(ZM_mul(x, y), p);
530 1634305 : return gc_upto(av, z);
531 : }
532 :
533 : static GEN
534 28037968 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
535 : {
536 : long j;
537 28037968 : GEN z = cgetg(ly, t_MAT);
538 28037190 : if (SMALL_ULONG(p))
539 102481849 : for (j=1; j<ly; j++)
540 82450414 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
541 : else
542 98649006 : for (j=1; j<ly; j++)
543 90645417 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
544 28035024 : return z;
545 : }
546 :
547 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
548 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
549 : static void
550 66304 : add_slices_ip(long m, long n,
551 : GEN A, long ma, long da, long na, long ea,
552 : GEN B, long mb, long db, long nb, long eb,
553 : GEN M, long dy, long dx, ulong p)
554 : {
555 66304 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
556 : GEN C;
557 :
558 2782374 : for (j = 1; j <= min_e; j++) {
559 2716070 : C = gel(M, j + dx) + dy;
560 123921128 : for (i = 1; i <= min_d; i++)
561 121205058 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
562 121204921 : ucoeff(B, mb + i, nb + j), p);
563 2780767 : for (; i <= da; i++)
564 64560 : uel(C, i) = ucoeff(A, ma + i, na + j);
565 2716207 : for (; i <= db; i++)
566 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
567 2716207 : for (; i <= m; i++)
568 0 : uel(C, i) = 0;
569 : }
570 69917 : for (; j <= ea; j++) {
571 3613 : C = gel(M, j + dx) + dy;
572 195159 : for (i = 1; i <= da; i++)
573 191546 : uel(C, i) = ucoeff(A, ma + i, na + j);
574 3613 : for (; i <= m; i++)
575 0 : uel(C, i) = 0;
576 : }
577 66304 : for (; j <= eb; j++) {
578 0 : C = gel(M, j + dx) + dy;
579 0 : for (i = 1; i <= db; i++)
580 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
581 0 : for (; i <= m; i++)
582 0 : uel(C, i) = 0;
583 : }
584 66304 : for (; j <= n; j++) {
585 0 : C = gel(M, j + dx) + dy;
586 0 : for (i = 1; i <= m; i++)
587 0 : uel(C, i) = 0;
588 : }
589 66304 : }
590 :
591 : static GEN
592 33152 : add_slices(long m, long n,
593 : GEN A, long ma, long da, long na, long ea,
594 : GEN B, long mb, long db, long nb, long eb, ulong p)
595 : {
596 : GEN M;
597 : long j;
598 33152 : M = cgetg(n + 1, t_MAT);
599 1415694 : for (j = 1; j <= n; j++)
600 1382543 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
601 33151 : add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
602 33152 : return M;
603 : }
604 :
605 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
606 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
607 : static GEN
608 58016 : subtract_slices(long m, long n,
609 : GEN A, long ma, long da, long na, long ea,
610 : GEN B, long mb, long db, long nb, long eb, ulong p)
611 : {
612 58016 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
613 58015 : GEN M = cgetg(n + 1, t_MAT), C;
614 :
615 2527861 : for (j = 1; j <= min_e; j++) {
616 2469845 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
617 116128048 : for (i = 1; i <= min_d; i++)
618 113658181 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
619 113658267 : ucoeff(B, mb + i, nb + j), p);
620 2570226 : for (; i <= da; i++)
621 100445 : uel(C, i) = ucoeff(A, ma + i, na + j);
622 2544825 : for (; i <= db; i++)
623 75044 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
624 2469781 : for (; i <= m; i++)
625 0 : uel(C, i) = 0;
626 : }
627 58016 : for (; j <= ea; j++) {
628 0 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
629 0 : for (i = 1; i <= da; i++)
630 0 : uel(C, i) = ucoeff(A, ma + i, na + j);
631 0 : for (; i <= m; i++)
632 0 : uel(C, i) = 0;
633 : }
634 59435 : for (; j <= eb; j++) {
635 1419 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
636 82481 : for (i = 1; i <= db; i++)
637 81062 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
638 1419 : for (; i <= m; i++)
639 0 : uel(C, i) = 0;
640 : }
641 59435 : for (; j <= n; j++)
642 1419 : gel(M, j) = zero_Flv(m);
643 58016 : return M;
644 : }
645 :
646 : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
647 :
648 : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
649 : static GEN
650 8288 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
651 : {
652 : pari_sp av;
653 8288 : long m1 = (m + 1)/2, m2 = m/2,
654 8288 : n1 = (n + 1)/2, n2 = n/2,
655 8288 : p1 = (p + 1)/2, p2 = p/2;
656 : GEN A11, A12, A22, B11, B21, B22,
657 : S1, S2, S3, S4, T1, T2, T3, T4,
658 : M1, M2, M3, M4, M5, M6, M7,
659 : V1, V2, V3, C;
660 : long j;
661 8288 : C = cgetg(p + 1, t_MAT);
662 677039 : for (j = 1; j <= p; j++)
663 668751 : gel(C, j) = cgetg(m + 1, t_VECSMALL);
664 8288 : av = avma;
665 8288 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
666 8288 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
667 8288 : M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
668 8288 : if (gc_needed(av, 1))
669 0 : (void)gc_all(av, 2, &T2, &M2); /* destroy S1 */
670 8288 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
671 8288 : if (gc_needed(av, 1))
672 0 : (void)gc_all(av, 2, &M2, &T3); /* destroy T2 */
673 8288 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
674 8288 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
675 8288 : M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
676 8288 : if (gc_needed(av, 1))
677 0 : (void)gc_all(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
678 8288 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
679 8288 : if (gc_needed(av, 1))
680 0 : (void)gc_all(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
681 8288 : A11 = matslice(A, 1, m1, 1, n1);
682 8288 : B11 = matslice(B, 1, n1, 1, p1);
683 8288 : M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
684 8288 : if (gc_needed(av, 1))
685 0 : (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
686 8288 : A12 = matslice(A, 1, m1, n1 + 1, n);
687 8288 : B21 = matslice(B, n1 + 1, n, 1, p1);
688 8288 : M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
689 8288 : if (gc_needed(av, 1))
690 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
691 8288 : add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
692 8288 : if (gc_needed(av, 1))
693 0 : (void)gc_all(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy M4 */
694 8288 : M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
695 8288 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
696 8288 : if (gc_needed(av, 1))
697 0 : (void)gc_all(av, 6, &M2, &T3, &M3, &M1, &M5, &S4); /* destroy S3 */
698 8288 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
699 8288 : if (gc_needed(av, 1))
700 0 : (void)gc_all(av, 6, &M2, &M3, &M1, &M5, &S4, &T4); /* destroy T3 */
701 8288 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
702 8288 : if (gc_needed(av, 1))
703 0 : (void)gc_all(av, 5, &M2, &M3, &S4, &T4, &V1); /* destroy M1, M5 */
704 8288 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
705 8288 : M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
706 8288 : if (gc_needed(av, 1))
707 0 : (void)gc_all(av, 5, &M2, &M3, &T4, &V1, &M6); /* destroy S4, B22 */
708 8288 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
709 8288 : M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
710 8288 : if (gc_needed(av, 1))
711 0 : (void)gc_all(av, 5, &M2, &M3, &V1, &M6, &M7); /* destroy A22, T4 */
712 8288 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
713 8288 : add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
714 8288 : if (gc_needed(av, 1))
715 0 : (void)gc_all(av, 4, &M2, &M3, &V1, &M7); /* destroy V3, M6 */
716 8288 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
717 8288 : if (gc_needed(av, 1))
718 0 : (void)gc_all(av, 3, &M3, &M7, &V2); /* destroy V1, M2 */
719 8288 : add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
720 8288 : add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
721 8288 : set_avma(av); return C;
722 : }
723 :
724 : /* Strassen-Winograd used for dim >= ZM_sw_bound */
725 : static GEN
726 28044973 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
727 : {
728 28044973 : ulong e = expu(p);
729 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
730 23551993 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
731 : #else
732 4493486 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
733 : #endif
734 28045479 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
735 28037191 : return Flm_mul_classical(x, y, l, lx, ly, p, pi);
736 : else
737 8288 : return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
738 : }
739 :
740 : GEN
741 13798750 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
742 : {
743 13798750 : long lx=lg(x), ly=lg(y);
744 13798750 : if (ly==1) return cgetg(1,t_MAT);
745 13798750 : if (lx==1) return zero_Flm(0, ly-1);
746 13798750 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
747 : }
748 :
749 : GEN
750 14094637 : Flm_mul(GEN x, GEN y, ulong p)
751 : {
752 14094637 : long lx=lg(x), ly=lg(y);
753 14094637 : if (ly==1) return cgetg(1,t_MAT);
754 14094350 : if (lx==1) return zero_Flm(0, ly-1);
755 14094350 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
756 : }
757 :
758 : GEN
759 92727 : Flm_sqr(GEN x, ulong p)
760 : {
761 92727 : long lx = lg(x);
762 92727 : if (lx==1) return cgetg(1,t_MAT);
763 92727 : return Flm_mul_i(x, x, lx, lx, lx, p, get_Fl_red(p));
764 : }
765 :
766 : struct _Flm
767 : {
768 : ulong p;
769 : long n;
770 : };
771 :
772 : static GEN
773 21403 : _Flm_mul(void *E , GEN x, GEN y)
774 21403 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
775 : static GEN
776 92727 : _Flm_sqr(void *E, GEN x)
777 92727 : { return Flm_sqr(x,((struct _Flm*)E)->p); }
778 : static GEN
779 847 : _Flm_one(void *E)
780 847 : { return matid_Flm(((struct _Flm*)E)->n); }
781 : GEN
782 51360 : Flm_powu(GEN x, ulong n, ulong p)
783 : {
784 : struct _Flm d;
785 51360 : if (!n) return matid(lg(x)-1);
786 51360 : d.p = p;
787 51360 : return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
788 : }
789 : GEN
790 847 : Flm_powers(GEN x, ulong n, ulong p)
791 : {
792 : struct _Flm d;
793 847 : d.p = p;
794 847 : d.n = lg(x)-1;
795 847 : return gen_powers(x, n, 1, (void*)&d,
796 : &_Flm_sqr, &_Flm_mul, &_Flm_one);
797 : }
798 :
799 : static GEN
800 0 : _FpM_mul(void *p , GEN x, GEN y)
801 0 : { return FpM_mul(x,y,(GEN)p); }
802 : static GEN
803 0 : _FpM_sqr(void *p, GEN x)
804 0 : { return FpM_mul(x,x,(GEN)p); }
805 : GEN
806 0 : FpM_powu(GEN x, ulong n, GEN p)
807 : {
808 0 : if (!n) return matid(lg(x)-1);
809 0 : if (lgefint(p) == 3)
810 : {
811 0 : pari_sp av = avma;
812 0 : ulong pp = uel(p,2);
813 : GEN z;
814 0 : if (pp == 2)
815 0 : z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
816 : else
817 0 : z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
818 0 : return gc_upto(av, z);
819 : }
820 0 : return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
821 : }
822 :
823 : /*Multiple a column vector by a line vector to make a matrix*/
824 : GEN
825 14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
826 : {
827 14 : long i,j, lx=lg(x), ly=lg(y);
828 : GEN z;
829 14 : if (ly==1) return cgetg(1,t_MAT);
830 14 : z = cgetg(ly, t_MAT);
831 49 : for (j=1; j < ly; j++)
832 : {
833 35 : GEN zj = cgetg(lx,t_VECSMALL);
834 112 : for (i=1; i<lx; i++)
835 77 : uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
836 35 : gel(z,j) = zj;
837 : }
838 14 : return z;
839 : }
840 :
841 : /*Multiple a column vector by a line vector to make a matrix*/
842 : GEN
843 5651 : FpC_FpV_mul(GEN x, GEN y, GEN p)
844 : {
845 5651 : long i,j, lx=lg(x), ly=lg(y);
846 : GEN z;
847 5651 : if (ly==1) return cgetg(1,t_MAT);
848 5651 : z = cgetg(ly,t_MAT);
849 63794 : for (j=1; j < ly; j++)
850 : {
851 58143 : GEN zj = cgetg(lx,t_COL);
852 1379166 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
853 58143 : gel(z, j) = zj;
854 : }
855 5651 : return z;
856 : }
857 :
858 : /* Multiply a line vector by a column and return a scalar (t_INT) */
859 : GEN
860 8509289 : FpV_dotproduct(GEN x, GEN y, GEN p)
861 : {
862 8509289 : long i, lx = lg(x);
863 : pari_sp av;
864 : GEN c;
865 8509289 : if (lx == 1) return gen_0;
866 8506783 : av = avma; c = mulii(gel(x,1),gel(y,1));
867 36654878 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
868 8506744 : return gc_INT(av, modii(c,p));
869 : }
870 : GEN
871 0 : FpV_dotsquare(GEN x, GEN p)
872 : {
873 0 : long i, lx = lg(x);
874 : pari_sp av;
875 : GEN c;
876 0 : if (lx == 1) return gen_0;
877 0 : av = avma; c = sqri(gel(x,1));
878 0 : for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
879 0 : return gc_INT(av, modii(c,p));
880 : }
881 :
882 : INLINE ulong
883 9394755 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
884 : {
885 9394755 : ulong c = uel(x,0) * uel(y,0);
886 : long k;
887 72129575 : for (k = 1; k < lx; k++) {
888 62734820 : c += uel(x,k) * uel(y,k);
889 62734820 : if (c & HIGHBIT) c %= p;
890 : }
891 9394755 : return c % p;
892 : }
893 :
894 : INLINE ulong
895 741308 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
896 : {
897 : ulong l0, l1, v1, h0, h1;
898 741308 : long i = 0;
899 : LOCAL_OVERFLOW;
900 : LOCAL_HIREMAINDER;
901 741308 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
902 15418849 : while (++i < lx) {
903 14677541 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
904 14677541 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
905 : }
906 741308 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
907 464137 : else return remlll_pre(v1, h1, l1, p, pi);
908 : }
909 :
910 : ulong
911 730736 : Flv_dotproduct(GEN x, GEN y, ulong p)
912 : {
913 730736 : long lx = lg(x)-1;
914 730736 : if (lx == 0) return 0;
915 730736 : if (SMALL_ULONG(p))
916 730736 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
917 : else
918 0 : return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
919 : }
920 :
921 : ulong
922 60008 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
923 : {
924 60008 : long lx = lg(x)-1;
925 60008 : if (lx == 0) return 0;
926 60008 : if (SMALL_ULONG(p))
927 54348 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
928 : else
929 5660 : return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
930 : }
931 :
932 : ulong
933 9412774 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
934 : {
935 9412774 : long lx = minss(lgpol(x), lgpol(y));
936 9412832 : if (lx == 0) return 0;
937 9345326 : if (pi)
938 735648 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
939 : else
940 8609678 : return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
941 : }
942 : ulong
943 0 : Flx_dotproduct(GEN x, GEN y, ulong p)
944 0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
945 :
946 : GEN
947 1754821 : FpM_FpC_mul(GEN x, GEN y, GEN p)
948 : {
949 1754821 : long lx = lg(x);
950 1754821 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
951 : }
952 : GEN
953 1077194 : Flm_Flc_mul(GEN x, GEN y, ulong p)
954 : {
955 1077194 : long l, lx = lg(x);
956 1077194 : if (lx==1) return cgetg(1,t_VECSMALL);
957 1077194 : l = lgcols(x);
958 1077194 : if (p==2)
959 329856 : return Flm_Flc_mul_i_2(x, y, lx, l);
960 747338 : else if (SMALL_ULONG(p))
961 510000 : return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
962 : else
963 237338 : return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
964 : }
965 :
966 : /* allow pi = 0 */
967 : GEN
968 6789 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
969 : {
970 6789 : long l, lx = lg(x);
971 : GEN z;
972 6789 : if (lx==1) return cgetg(1,t_VECSMALL);
973 6789 : l = lgcols(x);
974 6789 : z = cgetg(l, t_VECSMALL);
975 6789 : if (SMALL_ULONG(p))
976 4405 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
977 : else
978 2384 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
979 6790 : return z;
980 : }
981 :
982 : /* allow pi = 0 */
983 : GEN
984 7541394 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
985 : {
986 7541394 : long l, lx = lg(x);
987 : GEN z;
988 7541394 : if (lx==1) return pol0_Flx(sv);
989 7541394 : l = lgcols(x);
990 7540034 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
991 7541981 : if (SMALL_ULONG(p))
992 3286448 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
993 : else
994 4255533 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
995 7551519 : return Flx_renormalize(z, l + 1);
996 : }
997 :
998 : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
999 : GEN
1000 1234390 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
1001 : {
1002 1234390 : long i, l, lx = lg(x);
1003 : GEN z;
1004 1234390 : if (lx==1) return pol_0(v);
1005 1234390 : l = lgcols(x);
1006 1234390 : z = new_chunk(l+1);
1007 1688748 : for (i=l-1; i; i--)
1008 : {
1009 1631762 : pari_sp av = avma;
1010 1631762 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1011 1631760 : p1 = modii(p1, p);
1012 1631759 : if (signe(p1))
1013 : {
1014 1177401 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
1015 1177401 : gel(z,i+1) = gc_INT(av, p1);
1016 1177401 : break;
1017 : }
1018 454358 : set_avma(av);
1019 : }
1020 1234387 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1021 1177401 : z[0] = evaltyp(t_POL) | _evallg(i+2);
1022 1177401 : z[1] = evalsigne(1) | evalvarn(v);
1023 3437149 : for (; i; i--)
1024 : {
1025 2259745 : pari_sp av = avma;
1026 2259745 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1027 2259744 : gel(z,i+1) = gc_INT(av, modii(p1,p));
1028 : }
1029 1177404 : return z;
1030 : }
1031 :
1032 : /********************************************************************/
1033 : /** **/
1034 : /** TRANSPOSITION **/
1035 : /** **/
1036 : /********************************************************************/
1037 :
1038 : /* == zm_transpose */
1039 : GEN
1040 735386 : Flm_transpose(GEN x)
1041 : {
1042 735386 : long i, dx, lx = lg(x);
1043 : GEN y;
1044 735386 : if (lx == 1) return cgetg(1,t_MAT);
1045 735253 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1046 9445745 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1047 735256 : return y;
1048 : }
1049 :
1050 : /********************************************************************/
1051 : /** **/
1052 : /** SCALAR MATRICES **/
1053 : /** **/
1054 : /********************************************************************/
1055 :
1056 : GEN
1057 2464 : gen_matid(long n, void *E, const struct bb_field *S)
1058 : {
1059 2464 : GEN y = cgetg(n+1,t_MAT), _0, _1;
1060 : long i;
1061 2464 : if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
1062 2464 : _0 = S->s(E,0);
1063 2464 : _1 = S->s(E,1);
1064 11515 : for (i=1; i<=n; i++)
1065 : {
1066 9051 : GEN z = const_col(n, _0); gel(z,i) = _1;
1067 9051 : gel(y, i) = z;
1068 : }
1069 2464 : return y;
1070 : }
1071 :
1072 : GEN
1073 35 : matid_F2xqM(long n, GEN T)
1074 : {
1075 : void *E;
1076 35 : const struct bb_field *S = get_F2xq_field(&E, T);
1077 35 : return gen_matid(n, E, S);
1078 : }
1079 : GEN
1080 56 : matid_FlxqM(long n, GEN T, ulong p)
1081 : {
1082 : void *E;
1083 56 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1084 56 : return gen_matid(n, E, S);
1085 : }
1086 :
1087 : GEN
1088 1435456 : matid_Flm(long n)
1089 : {
1090 1435456 : GEN y = cgetg(n+1,t_MAT);
1091 : long i;
1092 1435457 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1093 9110288 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1094 1435454 : return y;
1095 : }
1096 :
1097 : GEN
1098 42 : scalar_Flm(long s, long n)
1099 : {
1100 : long i;
1101 42 : GEN y = cgetg(n+1,t_MAT);
1102 560 : for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
1103 42 : return y;
1104 : }
1105 :
1106 : /********************************************************************/
1107 : /** **/
1108 : /** CONVERSIONS **/
1109 : /** **/
1110 : /********************************************************************/
1111 : GEN
1112 56245293 : ZV_to_Flv(GEN x, ulong p)
1113 709187827 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1114 :
1115 : GEN
1116 12899850 : ZM_to_Flm(GEN x, ulong p)
1117 68041387 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1118 :
1119 : GEN
1120 14722 : ZVV_to_FlvV(GEN x, ulong p)
1121 112431 : { pari_APPLY_type(t_VEC, ZV_to_Flv(gel(x,i), p)) }
1122 :
1123 : GEN
1124 4940 : ZMV_to_FlmV(GEN x, ulong m)
1125 41848 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1126 :
1127 : /* TO INTMOD */
1128 : static GEN
1129 18978230 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1130 : static GEN
1131 578850 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1132 :
1133 : GEN
1134 4013 : Fp_to_mod(GEN z, GEN p)
1135 : {
1136 4013 : retmkintmod(modii(z, p), icopy(p));
1137 : }
1138 :
1139 : GEN
1140 997241 : FpX_to_mod_raw(GEN z, GEN p)
1141 : {
1142 997241 : long i, l = lg(z);
1143 : GEN x;
1144 997241 : if (l == 2)
1145 : {
1146 522424 : x = cgetg(3,t_POL); x[1] = z[1];
1147 522424 : gel(x,2) = mkintmod(gen_0,p); return x;
1148 : }
1149 474817 : x = cgetg(l,t_POL); x[1] = z[1];
1150 3918607 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1151 474817 : return normalizepol_lg(x,l);
1152 : }
1153 :
1154 : /* z in Z[X], return z * Mod(1,p), normalized*/
1155 : GEN
1156 1260819 : FpX_to_mod(GEN z, GEN p)
1157 : {
1158 1260819 : long i, l = lg(z);
1159 : GEN x;
1160 1260819 : if (l == 2)
1161 : {
1162 1197 : x = cgetg(3,t_POL); x[1] = z[1];
1163 1197 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1164 : }
1165 1259622 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1166 16690366 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1167 1259622 : return normalizepol_lg(x,l);
1168 : }
1169 :
1170 : GEN
1171 84021 : FpXC_to_mod(GEN z, GEN p)
1172 : {
1173 84021 : long i,l = lg(z);
1174 84021 : GEN x = cgetg(l, t_COL);
1175 84021 : p = icopy(p);
1176 474159 : for (i=1; i<l; i++)
1177 390138 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1178 84021 : return x;
1179 : }
1180 :
1181 : static GEN
1182 0 : FpXC_to_mod_raw(GEN x, GEN p)
1183 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1184 :
1185 : GEN
1186 0 : FpXM_to_mod(GEN z, GEN p)
1187 : {
1188 0 : long i,l = lg(z);
1189 0 : GEN x = cgetg(l, t_MAT);
1190 0 : p = icopy(p);
1191 0 : for (i=1; i<l; i++)
1192 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1193 0 : return x;
1194 : }
1195 :
1196 : /* z in Z^n, return z * Mod(1,p), normalized*/
1197 : GEN
1198 36495 : FpV_to_mod(GEN z, GEN p)
1199 : {
1200 36495 : long i,l = lg(z);
1201 36495 : GEN x = cgetg(l, t_VEC);
1202 36495 : if (l == 1) return x;
1203 36495 : p = icopy(p);
1204 73585 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1205 36495 : return x;
1206 : }
1207 : /* z in Z^n, return z * Mod(1,p), normalized*/
1208 : GEN
1209 105 : FpC_to_mod(GEN z, GEN p)
1210 : {
1211 105 : long i, l = lg(z);
1212 105 : GEN x = cgetg(l, t_COL);
1213 105 : if (l == 1) return x;
1214 91 : p = icopy(p);
1215 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1216 91 : return x;
1217 : }
1218 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1219 : GEN
1220 226 : FpM_to_mod(GEN z, GEN p)
1221 : {
1222 226 : long i, j, m, l = lg(z);
1223 226 : GEN x = cgetg(l,t_MAT), y, zi;
1224 226 : if (l == 1) return x;
1225 205 : m = lgcols(z);
1226 205 : p = icopy(p);
1227 2456 : for (i=1; i<l; i++)
1228 : {
1229 2251 : gel(x,i) = cgetg(m,t_COL);
1230 2251 : y = gel(x,i); zi= gel(z,i);
1231 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1232 : }
1233 205 : return x;
1234 : }
1235 : GEN
1236 28 : Flc_to_mod(GEN z, ulong pp)
1237 : {
1238 28 : long i, l = lg(z);
1239 28 : GEN p, x = cgetg(l, t_COL);
1240 28 : if (l == 1) return x;
1241 28 : p = utoipos(pp);
1242 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1243 28 : return x;
1244 : }
1245 : GEN
1246 59694 : Flm_to_mod(GEN z, ulong pp)
1247 : {
1248 59694 : long i, j, m, l = lg(z);
1249 59694 : GEN p, x = cgetg(l,t_MAT), y, zi;
1250 59694 : if (l == 1) return x;
1251 59666 : m = lgcols(z);
1252 59666 : p = utoipos(pp);
1253 201328 : for (i=1; i<l; i++)
1254 : {
1255 141662 : gel(x,i) = cgetg(m,t_COL);
1256 141662 : y = gel(x,i); zi= gel(z,i);
1257 720428 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1258 : }
1259 59666 : return x;
1260 : }
1261 :
1262 : GEN
1263 574 : FpVV_to_mod(GEN z, GEN p)
1264 : {
1265 574 : long i, j, m, l = lg(z);
1266 574 : GEN x = cgetg(l,t_VEC), y, zi;
1267 574 : if (l == 1) return x;
1268 574 : m = lgcols(z);
1269 574 : p = icopy(p);
1270 1246 : for (i=1; i<l; i++)
1271 : {
1272 672 : gel(x,i) = cgetg(m,t_VEC);
1273 672 : y = gel(x,i); zi= gel(z,i);
1274 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1275 : }
1276 574 : return x;
1277 : }
1278 :
1279 : /* z in Z^n, return z * Mod(1,p), normalized*/
1280 : GEN
1281 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1282 : {
1283 7 : long i,l = lg(z);
1284 7 : GEN x = cgetg(l, t_COL);
1285 7 : if (l == 1) return x;
1286 7 : p = icopy(p);
1287 7 : T = FpX_to_mod_raw(T, p);
1288 21 : for (i=1; i<l; i++)
1289 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1290 7 : return x;
1291 : }
1292 :
1293 : static GEN
1294 582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1295 : {
1296 582533 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1297 582533 : return mkpolmod(z, T);
1298 : }
1299 :
1300 : /* z in Z^n, return z * Mod(1,p), normalized*/
1301 : GEN
1302 28 : FqC_to_mod(GEN z, GEN T, GEN p)
1303 : {
1304 : GEN x;
1305 28 : long i,l = lg(z);
1306 28 : if (!T) return FpC_to_mod(z, p);
1307 28 : x = cgetg(l, t_COL);
1308 28 : if (l == 1) return x;
1309 28 : p = icopy(p);
1310 28 : T = FpX_to_mod_raw(T, p);
1311 504 : for (i=1; i<l; i++)
1312 476 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1313 28 : return x;
1314 : }
1315 :
1316 : /* z in Z^n, return z * Mod(1,p), normalized*/
1317 : static GEN
1318 107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1319 690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1320 :
1321 : /* z in Z^n, return z * Mod(1,p), normalized*/
1322 : GEN
1323 26145 : FqM_to_mod(GEN z, GEN T, GEN p)
1324 : {
1325 : GEN x;
1326 26145 : long i,l = lg(z);
1327 26145 : if (!T) return FpM_to_mod(z, p);
1328 26145 : x = cgetg(l, t_MAT);
1329 26145 : if (l == 1) return x;
1330 26117 : p = icopy(p);
1331 26117 : T = FpX_to_mod_raw(T, p);
1332 134092 : for (i=1; i<l; i++)
1333 107975 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1334 26117 : return x;
1335 : }
1336 :
1337 : /********************************************************************/
1338 : /* */
1339 : /* Blackbox linear algebra */
1340 : /* */
1341 : /********************************************************************/
1342 :
1343 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1344 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1345 : * (e_j) is the canonical basis.
1346 : * A sparse matrix (zMs) is a t_VEC of zCs */
1347 :
1348 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1349 : * integer representing an element of Fp. This is important since p can be
1350 : * large and p+E[i] would not fit in a C long. */
1351 :
1352 : /* RgCs and RgMs are similar, except that the type of the component is
1353 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1354 : * of E. */
1355 :
1356 : /* Most functions take an argument nbrow which is the number of lines of the
1357 : * column/matrix, which cannot be derived from the data. */
1358 :
1359 : GEN
1360 0 : zCs_to_ZC(GEN R, long nbrow)
1361 : {
1362 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1363 0 : long j, l = lg(C);
1364 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1365 0 : return c;
1366 : }
1367 :
1368 : GEN
1369 0 : zMs_to_ZM(GEN x, long nbrow)
1370 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1371 :
1372 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1373 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1374 : GEN
1375 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1376 : {
1377 147 : pari_sp ltop = avma;
1378 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1379 147 : if (ZV_equal0(B)) return zerocol(n);
1380 147 : while (++col <= n)
1381 : {
1382 147 : pari_sp btop = avma, av;
1383 : long i, lQ;
1384 147 : GEN V, Q, M, W = B;
1385 147 : GEN b = cgetg(m+2, t_POL);
1386 147 : b[1] = evalsigne(1)|evalvarn(0);
1387 147 : gel(b, 2) = gel(W, col);
1388 46219 : for (i = 3; i<m+2; i++)
1389 46072 : gel(b, i) = cgeti(lgefint(p));
1390 147 : av = avma;
1391 46219 : for (i = 3; i<m+2; i++)
1392 : {
1393 46072 : W = f(E, W);
1394 46072 : affii(gel(W, col),gel(b, i));
1395 46072 : if (gc_needed(av,1))
1396 : {
1397 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1398 72 : W = gc_upto(av, W);
1399 : }
1400 : }
1401 147 : b = FpX_renormalize(b, m+2);
1402 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1403 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1404 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1405 147 : W = B; lQ =lg(Q);
1406 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1407 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1408 147 : av = avma;
1409 22742 : for (i = lQ-3; i > 1; i--)
1410 : {
1411 22595 : W = f(E, W);
1412 22595 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1413 22595 : if (gc_needed(av,1))
1414 : {
1415 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1416 110 : (void)gc_all(av, 2, &V, &W);
1417 : }
1418 : }
1419 147 : V = FpC_red(V, p);
1420 147 : W = FpC_sub(f(E,V), B, p);
1421 147 : if (ZV_equal0(W)) return gc_GEN(ltop, V);
1422 0 : av = avma;
1423 0 : for (i = 1; i <= n; ++i)
1424 : {
1425 0 : V = W;
1426 0 : W = f(E, V);
1427 0 : if (ZV_equal0(W))
1428 0 : return gc_GEN(ltop, shallowtrans(V));
1429 0 : (void)gc_all(av, 2, &V, &W);
1430 : }
1431 0 : set_avma(btop);
1432 : }
1433 0 : return NULL;
1434 : }
1435 :
1436 : GEN
1437 0 : zMs_ZC_mul(GEN M, GEN B)
1438 : {
1439 : long i, j;
1440 0 : long n = lg(B)-1;
1441 0 : GEN V = zerocol(n);
1442 0 : for (i = 1; i <= n; ++i)
1443 0 : if (signe(gel(B, i)))
1444 : {
1445 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1446 0 : long l = lg(C);
1447 0 : for (j = 1; j < l; ++j)
1448 : {
1449 0 : long k = C[j];
1450 0 : switch(E[j])
1451 : {
1452 0 : case 1:
1453 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1454 0 : break;
1455 0 : case -1:
1456 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1457 0 : break;
1458 0 : default:
1459 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]));
1460 0 : break;
1461 : }
1462 : }
1463 : }
1464 0 : return V;
1465 : }
1466 :
1467 : GEN
1468 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1469 :
1470 : GEN
1471 68961 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1472 : {
1473 68961 : long i, j, lM = lg(M);
1474 68961 : GEN V = cgetg(lM,t_VEC);
1475 28332297 : for (i = 1; i < lM; ++i)
1476 : {
1477 28263336 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1478 28263336 : pari_sp av = avma;
1479 28263336 : long lC = lg(C);
1480 28263336 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1481 28263336 : z = mulis(gel(B, C[1]), E[1]);
1482 228730820 : for (j = 2; j < lC; ++j)
1483 : {
1484 200467484 : GEN b = gel(B, C[j]);
1485 200467484 : switch(E[j])
1486 : {
1487 131566519 : case 1: z = addii(z, b); break;
1488 57525441 : case -1: z = subii(z, b); break;
1489 11375524 : default: z = addii(z, mulis(b, E[j])); break;
1490 : }
1491 : }
1492 28263336 : gel(V,i) = gc_INT(av, p? Fp_red(z, p): z);
1493 : }
1494 68961 : return V;
1495 : }
1496 : GEN
1497 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1498 :
1499 : GEN
1500 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1501 : {
1502 0 : pari_sp av = avma, av2;
1503 0 : GEN xi, xb, pi = gen_1, P;
1504 : long i;
1505 0 : if (!C) {
1506 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1507 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1508 : }
1509 0 : P = utoipos(p);
1510 0 : av2 = avma;
1511 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1512 0 : xb = Flm_to_ZM(xi);
1513 0 : for (i = 2; i <= e; i++)
1514 : {
1515 0 : pi = muliu(pi, p); /* = p^(i-1) */
1516 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1517 0 : if (gc_needed(av,2))
1518 : {
1519 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1520 0 : (void)gc_all(av2,3, &pi,&b,&xb);
1521 : }
1522 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1523 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1524 : }
1525 0 : return gc_upto(av, xb);
1526 : }
1527 :
1528 : struct wrapper_modp_s {
1529 : GEN (*f)(void*, GEN);
1530 : void *E;
1531 : GEN p;
1532 : };
1533 :
1534 : /* compute f(x) mod p */
1535 : static GEN
1536 0 : wrap_relcomb_modp(void *E, GEN x)
1537 : {
1538 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1539 0 : return FpC_red(W->f(W->E, x), W->p);
1540 : }
1541 : static GEN
1542 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1543 :
1544 : static GEN
1545 68814 : wrap_relker(void*E, GEN x)
1546 : {
1547 68814 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1548 68814 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1549 : }
1550 :
1551 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1552 : GEN
1553 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1554 : {
1555 : struct wrapper_modp_s W;
1556 0 : pari_sp av = avma;
1557 0 : GEN xb, xi, pi = gen_1;
1558 : long i;
1559 0 : W.E = E;
1560 0 : W.f = f;
1561 0 : W.p = p;
1562 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1563 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1564 0 : xb = xi;
1565 0 : for (i = 2; i <= e; i++)
1566 : {
1567 0 : pi = mulii(pi, p); /* = p^(i-1) */
1568 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1569 0 : if (gc_needed(av,2))
1570 : {
1571 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1572 0 : (void)gc_all(av,3, &pi,&B,&xb);
1573 : }
1574 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1575 0 : if (!xi) return NULL;
1576 0 : if (typ(xi) == t_VEC) return gc_upto(av, xi);
1577 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1578 : }
1579 0 : return gc_upto(av, xb);
1580 : }
1581 :
1582 : static GEN
1583 23036 : vecprow(GEN A, GEN prow)
1584 : {
1585 23036 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1586 : }
1587 :
1588 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1589 : * or the index of a column which is linearly dependent from the others as a
1590 : * t_VECSMALL with a single component. */
1591 : GEN
1592 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1593 : {
1594 0 : pari_sp av = avma;
1595 : GEN pcol, prow;
1596 0 : long nbi=lg(M)-1, lR;
1597 : long i, n;
1598 : GEN Mp, Ap, Rp;
1599 : pari_timer ti;
1600 0 : if (DEBUGLEVEL) timer_start(&ti);
1601 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1602 0 : if (!pcol) return gc_NULL(av);
1603 0 : if (DEBUGLEVEL)
1604 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1605 0 : n = lg(pcol)-1;
1606 0 : Mp = cgetg(n+1, t_MAT);
1607 0 : for(i=1; i<=n; i++)
1608 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1609 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1610 0 : if (DEBUGLEVEL) timer_start(&ti);
1611 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1612 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1613 0 : if (!Rp) return gc_NULL(av);
1614 0 : lR = lg(Rp)-1;
1615 0 : if (typ(Rp) == t_COL)
1616 : {
1617 0 : GEN R = zerocol(nbi+1);
1618 0 : for (i=1; i<=lR; i++)
1619 0 : gel(R,pcol[i]) = gel(Rp,i);
1620 0 : return gc_GEN(av, R);
1621 : }
1622 0 : for (i = 1; i <= lR; ++i)
1623 0 : if (signe(gel(Rp, i)))
1624 0 : return gc_leaf(av, mkvecsmall(pcol[i]));
1625 : return NULL; /* LCOV_EXCL_LINE */
1626 : }
1627 :
1628 : GEN
1629 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1630 : {
1631 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1632 : }
1633 :
1634 : GEN
1635 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1636 : {
1637 : GEN res;
1638 0 : pari_CATCH(e_INV)
1639 : {
1640 0 : GEN E = pari_err_last();
1641 0 : GEN x = err_get_compo(E,2);
1642 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1643 0 : if (DEBUGLEVEL)
1644 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1645 0 : res = NULL;
1646 : } pari_TRY {
1647 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1648 0 : } pari_ENDCATCH
1649 0 : return res;
1650 : }
1651 :
1652 : static GEN
1653 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1654 : {
1655 147 : long i, j, oldf = 0, f = 0;
1656 147 : long lrow = lg(prow), lM = lg(M);
1657 147 : GEN W = const_vecsmall(lM-1,1);
1658 147 : GEN R = cgetg(lrow, t_VEC);
1659 228354 : for (i=1; i<lrow; i++)
1660 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1661 : do
1662 : {
1663 442 : oldf = f;
1664 374995 : for(i=1; i<lM; i++)
1665 374553 : if (W[i])
1666 : {
1667 131733 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1668 131733 : long c=0, cj=0, lF = lg(F);
1669 1093676 : for(j=1; j<lF; j++)
1670 961943 : if (!gel(R,F[j])) { c++; cj=j; }
1671 131733 : if (c>=2) continue;
1672 112310 : if (c==1)
1673 : {
1674 32898 : pari_sp av = avma;
1675 32898 : GEN s = gen_0;
1676 272959 : for(j=1; j<lF; j++)
1677 240061 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1678 : /* s /= -E[cj] mod p */
1679 32898 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1680 32898 : gel(R,F[cj]) = gc_INT(av, s);
1681 32898 : f++;
1682 : }
1683 112310 : W[i]=0;
1684 : }
1685 442 : } while(oldf!=f);
1686 228354 : for (i=1; i<lrow; i++)
1687 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1688 147 : return R;
1689 : }
1690 :
1691 : /* Return a linear form Y such that YM is essentially 0 */
1692 : GEN
1693 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1694 : {
1695 147 : pari_sp av = avma, av2;
1696 : GEN Mp, pcol, prow;
1697 : long i, n;
1698 : pari_timer ti;
1699 : struct wrapper_modp_s W;
1700 147 : if (DEBUGLEVEL) timer_start(&ti);
1701 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1702 147 : if (!pcol) return gc_NULL(av);
1703 147 : if (DEBUGLEVEL)
1704 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1705 147 : n = lg(pcol)-1;
1706 147 : Mp = cgetg(n+1, t_MAT);
1707 23183 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1708 147 : W.E = (void*)Mp;
1709 147 : W.p = p;
1710 147 : av2 = avma;
1711 0 : for(;; set_avma(av2))
1712 0 : {
1713 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1714 147 : if (DEBUGLEVEL) timer_start(&ti);
1715 147 : pari_CATCH(e_INV)
1716 : {
1717 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1718 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1719 0 : if (DEBUGLEVEL)
1720 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1721 0 : Rp = NULL;
1722 : } pari_TRY {
1723 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1724 147 : } pari_ENDCATCH
1725 147 : if (!Rp) continue;
1726 147 : if (typ(Rp)==t_COL)
1727 : {
1728 147 : Rp = FpC_sub(Rp,B,p);
1729 147 : if (ZV_equal0(Rp)) continue;
1730 : }
1731 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1732 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1733 147 : return gc_GEN(av, R);
1734 : }
1735 : }
1736 :
1737 : GEN
1738 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1739 : {
1740 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1741 : }
|
