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 29387296 : FpC_red(GEN x, GEN p)
28 252210552 : { 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 6865159 : FpC_center(GEN x, GEN p, GEN pov2)
37 43234208 : { 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 5156818 : FpM_red(GEN x, GEN p)
46 23836581 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2756195 : FpM_center(GEN x, GEN p, GEN pov2)
50 8300464 : { 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 23181 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
98 147 : return y;
99 : }
100 :
101 : GEN
102 924 : random_FpC(long d, GEN p)
103 : {
104 : long i;
105 924 : GEN y = cgetg(d+1,t_COL);
106 6188 : for (i=1; i<=d; i++) gel(y,i) = randomi(p);
107 924 : return y;
108 : }
109 :
110 : GEN
111 41334 : random_Flv(long n, ulong p)
112 : {
113 41334 : GEN y = cgetg(n+1, t_VECSMALL);
114 : long i;
115 261418 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
116 41334 : return y;
117 : }
118 :
119 : /********************************************************************/
120 : /** **/
121 : /** ADD, SUB **/
122 : /** **/
123 : /********************************************************************/
124 : GEN
125 279817 : FpC_add(GEN x, GEN y, GEN p)
126 3620767 : { 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 1663859 : Flv_add(GEN x, GEN y, ulong p)
138 : {
139 1663859 : if (p==2)
140 6520848 : pari_APPLY_ulong(x[i]^y[i])
141 : else
142 38489403 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
143 : }
144 :
145 : void
146 441884 : Flv_add_inplace(GEN x, GEN y, ulong p)
147 : {
148 441884 : long i, l = lg(x);
149 441884 : if (p==2)
150 1373290 : 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 441884 : }
154 :
155 : ulong
156 3871 : Flv_sum(GEN x, ulong p)
157 : {
158 3871 : long i, l = lg(x);
159 3871 : ulong s = 0;
160 3871 : if (p==2)
161 17689 : 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 3871 : return s;
165 : }
166 :
167 : GEN
168 932749 : FpC_sub(GEN x, GEN y, GEN p)
169 18885838 : { 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 227952732 : Flv_sub(GEN x, GEN y, ulong p)
181 1024006155 : { 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 51129 : Flm_Fl_add(GEN x, ulong y, ulong p)
192 : {
193 51129 : long l = lg(x), i, j;
194 51129 : GEN z = cgetg(l,t_MAT);
195 :
196 51129 : if (l==1) return z;
197 51129 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
198 240844 : for (i=1; i<l; i++)
199 : {
200 189715 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
201 189715 : gel(z,i) = zi;
202 1318352 : for (j=1; j<l; j++) zi[j] = xi[j];
203 189715 : zi[i] = Fl_add(zi[i], y, p);
204 : }
205 51129 : return z;
206 : }
207 : GEN
208 3745 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
209 :
210 : GEN
211 25350 : Flm_add(GEN x, GEN y, ulong p)
212 618167 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
213 :
214 : GEN
215 25481308 : Flm_sub(GEN x, GEN y, ulong p)
216 253417955 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
217 :
218 : /********************************************************************/
219 : /** **/
220 : /** MULTIPLICATION **/
221 : /** **/
222 : /********************************************************************/
223 : GEN
224 948116 : FpC_Fp_mul(GEN x, GEN y, GEN p)
225 11733238 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
226 :
227 : GEN
228 1579384 : Flv_Fl_mul(GEN x, ulong y, ulong p)
229 37251172 : { pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
230 :
231 : GEN
232 4692 : Flv_Fl_div(GEN x, ulong y, ulong p)
233 : {
234 4692 : 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 762055 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
243 762055 : long i, j, h, l = lg(X);
244 762055 : GEN A = cgetg(l, t_MAT);
245 762056 : if (l == 1) return A;
246 762056 : h = lgcols(X);
247 4039508 : for (j=1; j<l; j++)
248 : {
249 3277563 : GEN a = cgetg(h, t_COL), x = gel(X, j);
250 23869531 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
251 3277451 : gel(A,j) = a;
252 : }
253 761945 : return A;
254 : }
255 :
256 : /* x *= y */
257 : void
258 7733204 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
259 : {
260 : long i;
261 7733204 : if (HIGHWORD(y | p))
262 8555930 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
263 : else
264 31954696 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
265 7733203 : }
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 16981150 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
288 : {
289 16981150 : long i, j, m, l = lg(y);
290 16981150 : GEN z = cgetg(l, t_MAT);
291 16979118 : if (l == 1) return z;
292 16979118 : m = lgcols(y);
293 152000442 : for(j=1; j<l; j++) {
294 134947079 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
295 362828502 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
296 : }
297 17053363 : return z;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 8517335 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
303 : {
304 8517335 : long i, j, m, l = lg(y);
305 8517335 : GEN z = cgetg(l, t_MAT);
306 8517335 : if (l == 1) return z;
307 8517335 : m = lgcols(y);
308 52657465 : for(j=1; j<l; j++) {
309 44140173 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 157402218 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
311 : }
312 8517292 : return z;
313 : }
314 :
315 : /* return x * y */
316 : GEN
317 25429768 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
318 : {
319 25429768 : if (HIGHWORD(p))
320 16981172 : return Flm_Fl_mul_pre_i(y, x, p, pi);
321 : else
322 8448596 : return Flm_Fl_mul_OK(y, x, p);
323 : }
324 :
325 : /* return x * y */
326 : GEN
327 68768 : Flm_Fl_mul(GEN y, ulong x, ulong p)
328 : {
329 68768 : if (HIGHWORD(p))
330 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
331 : else
332 68694 : return Flm_Fl_mul_OK(y, x, p);
333 : }
334 :
335 : GEN
336 4064703 : Flv_neg(GEN x, ulong p)
337 31146870 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
338 :
339 : void
340 6184 : Flv_neg_inplace(GEN v, ulong p)
341 : {
342 : long i;
343 186524 : for (i = 1; i < lg(v); ++i)
344 180340 : v[i] = Fl_neg(v[i], p);
345 6184 : }
346 :
347 : GEN
348 320494 : Flm_neg(GEN x, ulong p)
349 4354399 : { 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 19390632 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
354 : {
355 19390632 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
356 : long k;
357 233107375 : for (k = 2; k < lx; k++)
358 : {
359 213718461 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
360 213715557 : if (signe(t)) c = addii(c, t);
361 : }
362 19388914 : return c;
363 : }
364 :
365 : static long
366 97688493 : zmrow_zc_mul(GEN x, GEN y, long lx, long i)
367 : {
368 : long k;
369 97688493 : long c = coeff(x,i,1) * y[1];
370 1498281456 : for (k = 2; k < lx; k++)
371 1400592963 : c += coeff(x,i,k) * y[k];
372 97688493 : return c;
373 : }
374 :
375 : GEN
376 6450479 : zm_zc_mul(GEN x, GEN y)
377 : {
378 6450479 : long lx = lg(x), l, i;
379 : GEN z;
380 6450479 : if (lx == 1) return cgetg(1, t_VECSMALL);
381 6450479 : l = lg(gel(x,1));
382 6450479 : z = cgetg(l,t_VECSMALL);
383 104138972 : for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
384 6450479 : 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 759491183 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
406 : {
407 759491183 : ulong c = ucoeff(x,i,1) * uel(y,1);
408 : long k;
409 11616448114 : for (k = 2; k < lx; k++) {
410 10856956931 : c += ucoeff(x,i,k) * uel(y,k);
411 10856956931 : if (c & HIGHBIT) c %= p;
412 : }
413 759491183 : return c % p;
414 : }
415 :
416 : static ulong
417 531365250 : 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 531365250 : long k = 1;
421 : LOCAL_OVERFLOW;
422 : LOCAL_HIREMAINDER;
423 531365250 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
424 8952388011 : while (++k < lx) {
425 8421022761 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
426 8421022761 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
427 : }
428 531365250 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
429 54537949 : else return remlll_pre(v1, h1, l1, p, pi);
430 : }
431 :
432 : static GEN
433 276443 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
434 : {
435 : long i,j;
436 276443 : GEN z = NULL;
437 :
438 3225730 : for (j=1; j<lx; j++)
439 : {
440 2949287 : if (!y[j]) continue;
441 1040152 : if (!z) z = Flv_copy(gel(x,j));
442 29407906 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
443 : }
444 276443 : if (!z) z = zero_zv(l-1);
445 276443 : return z;
446 : }
447 :
448 : static GEN
449 1737467 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
450 : {
451 1737467 : GEN z = cgetg(l,t_COL);
452 : long i;
453 16782559 : for (i = 1; i < l; i++)
454 : {
455 15045143 : pari_sp av = avma;
456 15045143 : GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
457 15044916 : gel(z,i) = gerepileuptoint(av, modii(c,p));
458 : }
459 1737416 : return z;
460 : }
461 :
462 : static void
463 85922243 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
464 : {
465 : long i;
466 845431665 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
467 85944891 : }
468 : static GEN
469 82613815 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
470 : {
471 82613815 : GEN z = cgetg(l,t_VECSMALL);
472 82612496 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
473 82614097 : return z;
474 : }
475 :
476 : static void
477 94027586 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
478 : {
479 : long i;
480 628136890 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
481 94495529 : }
482 : static GEN
483 89903729 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
484 : {
485 89903729 : GEN z = cgetg(l,t_VECSMALL);
486 89788468 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
487 89953103 : return z;
488 : }
489 :
490 : GEN
491 1545287 : FpM_mul(GEN x, GEN y, GEN p)
492 : {
493 1545287 : long lx=lg(x), ly=lg(y);
494 : GEN z;
495 1545287 : pari_sp av = avma;
496 1545287 : if (ly==1) return cgetg(1,t_MAT);
497 1545287 : if (lx==1) return zeromat(0, ly-1);
498 1545287 : if (lgefint(p) == 3)
499 : {
500 1499391 : ulong pp = uel(p,2);
501 1499391 : if (pp == 2)
502 : {
503 698021 : x = ZM_to_F2m(x);
504 698076 : y = ZM_to_F2m(y);
505 698082 : z = F2m_to_ZM(F2m_mul(x,y));
506 : }
507 : else
508 : {
509 801370 : x = ZM_to_Flm(x, pp);
510 801383 : y = ZM_to_Flm(y, pp);
511 801383 : z = Flm_to_ZM(Flm_mul(x,y, pp));
512 : }
513 : } else
514 45896 : z = FpM_red(ZM_mul(x, y), p);
515 1545363 : return gerepileupto(av, z);
516 : }
517 :
518 : static GEN
519 27994132 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
520 : {
521 : long j;
522 27994132 : GEN z = cgetg(ly, t_MAT);
523 27992919 : if (SMALL_ULONG(p))
524 102233490 : for (j=1; j<ly; j++)
525 82150172 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
526 : else
527 97608503 : for (j=1; j<ly; j++)
528 89698934 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
529 27992887 : return z;
530 : }
531 :
532 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
533 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
534 : static void
535 66664 : add_slices_ip(long m, long n,
536 : GEN A, long ma, long da, long na, long ea,
537 : GEN B, long mb, long db, long nb, long eb,
538 : GEN M, long dy, long dx, ulong p)
539 : {
540 66664 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
541 : GEN C;
542 :
543 2807329 : for (j = 1; j <= min_e; j++) {
544 2740665 : C = gel(M, j + dx) + dy;
545 125602844 : for (i = 1; i <= min_d; i++)
546 122862179 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
547 122861856 : ucoeff(B, mb + i, nb + j), p);
548 2808608 : for (; i <= da; i++)
549 67620 : uel(C, i) = ucoeff(A, ma + i, na + j);
550 2740988 : for (; i <= db; i++)
551 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
552 2740988 : for (; i <= m; i++)
553 0 : uel(C, i) = 0;
554 : }
555 70367 : for (; j <= ea; j++) {
556 3703 : C = gel(M, j + dx) + dy;
557 201369 : for (i = 1; i <= da; i++)
558 197666 : uel(C, i) = ucoeff(A, ma + i, na + j);
559 3703 : for (; i <= m; i++)
560 0 : uel(C, i) = 0;
561 : }
562 66664 : for (; j <= eb; j++) {
563 0 : C = gel(M, j + dx) + dy;
564 0 : for (i = 1; i <= db; i++)
565 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
566 0 : for (; i <= m; i++)
567 0 : uel(C, i) = 0;
568 : }
569 66664 : for (; j <= n; j++) {
570 0 : C = gel(M, j + dx) + dy;
571 0 : for (i = 1; i <= m; i++)
572 0 : uel(C, i) = 0;
573 : }
574 66664 : }
575 :
576 : static GEN
577 33332 : add_slices(long m, long n,
578 : GEN A, long ma, long da, long na, long ea,
579 : GEN B, long mb, long db, long nb, long eb, ulong p)
580 : {
581 : GEN M;
582 : long j;
583 33332 : M = cgetg(n + 1, t_MAT);
584 1428233 : for (j = 1; j <= n; j++)
585 1394905 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
586 33328 : add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
587 33332 : return M;
588 : }
589 :
590 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
591 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
592 : static GEN
593 58330 : subtract_slices(long m, long n,
594 : GEN A, long ma, long da, long na, long ea,
595 : GEN B, long mb, long db, long nb, long eb, ulong p)
596 : {
597 58330 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
598 58330 : GEN M = cgetg(n + 1, t_MAT), C;
599 :
600 2549732 : for (j = 1; j <= min_e; j++) {
601 2491401 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
602 117622833 : for (i = 1; i <= min_d; i++)
603 115131411 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
604 115131433 : ucoeff(B, mb + i, nb + j), p);
605 2594897 : for (; i <= da; i++)
606 103497 : uel(C, i) = ucoeff(A, ma + i, na + j);
607 2569549 : for (; i <= db; i++)
608 78149 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
609 2491400 : for (; i <= m; i++)
610 0 : uel(C, i) = 0;
611 : }
612 58331 : for (; j <= ea; j++) {
613 0 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
614 0 : for (i = 1; i <= da; i++)
615 0 : uel(C, i) = ucoeff(A, ma + i, na + j);
616 0 : for (; i <= m; i++)
617 0 : uel(C, i) = 0;
618 : }
619 59795 : for (; j <= eb; j++) {
620 1464 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
621 85631 : for (i = 1; i <= db; i++)
622 84167 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
623 1464 : for (; i <= m; i++)
624 0 : uel(C, i) = 0;
625 : }
626 59795 : for (; j <= n; j++)
627 1464 : gel(M, j) = zero_Flv(m);
628 58331 : return M;
629 : }
630 :
631 : static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
632 :
633 : /* Strassen-Winograd matrix product A (m x n) * B (n x p) */
634 : static GEN
635 8333 : Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
636 : {
637 : pari_sp av;
638 8333 : long m1 = (m + 1)/2, m2 = m/2,
639 8333 : n1 = (n + 1)/2, n2 = n/2,
640 8333 : p1 = (p + 1)/2, p2 = p/2;
641 : GEN A11, A12, A22, B11, B21, B22,
642 : S1, S2, S3, S4, T1, T2, T3, T4,
643 : M1, M2, M3, M4, M5, M6, M7,
644 : V1, V2, V3, C;
645 : long j;
646 8333 : C = cgetg(p + 1, t_MAT);
647 683199 : for (j = 1; j <= p; j++)
648 674867 : gel(C, j) = cgetg(m + 1, t_VECSMALL);
649 8332 : av = avma;
650 8332 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
651 8333 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
652 8333 : M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
653 8333 : if (gc_needed(av, 1))
654 0 : gerepileall(av, 2, &T2, &M2); /* destroy S1 */
655 8333 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
656 8333 : if (gc_needed(av, 1))
657 0 : gerepileall(av, 2, &M2, &T3); /* destroy T2 */
658 8333 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
659 8333 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
660 8333 : M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
661 8332 : if (gc_needed(av, 1))
662 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
663 8332 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
664 8333 : if (gc_needed(av, 1))
665 0 : gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
666 8333 : A11 = matslice(A, 1, m1, 1, n1);
667 8333 : B11 = matslice(B, 1, n1, 1, p1);
668 8333 : M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
669 8333 : if (gc_needed(av, 1))
670 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
671 8333 : A12 = matslice(A, 1, m1, n1 + 1, n);
672 8333 : B21 = matslice(B, n1 + 1, n, 1, p1);
673 8333 : M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
674 8333 : if (gc_needed(av, 1))
675 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
676 8333 : add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
677 8333 : if (gc_needed(av, 1))
678 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy M4 */
679 8333 : M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
680 8333 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
681 8333 : if (gc_needed(av, 1))
682 0 : gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4); /* destroy S3 */
683 8333 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
684 8333 : if (gc_needed(av, 1))
685 0 : gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4); /* destroy T3 */
686 8333 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
687 8333 : if (gc_needed(av, 1))
688 0 : gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1); /* destroy M1, M5 */
689 8333 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
690 8333 : M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
691 8333 : if (gc_needed(av, 1))
692 0 : gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6); /* destroy S4, B22 */
693 8333 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
694 8333 : M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
695 8333 : if (gc_needed(av, 1))
696 0 : gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7); /* destroy A22, T4 */
697 8333 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
698 8333 : add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
699 8333 : if (gc_needed(av, 1))
700 0 : gerepileall(av, 4, &M2, &M3, &V1, &M7); /* destroy V3, M6 */
701 8333 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
702 8333 : if (gc_needed(av, 1))
703 0 : gerepileall(av, 3, &M3, &M7, &V2); /* destroy V1, M2 */
704 8333 : add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
705 8333 : add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
706 8333 : set_avma(av); return C;
707 : }
708 :
709 : /* Strassen-Winograd used for dim >= ZM_sw_bound */
710 : static GEN
711 28000990 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
712 : {
713 28000990 : ulong e = expu(p);
714 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
715 23528407 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
716 : #else
717 4473242 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
718 : #endif
719 28001649 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
720 27993316 : return Flm_mul_classical(x, y, l, lx, ly, p, pi);
721 : else
722 8333 : return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
723 : }
724 :
725 : GEN
726 13637252 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
727 : {
728 13637252 : long lx=lg(x), ly=lg(y);
729 13637252 : if (ly==1) return cgetg(1,t_MAT);
730 13637252 : if (lx==1) return zero_Flm(0, ly-1);
731 13637252 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
732 : }
733 :
734 : GEN
735 14229329 : Flm_mul(GEN x, GEN y, ulong p)
736 : {
737 14229329 : long lx=lg(x), ly=lg(y);
738 14229329 : if (ly==1) return cgetg(1,t_MAT);
739 14229042 : if (lx==1) return zero_Flm(0, ly-1);
740 14229042 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
741 : }
742 :
743 : GEN
744 74945 : Flm_sqr(GEN x, ulong p)
745 : {
746 74945 : long lx = lg(x);
747 74945 : if (lx==1) return cgetg(1,t_MAT);
748 74945 : return Flm_mul_i(x, x, lx, lx, lx, p, get_Fl_red(p));
749 : }
750 :
751 : struct _Flm
752 : {
753 : ulong p;
754 : long n;
755 : };
756 :
757 : static GEN
758 18069 : _Flm_mul(void *E , GEN x, GEN y)
759 18069 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
760 : static GEN
761 74945 : _Flm_sqr(void *E, GEN x)
762 74945 : { return Flm_sqr(x,((struct _Flm*)E)->p); }
763 : static GEN
764 868 : _Flm_one(void *E)
765 868 : { return matid_Flm(((struct _Flm*)E)->n); }
766 : GEN
767 41430 : Flm_powu(GEN x, ulong n, ulong p)
768 : {
769 : struct _Flm d;
770 41430 : if (!n) return matid(lg(x)-1);
771 41430 : d.p = p;
772 41430 : return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
773 : }
774 : GEN
775 868 : Flm_powers(GEN x, ulong n, ulong p)
776 : {
777 : struct _Flm d;
778 868 : d.p = p;
779 868 : d.n = lg(x)-1;
780 868 : return gen_powers(x, n, 1, (void*)&d,
781 : &_Flm_sqr, &_Flm_mul, &_Flm_one);
782 : }
783 :
784 : static GEN
785 0 : _FpM_mul(void *p , GEN x, GEN y)
786 0 : { return FpM_mul(x,y,(GEN)p); }
787 : static GEN
788 0 : _FpM_sqr(void *p, GEN x)
789 0 : { return FpM_mul(x,x,(GEN)p); }
790 : GEN
791 0 : FpM_powu(GEN x, ulong n, GEN p)
792 : {
793 0 : if (!n) return matid(lg(x)-1);
794 0 : if (lgefint(p) == 3)
795 : {
796 0 : pari_sp av = avma;
797 0 : ulong pp = uel(p,2);
798 : GEN z;
799 0 : if (pp == 2)
800 0 : z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
801 : else
802 0 : z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
803 0 : return gerepileupto(av, z);
804 : }
805 0 : return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
806 : }
807 :
808 : /*Multiple a column vector by a line vector to make a matrix*/
809 : GEN
810 14 : Flc_Flv_mul(GEN x, GEN y, ulong p)
811 : {
812 14 : long i,j, lx=lg(x), ly=lg(y);
813 : GEN z;
814 14 : if (ly==1) return cgetg(1,t_MAT);
815 14 : z = cgetg(ly, t_MAT);
816 49 : for (j=1; j < ly; j++)
817 : {
818 35 : GEN zj = cgetg(lx,t_VECSMALL);
819 112 : for (i=1; i<lx; i++)
820 77 : uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
821 35 : gel(z,j) = zj;
822 : }
823 14 : return z;
824 : }
825 :
826 : /*Multiple a column vector by a line vector to make a matrix*/
827 : GEN
828 5669 : FpC_FpV_mul(GEN x, GEN y, GEN p)
829 : {
830 5669 : long i,j, lx=lg(x), ly=lg(y);
831 : GEN z;
832 5669 : if (ly==1) return cgetg(1,t_MAT);
833 5669 : z = cgetg(ly,t_MAT);
834 63924 : for (j=1; j < ly; j++)
835 : {
836 58257 : GEN zj = cgetg(lx,t_COL);
837 1379985 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
838 58255 : gel(z, j) = zj;
839 : }
840 5667 : return z;
841 : }
842 :
843 : /* Multiply a line vector by a column and return a scalar (t_INT) */
844 : GEN
845 8982007 : FpV_dotproduct(GEN x, GEN y, GEN p)
846 : {
847 8982007 : long i, lx = lg(x);
848 : pari_sp av;
849 : GEN c;
850 8982007 : if (lx == 1) return gen_0;
851 8979501 : av = avma; c = mulii(gel(x,1),gel(y,1));
852 37006213 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
853 8979476 : return gerepileuptoint(av, modii(c,p));
854 : }
855 : GEN
856 0 : FpV_dotsquare(GEN x, GEN p)
857 : {
858 0 : long i, lx = lg(x);
859 : pari_sp av;
860 : GEN c;
861 0 : if (lx == 1) return gen_0;
862 0 : av = avma; c = sqri(gel(x,1));
863 0 : for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
864 0 : return gerepileuptoint(av, modii(c,p));
865 : }
866 :
867 : INLINE ulong
868 9291126 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
869 : {
870 9291126 : ulong c = uel(x,0) * uel(y,0);
871 : long k;
872 70746994 : for (k = 1; k < lx; k++) {
873 61455868 : c += uel(x,k) * uel(y,k);
874 61455868 : if (c & HIGHBIT) c %= p;
875 : }
876 9291126 : return c % p;
877 : }
878 :
879 : INLINE ulong
880 740062 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
881 : {
882 : ulong l0, l1, v1, h0, h1;
883 740062 : long i = 0;
884 : LOCAL_OVERFLOW;
885 : LOCAL_HIREMAINDER;
886 740062 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
887 15410470 : while (++i < lx) {
888 14670408 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
889 14670408 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
890 : }
891 740062 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
892 464435 : else return remlll_pre(v1, h1, l1, p, pi);
893 : }
894 :
895 : ulong
896 631628 : Flv_dotproduct(GEN x, GEN y, ulong p)
897 : {
898 631628 : long lx = lg(x)-1;
899 631628 : if (lx == 0) return 0;
900 631628 : if (SMALL_ULONG(p))
901 631628 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
902 : else
903 0 : return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
904 : }
905 :
906 : ulong
907 58818 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
908 : {
909 58818 : long lx = lg(x)-1;
910 58818 : if (lx == 0) return 0;
911 58818 : if (SMALL_ULONG(p))
912 53260 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
913 : else
914 5558 : return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
915 : }
916 :
917 : ulong
918 9408465 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
919 : {
920 9408465 : long lx = minss(lgpol(x), lgpol(y));
921 9408478 : if (lx == 0) return 0;
922 9340738 : if (pi)
923 734504 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
924 : else
925 8606234 : return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
926 : }
927 : ulong
928 0 : Flx_dotproduct(GEN x, GEN y, ulong p)
929 0 : { return Flx_dotproduct_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
930 :
931 : GEN
932 1737468 : FpM_FpC_mul(GEN x, GEN y, GEN p)
933 : {
934 1737468 : long lx = lg(x);
935 1737468 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
936 : }
937 : GEN
938 977524 : Flm_Flc_mul(GEN x, GEN y, ulong p)
939 : {
940 977524 : long l, lx = lg(x);
941 977524 : if (lx==1) return cgetg(1,t_VECSMALL);
942 977524 : l = lgcols(x);
943 977523 : if (p==2)
944 276443 : return Flm_Flc_mul_i_2(x, y, lx, l);
945 701080 : else if (SMALL_ULONG(p))
946 463742 : return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
947 : else
948 237338 : return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
949 : }
950 :
951 : /* allow pi = 0 */
952 : GEN
953 6693 : Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
954 : {
955 6693 : long l, lx = lg(x);
956 : GEN z;
957 6693 : if (lx==1) return cgetg(1,t_VECSMALL);
958 6693 : l = lgcols(x);
959 6693 : z = cgetg(l, t_VECSMALL);
960 6693 : if (SMALL_ULONG(p))
961 4407 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
962 : else
963 2286 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
964 6693 : return z;
965 : }
966 :
967 : /* allow pi = 0 */
968 : GEN
969 7535246 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
970 : {
971 7535246 : long l, lx = lg(x);
972 : GEN z;
973 7535246 : if (lx==1) return pol0_Flx(sv);
974 7535246 : l = lgcols(x);
975 7533903 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
976 7533429 : if (SMALL_ULONG(p))
977 3304142 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
978 : else
979 4229287 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
980 7550345 : return Flx_renormalize(z, l + 1);
981 : }
982 :
983 : /* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
984 : GEN
985 1420753 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
986 : {
987 1420753 : long i, l, lx = lg(x);
988 : GEN z;
989 1420753 : if (lx==1) return pol_0(v);
990 1420753 : l = lgcols(x);
991 1420753 : z = new_chunk(l+1);
992 2162526 : for (i=l-1; i; i--)
993 : {
994 1979340 : pari_sp av = avma;
995 1979340 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
996 1979339 : p1 = modii(p1, p);
997 1979335 : if (signe(p1))
998 : {
999 1237562 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
1000 1237562 : gel(z,i+1) = gerepileuptoint(av, p1);
1001 1237567 : break;
1002 : }
1003 741773 : set_avma(av);
1004 : }
1005 1420753 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1006 1237567 : z[0] = evaltyp(t_POL) | _evallg(i+2);
1007 1237567 : z[1] = evalsigne(1) | evalvarn(v);
1008 3603698 : for (; i; i--)
1009 : {
1010 2366135 : pari_sp av = avma;
1011 2366135 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1012 2366132 : gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
1013 : }
1014 1237563 : return z;
1015 : }
1016 :
1017 : /********************************************************************/
1018 : /** **/
1019 : /** TRANSPOSITION **/
1020 : /** **/
1021 : /********************************************************************/
1022 :
1023 : /* == zm_transpose */
1024 : GEN
1025 726598 : Flm_transpose(GEN x)
1026 : {
1027 726598 : long i, dx, lx = lg(x);
1028 : GEN y;
1029 726598 : if (lx == 1) return cgetg(1,t_MAT);
1030 726465 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1031 9163109 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1032 726467 : return y;
1033 : }
1034 :
1035 : /********************************************************************/
1036 : /** **/
1037 : /** SCALAR MATRICES **/
1038 : /** **/
1039 : /********************************************************************/
1040 :
1041 : GEN
1042 4029 : gen_matid(long n, void *E, const struct bb_field *S)
1043 : {
1044 4029 : GEN y = cgetg(n+1,t_MAT), _0, _1;
1045 : long i;
1046 4029 : if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
1047 4029 : _0 = S->s(E,0);
1048 4029 : _1 = S->s(E,1);
1049 17952 : for (i=1; i<=n; i++)
1050 : {
1051 13923 : GEN z = const_col(n, _0); gel(z,i) = _1;
1052 13923 : gel(y, i) = z;
1053 : }
1054 4029 : return y;
1055 : }
1056 :
1057 : GEN
1058 35 : matid_F2xqM(long n, GEN T)
1059 : {
1060 : void *E;
1061 35 : const struct bb_field *S = get_F2xq_field(&E, T);
1062 35 : return gen_matid(n, E, S);
1063 : }
1064 : GEN
1065 56 : matid_FlxqM(long n, GEN T, ulong p)
1066 : {
1067 : void *E;
1068 56 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1069 56 : return gen_matid(n, E, S);
1070 : }
1071 :
1072 : GEN
1073 1421425 : matid_Flm(long n)
1074 : {
1075 1421425 : GEN y = cgetg(n+1,t_MAT);
1076 : long i;
1077 1421419 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1078 8995730 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1079 1421406 : return y;
1080 : }
1081 :
1082 : GEN
1083 42 : scalar_Flm(long s, long n)
1084 : {
1085 : long i;
1086 42 : GEN y = cgetg(n+1,t_MAT);
1087 560 : for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
1088 42 : return y;
1089 : }
1090 :
1091 : /********************************************************************/
1092 : /** **/
1093 : /** CONVERSIONS **/
1094 : /** **/
1095 : /********************************************************************/
1096 : GEN
1097 54538366 : ZV_to_Flv(GEN x, ulong p)
1098 677400684 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1099 :
1100 : GEN
1101 12777848 : ZM_to_Flm(GEN x, ulong p)
1102 66216860 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1103 :
1104 : GEN
1105 14722 : ZVV_to_FlvV(GEN x, ulong p)
1106 112431 : { pari_APPLY_type(t_VEC, ZV_to_Flv(gel(x,i), p)) }
1107 :
1108 : GEN
1109 3944 : ZMV_to_FlmV(GEN x, ulong m)
1110 34243 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1111 :
1112 : /* TO INTMOD */
1113 : static GEN
1114 20251575 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1115 : static GEN
1116 578850 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1117 :
1118 : GEN
1119 3603 : Fp_to_mod(GEN z, GEN p)
1120 : {
1121 3603 : retmkintmod(modii(z, p), icopy(p));
1122 : }
1123 :
1124 : GEN
1125 997248 : FpX_to_mod_raw(GEN z, GEN p)
1126 : {
1127 997248 : long i, l = lg(z);
1128 : GEN x;
1129 997248 : if (l == 2)
1130 : {
1131 522424 : x = cgetg(3,t_POL); x[1] = z[1];
1132 522424 : gel(x,2) = mkintmod(gen_0,p); return x;
1133 : }
1134 474824 : x = cgetg(l,t_POL); x[1] = z[1];
1135 3918509 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1136 474824 : return normalizepol_lg(x,l);
1137 : }
1138 :
1139 : /* z in Z[X], return z * Mod(1,p), normalized*/
1140 : GEN
1141 1669085 : FpX_to_mod(GEN z, GEN p)
1142 : {
1143 1669085 : long i, l = lg(z);
1144 : GEN x;
1145 1669085 : if (l == 2)
1146 : {
1147 3171 : x = cgetg(3,t_POL); x[1] = z[1];
1148 3171 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1149 : }
1150 1665914 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1151 18370815 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1152 1665914 : return normalizepol_lg(x,l);
1153 : }
1154 :
1155 : GEN
1156 84021 : FpXC_to_mod(GEN z, GEN p)
1157 : {
1158 84021 : long i,l = lg(z);
1159 84021 : GEN x = cgetg(l, t_COL);
1160 84021 : p = icopy(p);
1161 474166 : for (i=1; i<l; i++)
1162 390145 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1163 84021 : return x;
1164 : }
1165 :
1166 : static GEN
1167 0 : FpXC_to_mod_raw(GEN x, GEN p)
1168 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1169 :
1170 : GEN
1171 0 : FpXM_to_mod(GEN z, GEN p)
1172 : {
1173 0 : long i,l = lg(z);
1174 0 : GEN x = cgetg(l, t_MAT);
1175 0 : p = icopy(p);
1176 0 : for (i=1; i<l; i++)
1177 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1178 0 : return x;
1179 : }
1180 :
1181 : /* z in Z^n, return z * Mod(1,p), normalized*/
1182 : GEN
1183 35984 : FpV_to_mod(GEN z, GEN p)
1184 : {
1185 35984 : long i,l = lg(z);
1186 35984 : GEN x = cgetg(l, t_VEC);
1187 35984 : if (l == 1) return x;
1188 35984 : p = icopy(p);
1189 72367 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1190 35984 : return x;
1191 : }
1192 : /* z in Z^n, return z * Mod(1,p), normalized*/
1193 : GEN
1194 105 : FpC_to_mod(GEN z, GEN p)
1195 : {
1196 105 : long i, l = lg(z);
1197 105 : GEN x = cgetg(l, t_COL);
1198 105 : if (l == 1) return x;
1199 91 : p = icopy(p);
1200 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1201 91 : return x;
1202 : }
1203 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1204 : GEN
1205 226 : FpM_to_mod(GEN z, GEN p)
1206 : {
1207 226 : long i, j, m, l = lg(z);
1208 226 : GEN x = cgetg(l,t_MAT), y, zi;
1209 226 : if (l == 1) return x;
1210 205 : m = lgcols(z);
1211 205 : p = icopy(p);
1212 2456 : for (i=1; i<l; i++)
1213 : {
1214 2251 : gel(x,i) = cgetg(m,t_COL);
1215 2251 : y = gel(x,i); zi= gel(z,i);
1216 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1217 : }
1218 205 : return x;
1219 : }
1220 : GEN
1221 28 : Flc_to_mod(GEN z, ulong pp)
1222 : {
1223 28 : long i, l = lg(z);
1224 28 : GEN p, x = cgetg(l, t_COL);
1225 28 : if (l == 1) return x;
1226 28 : p = utoipos(pp);
1227 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1228 28 : return x;
1229 : }
1230 : GEN
1231 59694 : Flm_to_mod(GEN z, ulong pp)
1232 : {
1233 59694 : long i, j, m, l = lg(z);
1234 59694 : GEN p, x = cgetg(l,t_MAT), y, zi;
1235 59694 : if (l == 1) return x;
1236 59666 : m = lgcols(z);
1237 59666 : p = utoipos(pp);
1238 201328 : for (i=1; i<l; i++)
1239 : {
1240 141662 : gel(x,i) = cgetg(m,t_COL);
1241 141662 : y = gel(x,i); zi= gel(z,i);
1242 720428 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1243 : }
1244 59666 : return x;
1245 : }
1246 :
1247 : GEN
1248 574 : FpVV_to_mod(GEN z, GEN p)
1249 : {
1250 574 : long i, j, m, l = lg(z);
1251 574 : GEN x = cgetg(l,t_VEC), y, zi;
1252 574 : if (l == 1) return x;
1253 574 : m = lgcols(z);
1254 574 : p = icopy(p);
1255 1246 : for (i=1; i<l; i++)
1256 : {
1257 672 : gel(x,i) = cgetg(m,t_VEC);
1258 672 : y = gel(x,i); zi= gel(z,i);
1259 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1260 : }
1261 574 : return x;
1262 : }
1263 :
1264 : /* z in Z^n, return z * Mod(1,p), normalized*/
1265 : GEN
1266 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1267 : {
1268 7 : long i,l = lg(z);
1269 7 : GEN x = cgetg(l, t_COL);
1270 7 : if (l == 1) return x;
1271 7 : p = icopy(p);
1272 7 : T = FpX_to_mod_raw(T, p);
1273 21 : for (i=1; i<l; i++)
1274 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1275 7 : return x;
1276 : }
1277 :
1278 : static GEN
1279 582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1280 : {
1281 582533 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1282 582533 : return mkpolmod(z, T);
1283 : }
1284 :
1285 : /* z in Z^n, return z * Mod(1,p), normalized*/
1286 : GEN
1287 28 : FqC_to_mod(GEN z, GEN T, GEN p)
1288 : {
1289 : GEN x;
1290 28 : long i,l = lg(z);
1291 28 : if (!T) return FpC_to_mod(z, p);
1292 28 : x = cgetg(l, t_COL);
1293 28 : if (l == 1) return x;
1294 28 : p = icopy(p);
1295 28 : T = FpX_to_mod_raw(T, p);
1296 504 : for (i=1; i<l; i++)
1297 476 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1298 28 : return x;
1299 : }
1300 :
1301 : /* z in Z^n, return z * Mod(1,p), normalized*/
1302 : static GEN
1303 107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1304 690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1305 :
1306 : /* z in Z^n, return z * Mod(1,p), normalized*/
1307 : GEN
1308 26145 : FqM_to_mod(GEN z, GEN T, GEN p)
1309 : {
1310 : GEN x;
1311 26145 : long i,l = lg(z);
1312 26145 : if (!T) return FpM_to_mod(z, p);
1313 26145 : x = cgetg(l, t_MAT);
1314 26145 : if (l == 1) return x;
1315 26117 : p = icopy(p);
1316 26117 : T = FpX_to_mod_raw(T, p);
1317 134092 : for (i=1; i<l; i++)
1318 107975 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1319 26117 : return x;
1320 : }
1321 :
1322 : /********************************************************************/
1323 : /* */
1324 : /* Blackbox linear algebra */
1325 : /* */
1326 : /********************************************************************/
1327 :
1328 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1329 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1330 : * (e_j) is the canonical basis.
1331 : * A sparse matrix (zMs) is a t_VEC of zCs */
1332 :
1333 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1334 : * integer representing an element of Fp. This is important since p can be
1335 : * large and p+E[i] would not fit in a C long. */
1336 :
1337 : /* RgCs and RgMs are similar, except that the type of the component is
1338 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1339 : * of E. */
1340 :
1341 : /* Most functions take an argument nbrow which is the number of lines of the
1342 : * column/matrix, which cannot be derived from the data. */
1343 :
1344 : GEN
1345 0 : zCs_to_ZC(GEN R, long nbrow)
1346 : {
1347 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1348 0 : long j, l = lg(C);
1349 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1350 0 : return c;
1351 : }
1352 :
1353 : GEN
1354 0 : zMs_to_ZM(GEN x, long nbrow)
1355 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1356 :
1357 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1358 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1359 : GEN
1360 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1361 : {
1362 147 : pari_sp ltop = avma;
1363 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1364 147 : if (ZV_equal0(B)) return zerocol(n);
1365 147 : while (++col <= n)
1366 : {
1367 147 : pari_sp btop = avma, av;
1368 : long i, lQ;
1369 147 : GEN V, Q, M, W = B;
1370 147 : GEN b = cgetg(m+2, t_POL);
1371 147 : b[1] = evalsigne(1)|evalvarn(0);
1372 147 : gel(b, 2) = gel(W, col);
1373 46215 : for (i = 3; i<m+2; i++)
1374 46068 : gel(b, i) = cgeti(lgefint(p));
1375 147 : av = avma;
1376 46215 : for (i = 3; i<m+2; i++)
1377 : {
1378 46068 : W = f(E, W);
1379 46068 : affii(gel(W, col),gel(b, i));
1380 46068 : if (gc_needed(av,1))
1381 : {
1382 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1383 72 : W = gerepileupto(av, W);
1384 : }
1385 : }
1386 147 : b = FpX_renormalize(b, m+2);
1387 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1388 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1389 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1390 147 : W = B; lQ =lg(Q);
1391 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1392 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1393 147 : av = avma;
1394 22740 : for (i = lQ-3; i > 1; i--)
1395 : {
1396 22593 : W = f(E, W);
1397 22593 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1398 22593 : if (gc_needed(av,1))
1399 : {
1400 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1401 110 : gerepileall(av, 2, &V, &W);
1402 : }
1403 : }
1404 147 : V = FpC_red(V, p);
1405 147 : W = FpC_sub(f(E,V), B, p);
1406 147 : if (ZV_equal0(W)) return gerepilecopy(ltop, V);
1407 0 : av = avma;
1408 0 : for (i = 1; i <= n; ++i)
1409 : {
1410 0 : V = W;
1411 0 : W = f(E, V);
1412 0 : if (ZV_equal0(W))
1413 0 : return gerepilecopy(ltop, shallowtrans(V));
1414 0 : gerepileall(av, 2, &V, &W);
1415 : }
1416 0 : set_avma(btop);
1417 : }
1418 0 : return NULL;
1419 : }
1420 :
1421 : GEN
1422 0 : zMs_ZC_mul(GEN M, GEN B)
1423 : {
1424 : long i, j;
1425 0 : long n = lg(B)-1;
1426 0 : GEN V = zerocol(n);
1427 0 : for (i = 1; i <= n; ++i)
1428 0 : if (signe(gel(B, i)))
1429 : {
1430 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1431 0 : long l = lg(C);
1432 0 : for (j = 1; j < l; ++j)
1433 : {
1434 0 : long k = C[j];
1435 0 : switch(E[j])
1436 : {
1437 0 : case 1:
1438 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1439 0 : break;
1440 0 : case -1:
1441 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1442 0 : break;
1443 0 : default:
1444 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]));
1445 0 : break;
1446 : }
1447 : }
1448 : }
1449 0 : return V;
1450 : }
1451 :
1452 : GEN
1453 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1454 :
1455 : GEN
1456 68955 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1457 : {
1458 68955 : long i, j, lM = lg(M);
1459 68955 : GEN V = cgetg(lM,t_VEC);
1460 28331169 : for (i = 1; i < lM; ++i)
1461 : {
1462 28262214 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1463 28262214 : pari_sp av = avma;
1464 28262214 : long lC = lg(C);
1465 28262214 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1466 28262214 : z = mulis(gel(B, C[1]), E[1]);
1467 228722330 : for (j = 2; j < lC; ++j)
1468 : {
1469 200460116 : GEN b = gel(B, C[j]);
1470 200460116 : switch(E[j])
1471 : {
1472 131566519 : case 1: z = addii(z, b); break;
1473 57521817 : case -1: z = subii(z, b); break;
1474 11371780 : default: z = addii(z, mulis(b, E[j])); break;
1475 : }
1476 : }
1477 28262214 : gel(V,i) = gerepileuptoint(av, p? Fp_red(z, p): z);
1478 : }
1479 68955 : return V;
1480 : }
1481 : GEN
1482 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1483 :
1484 : GEN
1485 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1486 : {
1487 0 : pari_sp av = avma, av2;
1488 0 : GEN xi, xb, pi = gen_1, P;
1489 : long i;
1490 0 : if (!C) {
1491 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1492 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1493 : }
1494 0 : P = utoipos(p);
1495 0 : av2 = avma;
1496 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1497 0 : xb = Flm_to_ZM(xi);
1498 0 : for (i = 2; i <= e; i++)
1499 : {
1500 0 : pi = muliu(pi, p); /* = p^(i-1) */
1501 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1502 0 : if (gc_needed(av,2))
1503 : {
1504 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1505 0 : gerepileall(av2,3, &pi,&b,&xb);
1506 : }
1507 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1508 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1509 : }
1510 0 : return gerepileupto(av, xb);
1511 : }
1512 :
1513 : struct wrapper_modp_s {
1514 : GEN (*f)(void*, GEN);
1515 : void *E;
1516 : GEN p;
1517 : };
1518 :
1519 : /* compute f(x) mod p */
1520 : static GEN
1521 0 : wrap_relcomb_modp(void *E, GEN x)
1522 : {
1523 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1524 0 : return FpC_red(W->f(W->E, x), W->p);
1525 : }
1526 : static GEN
1527 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1528 :
1529 : static GEN
1530 68808 : wrap_relker(void*E, GEN x)
1531 : {
1532 68808 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1533 68808 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1534 : }
1535 :
1536 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1537 : GEN
1538 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1539 : {
1540 : struct wrapper_modp_s W;
1541 0 : pari_sp av = avma;
1542 0 : GEN xb, xi, pi = gen_1;
1543 : long i;
1544 0 : W.E = E;
1545 0 : W.f = f;
1546 0 : W.p = p;
1547 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1548 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1549 0 : xb = xi;
1550 0 : for (i = 2; i <= e; i++)
1551 : {
1552 0 : pi = mulii(pi, p); /* = p^(i-1) */
1553 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1554 0 : if (gc_needed(av,2))
1555 : {
1556 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1557 0 : gerepileall(av,3, &pi,&B,&xb);
1558 : }
1559 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1560 0 : if (!xi) return NULL;
1561 0 : if (typ(xi) == t_VEC) return gerepileupto(av, xi);
1562 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1563 : }
1564 0 : return gerepileupto(av, xb);
1565 : }
1566 :
1567 : static GEN
1568 23034 : vecprow(GEN A, GEN prow)
1569 : {
1570 23034 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1571 : }
1572 :
1573 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1574 : * or the index of a column which is linearly dependent from the others as a
1575 : * t_VECSMALL with a single component. */
1576 : GEN
1577 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1578 : {
1579 0 : pari_sp av = avma;
1580 : GEN pcol, prow;
1581 0 : long nbi=lg(M)-1, lR;
1582 : long i, n;
1583 : GEN Mp, Ap, Rp;
1584 : pari_timer ti;
1585 0 : if (DEBUGLEVEL) timer_start(&ti);
1586 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1587 0 : if (!pcol) return gc_NULL(av);
1588 0 : if (DEBUGLEVEL)
1589 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1590 0 : n = lg(pcol)-1;
1591 0 : Mp = cgetg(n+1, t_MAT);
1592 0 : for(i=1; i<=n; i++)
1593 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1594 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1595 0 : if (DEBUGLEVEL) timer_start(&ti);
1596 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1597 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1598 0 : if (!Rp) return gc_NULL(av);
1599 0 : lR = lg(Rp)-1;
1600 0 : if (typ(Rp) == t_COL)
1601 : {
1602 0 : GEN R = zerocol(nbi+1);
1603 0 : for (i=1; i<=lR; i++)
1604 0 : gel(R,pcol[i]) = gel(Rp,i);
1605 0 : return gerepilecopy(av, R);
1606 : }
1607 0 : for (i = 1; i <= lR; ++i)
1608 0 : if (signe(gel(Rp, i)))
1609 0 : return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
1610 : return NULL; /* LCOV_EXCL_LINE */
1611 : }
1612 :
1613 : GEN
1614 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1615 : {
1616 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1617 : }
1618 :
1619 : GEN
1620 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1621 : {
1622 : GEN res;
1623 0 : pari_CATCH(e_INV)
1624 : {
1625 0 : GEN E = pari_err_last();
1626 0 : GEN x = err_get_compo(E,2);
1627 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1628 0 : if (DEBUGLEVEL)
1629 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1630 0 : res = NULL;
1631 : } pari_TRY {
1632 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1633 0 : } pari_ENDCATCH
1634 0 : return res;
1635 : }
1636 :
1637 : static GEN
1638 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1639 : {
1640 147 : long i, j, oldf = 0, f = 0;
1641 147 : long lrow = lg(prow), lM = lg(M);
1642 147 : GEN W = const_vecsmall(lM-1,1);
1643 147 : GEN R = cgetg(lrow, t_VEC);
1644 228354 : for (i=1; i<lrow; i++)
1645 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1646 : do
1647 : {
1648 441 : oldf = f;
1649 374479 : for(i=1; i<lM; i++)
1650 374038 : if (W[i])
1651 : {
1652 131759 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1653 131759 : long c=0, cj=0, lF = lg(F);
1654 1093879 : for(j=1; j<lF; j++)
1655 962120 : if (!gel(R,F[j])) { c++; cj=j; }
1656 131759 : if (c>=2) continue;
1657 112310 : if (c==1)
1658 : {
1659 32901 : pari_sp av = avma;
1660 32901 : GEN s = gen_0;
1661 272970 : for(j=1; j<lF; j++)
1662 240069 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1663 : /* s /= -E[cj] mod p */
1664 32901 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1665 32901 : gel(R,F[cj]) = gerepileuptoint(av, s);
1666 32901 : f++;
1667 : }
1668 112310 : W[i]=0;
1669 : }
1670 441 : } while(oldf!=f);
1671 228354 : for (i=1; i<lrow; i++)
1672 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1673 147 : return R;
1674 : }
1675 :
1676 : /* Return a linear form Y such that YM is essentially 0 */
1677 : GEN
1678 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1679 : {
1680 147 : pari_sp av = avma, av2;
1681 : GEN Mp, pcol, prow;
1682 : long i, n;
1683 : pari_timer ti;
1684 : struct wrapper_modp_s W;
1685 147 : if (DEBUGLEVEL) timer_start(&ti);
1686 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1687 147 : if (!pcol) return gc_NULL(av);
1688 147 : if (DEBUGLEVEL)
1689 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1690 147 : n = lg(pcol)-1;
1691 147 : Mp = cgetg(n+1, t_MAT);
1692 23181 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1693 147 : W.E = (void*)Mp;
1694 147 : W.p = p;
1695 147 : av2 = avma;
1696 0 : for(;; set_avma(av2))
1697 0 : {
1698 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1699 147 : if (DEBUGLEVEL) timer_start(&ti);
1700 147 : pari_CATCH(e_INV)
1701 : {
1702 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1703 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1704 0 : if (DEBUGLEVEL)
1705 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1706 0 : Rp = NULL;
1707 : } pari_TRY {
1708 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1709 147 : } pari_ENDCATCH
1710 147 : if (!Rp) continue;
1711 147 : if (typ(Rp)==t_COL)
1712 : {
1713 147 : Rp = FpC_sub(Rp,B,p);
1714 147 : if (ZV_equal0(Rp)) continue;
1715 : }
1716 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1717 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1718 147 : return gerepilecopy(av, R);
1719 : }
1720 : }
1721 :
1722 : GEN
1723 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1724 : {
1725 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1726 : }
|
