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 29384775 : FpC_red(GEN x, GEN p)
28 252199190 : { 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 6864881 : FpC_center(GEN x, GEN p, GEN pov2)
37 43233554 : { 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 5155895 : FpM_red(GEN x, GEN p)
46 23832678 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2756114 : FpM_center(GEN x, GEN p, GEN pov2)
50 8300137 : { 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 23183 : 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 41335 : random_Flv(long n, ulong p)
112 : {
113 41335 : GEN y = cgetg(n+1, t_VECSMALL);
114 : long i;
115 261421 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
116 41335 : return y;
117 : }
118 :
119 : /********************************************************************/
120 : /** **/
121 : /** ADD, SUB **/
122 : /** **/
123 : /********************************************************************/
124 : GEN
125 279684 : FpC_add(GEN x, GEN y, GEN p)
126 3620368 : { 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 1663577 : Flv_add(GEN x, GEN y, ulong p)
138 : {
139 1663577 : if (p==2)
140 6520572 : pari_APPLY_ulong(x[i]^y[i])
141 : else
142 38488279 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
143 : }
144 :
145 : void
146 441831 : Flv_add_inplace(GEN x, GEN y, ulong p)
147 : {
148 441831 : long i, l = lg(x);
149 441831 : if (p==2)
150 1373130 : 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 441831 : }
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 932161 : FpC_sub(GEN x, GEN y, GEN p)
169 18883406 : { 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 227845433 : Flv_sub(GEN x, GEN y, ulong p)
181 1023348862 : { 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 51126 : Flm_Fl_add(GEN x, ulong y, ulong p)
192 : {
193 51126 : long l = lg(x), i, j;
194 51126 : GEN z = cgetg(l,t_MAT);
195 :
196 51126 : if (l==1) return z;
197 51126 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
198 240837 : for (i=1; i<l; i++)
199 : {
200 189709 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
201 189711 : gel(z,i) = zi;
202 1318340 : for (j=1; j<l; j++) zi[j] = xi[j];
203 189711 : zi[i] = Fl_add(zi[i], y, p);
204 : }
205 51128 : 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 25315 : Flm_add(GEN x, GEN y, ulong p)
212 617992 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
213 :
214 : GEN
215 25481787 : Flm_sub(GEN x, GEN y, ulong p)
216 253306619 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
217 :
218 : /********************************************************************/
219 : /** **/
220 : /** MULTIPLICATION **/
221 : /** **/
222 : /********************************************************************/
223 : GEN
224 947403 : FpC_Fp_mul(GEN x, GEN y, GEN p)
225 11731062 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
226 :
227 : GEN
228 1578670 : Flv_Fl_mul(GEN x, ulong y, ulong p)
229 37247529 : { 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 761966 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
243 761966 : long i, j, h, l = lg(X);
244 761966 : GEN A = cgetg(l, t_MAT);
245 761968 : if (l == 1) return A;
246 761968 : h = lgcols(X);
247 4039158 : for (j=1; j<l; j++)
248 : {
249 3277345 : GEN a = cgetg(h, t_COL), x = gel(X, j);
250 23869350 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
251 3277191 : gel(A,j) = a;
252 : }
253 761813 : return A;
254 : }
255 :
256 : /* x *= y */
257 : void
258 7733562 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
259 : {
260 : long i;
261 7733562 : if (HIGHWORD(y | p))
262 8556849 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
263 : else
264 31957406 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
265 7733562 : }
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 16981058 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
288 : {
289 16981058 : long i, j, m, l = lg(y);
290 16981058 : GEN z = cgetg(l, t_MAT);
291 16979139 : if (l == 1) return z;
292 16979139 : m = lgcols(y);
293 151990572 : for(j=1; j<l; j++) {
294 134914879 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
295 362614529 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
296 : }
297 17075693 : return z;
298 : }
299 :
300 : /* return x * y */
301 : static GEN
302 8517845 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
303 : {
304 8517845 : long i, j, m, l = lg(y);
305 8517845 : GEN z = cgetg(l, t_MAT);
306 8517841 : if (l == 1) return z;
307 8517841 : m = lgcols(y);
308 52658103 : for(j=1; j<l; j++) {
309 44140275 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 157402627 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
311 : }
312 8517828 : return z;
313 : }
314 :
315 : /* return x * y */
316 : GEN
317 25430253 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
318 : {
319 25430253 : if (HIGHWORD(p))
320 16980956 : return Flm_Fl_mul_pre_i(y, x, p, pi);
321 : else
322 8449297 : return Flm_Fl_mul_OK(y, x, p);
323 : }
324 :
325 : /* return x * y */
326 : GEN
327 68593 : Flm_Fl_mul(GEN y, ulong x, ulong p)
328 : {
329 68593 : if (HIGHWORD(p))
330 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
331 : else
332 68519 : return Flm_Fl_mul_OK(y, x, p);
333 : }
334 :
335 : GEN
336 4064738 : Flv_neg(GEN x, ulong p)
337 31147143 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
338 :
339 : void
340 6193 : Flv_neg_inplace(GEN v, ulong p)
341 : {
342 : long i;
343 187636 : for (i = 1; i < lg(v); ++i)
344 181443 : v[i] = Fl_neg(v[i], p);
345 6193 : }
346 :
347 : GEN
348 320494 : Flm_neg(GEN x, ulong p)
349 4354436 : { 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 19389276 : ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
354 : {
355 19389276 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
356 : long k;
357 233106724 : for (k = 2; k < lx; k++)
358 : {
359 213719043 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
360 213715083 : if (signe(t)) c = addii(c, t);
361 : }
362 19387681 : 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 759520322 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
406 : {
407 759520322 : ulong c = ucoeff(x,i,1) * uel(y,1);
408 : long k;
409 11616671164 : for (k = 2; k < lx; k++) {
410 10857150842 : c += ucoeff(x,i,k) * uel(y,k);
411 10857150842 : if (c & HIGHBIT) c %= p;
412 : }
413 759520322 : return c % p;
414 : }
415 :
416 : static ulong
417 531118746 : 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 531118746 : long k = 1;
421 : LOCAL_OVERFLOW;
422 : LOCAL_HIREMAINDER;
423 531118746 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
424 8963584435 : while (++k < lx) {
425 8432465689 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
426 8432465689 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
427 : }
428 531118746 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
429 54150889 : else return remlll_pre(v1, h1, l1, p, pi);
430 : }
431 :
432 : static GEN
433 275941 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
434 : {
435 : long i,j;
436 275941 : GEN z = NULL;
437 :
438 3223328 : for (j=1; j<lx; j++)
439 : {
440 2947387 : if (!y[j]) continue;
441 1039496 : if (!z) z = Flv_copy(gel(x,j));
442 29407149 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
443 : }
444 275941 : if (!z) z = zero_zv(l-1);
445 275941 : return z;
446 : }
447 :
448 : static GEN
449 1736724 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
450 : {
451 1736724 : GEN z = cgetg(l,t_COL);
452 : long i;
453 16780449 : for (i = 1; i < l; i++)
454 : {
455 15043779 : pari_sp av = avma;
456 15043779 : GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
457 15043505 : gel(z,i) = gerepileuptoint(av, modii(c,p));
458 : }
459 1736670 : return z;
460 : }
461 :
462 : static void
463 85924115 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
464 : {
465 : long i;
466 845462407 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
467 85939874 : }
468 : static GEN
469 82614002 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
470 : {
471 82614002 : GEN z = cgetg(l,t_VECSMALL);
472 82613095 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
473 82613877 : return z;
474 : }
475 :
476 : static void
477 93999505 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
478 : {
479 : long i;
480 627856974 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
481 94355643 : }
482 : static GEN
483 89896118 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
484 : {
485 89896118 : GEN z = cgetg(l,t_VECSMALL);
486 89758446 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
487 89952014 : return z;
488 : }
489 :
490 : GEN
491 1543942 : FpM_mul(GEN x, GEN y, GEN p)
492 : {
493 1543942 : long lx=lg(x), ly=lg(y);
494 : GEN z;
495 1543942 : pari_sp av = avma;
496 1543942 : if (ly==1) return cgetg(1,t_MAT);
497 1543942 : if (lx==1) return zeromat(0, ly-1);
498 1543942 : if (lgefint(p) == 3)
499 : {
500 1498046 : ulong pp = uel(p,2);
501 1498046 : if (pp == 2)
502 : {
503 697904 : x = ZM_to_F2m(x);
504 697965 : y = ZM_to_F2m(y);
505 697969 : z = F2m_to_ZM(F2m_mul(x,y));
506 : }
507 : else
508 : {
509 800142 : x = ZM_to_Flm(x, pp);
510 800155 : y = ZM_to_Flm(y, pp);
511 800157 : z = Flm_to_ZM(Flm_mul(x,y, pp));
512 : }
513 : } else
514 45896 : z = FpM_red(ZM_mul(x, y), p);
515 1544009 : return gerepileupto(av, z);
516 : }
517 :
518 : static GEN
519 27994101 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
520 : {
521 : long j;
522 27994101 : GEN z = cgetg(ly, t_MAT);
523 27993101 : if (SMALL_ULONG(p))
524 102233448 : for (j=1; j<ly; j++)
525 82150400 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
526 : else
527 97608466 : for (j=1; j<ly; j++)
528 89696804 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
529 27994710 : 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 2807355 : for (j = 1; j <= min_e; j++) {
544 2740691 : C = gel(M, j + dx) + dy;
545 125621088 : for (i = 1; i <= min_d; i++)
546 122880397 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
547 122880453 : ucoeff(B, mb + i, nb + j), p);
548 2808255 : for (; i <= da; i++)
549 67620 : uel(C, i) = ucoeff(A, ma + i, na + j);
550 2740635 : for (; i <= db; i++)
551 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
552 2740635 : 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 1428286 : for (j = 1; j <= n; j++)
585 1394954 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
586 33332 : 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 2549646 : for (j = 1; j <= min_e; j++) {
601 2491315 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
602 117609735 : for (i = 1; i <= min_d; i++)
603 115118413 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
604 115118819 : ucoeff(B, mb + i, nb + j), p);
605 2594417 : for (; i <= da; i++)
606 103501 : uel(C, i) = ucoeff(A, ma + i, na + j);
607 2569060 : for (; i <= db; i++)
608 78144 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
609 2490916 : 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 683249 : for (j = 1; j <= p; j++)
648 674916 : gel(C, j) = cgetg(m + 1, t_VECSMALL);
649 8333 : av = avma;
650 8333 : 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 8333 : if (gc_needed(av, 1))
662 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
663 8333 : 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 28001205 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
712 : {
713 28001205 : ulong e = expu(p);
714 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
715 23528702 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
716 : #else
717 4473065 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
718 : #endif
719 28001767 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
720 27993434 : 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 13637360 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
727 : {
728 13637360 : long lx=lg(x), ly=lg(y);
729 13637360 : if (ly==1) return cgetg(1,t_MAT);
730 13637360 : if (lx==1) return zero_Flm(0, ly-1);
731 13637360 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
732 : }
733 :
734 : GEN
735 14229608 : Flm_mul(GEN x, GEN y, ulong p)
736 : {
737 14229608 : long lx=lg(x), ly=lg(y);
738 14229608 : if (ly==1) return cgetg(1,t_MAT);
739 14229321 : if (lx==1) return zero_Flm(0, ly-1);
740 14229321 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
741 : }
742 :
743 : GEN
744 74763 : Flm_sqr(GEN x, ulong p)
745 : {
746 74763 : long lx = lg(x);
747 74763 : if (lx==1) return cgetg(1,t_MAT);
748 74763 : 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 18055 : _Flm_mul(void *E , GEN x, GEN y)
759 18055 : { return Flm_mul(x,y,((struct _Flm*)E)->p); }
760 : static GEN
761 74763 : _Flm_sqr(void *E, GEN x)
762 74763 : { 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 41290 : Flm_powu(GEN x, ulong n, ulong p)
768 : {
769 : struct _Flm d;
770 41290 : if (!n) return matid(lg(x)-1);
771 41290 : d.p = p;
772 41290 : 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 63926 : for (j=1; j < ly; j++)
835 : {
836 58257 : GEN zj = cgetg(lx,t_COL);
837 1379976 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
838 58257 : gel(z, j) = zj;
839 : }
840 5669 : return z;
841 : }
842 :
843 : /* Multiply a line vector by a column and return a scalar (t_INT) */
844 : GEN
845 8980585 : FpV_dotproduct(GEN x, GEN y, GEN p)
846 : {
847 8980585 : long i, lx = lg(x);
848 : pari_sp av;
849 : GEN c;
850 8980585 : if (lx == 1) return gen_0;
851 8978079 : av = avma; c = mulii(gel(x,1),gel(y,1));
852 37001955 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
853 8978034 : 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 9290373 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
869 : {
870 9290373 : ulong c = uel(x,0) * uel(y,0);
871 : long k;
872 70742053 : for (k = 1; k < lx; k++) {
873 61451680 : c += uel(x,k) * uel(y,k);
874 61451680 : if (c & HIGHBIT) c %= p;
875 : }
876 9290373 : 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 631068 : Flv_dotproduct(GEN x, GEN y, ulong p)
897 : {
898 631068 : long lx = lg(x)-1;
899 631068 : if (lx == 0) return 0;
900 631068 : if (SMALL_ULONG(p))
901 631068 : 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 58817 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
908 : {
909 58817 : long lx = lg(x)-1;
910 58817 : if (lx == 0) return 0;
911 58817 : if (SMALL_ULONG(p))
912 53259 : 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 9408222 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
919 : {
920 9408222 : long lx = minss(lgpol(x), lgpol(y));
921 9408318 : if (lx == 0) return 0;
922 9340552 : if (pi)
923 734504 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
924 : else
925 8606048 : 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 1736731 : FpM_FpC_mul(GEN x, GEN y, GEN p)
933 : {
934 1736731 : long lx = lg(x);
935 1736731 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
936 : }
937 : GEN
938 976963 : Flm_Flc_mul(GEN x, GEN y, ulong p)
939 : {
940 976963 : long l, lx = lg(x);
941 976963 : if (lx==1) return cgetg(1,t_VECSMALL);
942 976963 : l = lgcols(x);
943 976963 : if (p==2)
944 275941 : return Flm_Flc_mul_i_2(x, y, lx, l);
945 701022 : else if (SMALL_ULONG(p))
946 463684 : 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 6692 : if (SMALL_ULONG(p))
961 4406 : __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 7539360 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
970 : {
971 7539360 : long l, lx = lg(x);
972 : GEN z;
973 7539360 : if (lx==1) return pol0_Flx(sv);
974 7539360 : l = lgcols(x);
975 7538322 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
976 7537171 : if (SMALL_ULONG(p))
977 3305389 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
978 : else
979 4231782 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
980 7547201 : 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 1420752 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
986 : {
987 1420752 : long i, l, lx = lg(x);
988 : GEN z;
989 1420752 : if (lx==1) return pol_0(v);
990 1420752 : l = lgcols(x);
991 1420753 : z = new_chunk(l+1);
992 2162530 : for (i=l-1; i; i--)
993 : {
994 1979344 : pari_sp av = avma;
995 1979344 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
996 1979344 : p1 = modii(p1, p);
997 1979345 : if (signe(p1))
998 : {
999 1237567 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
1000 1237567 : gel(z,i+1) = gerepileuptoint(av, p1);
1001 1237567 : break;
1002 : }
1003 741778 : 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 3603697 : for (; i; i--)
1009 : {
1010 2366130 : pari_sp av = avma;
1011 2366130 : GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
1012 2366125 : gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
1013 : }
1014 1237567 : return z;
1015 : }
1016 :
1017 : /********************************************************************/
1018 : /** **/
1019 : /** TRANSPOSITION **/
1020 : /** **/
1021 : /********************************************************************/
1022 :
1023 : /* == zm_transpose */
1024 : GEN
1025 726584 : Flm_transpose(GEN x)
1026 : {
1027 726584 : long i, dx, lx = lg(x);
1028 : GEN y;
1029 726584 : if (lx == 1) return cgetg(1,t_MAT);
1030 726451 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1031 9162954 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1032 726453 : 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 1421163 : matid_Flm(long n)
1074 : {
1075 1421163 : GEN y = cgetg(n+1,t_MAT);
1076 : long i;
1077 1421154 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1078 8994638 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1079 1421147 : 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 54531964 : ZV_to_Flv(GEN x, ulong p)
1098 677359596 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1099 :
1100 : GEN
1101 12789350 : ZM_to_Flm(GEN x, ulong p)
1102 66316595 : { pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1103 :
1104 : GEN
1105 3916 : ZMV_to_FlmV(GEN x, ulong m)
1106 34103 : { pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1107 :
1108 : /* TO INTMOD */
1109 : static GEN
1110 20251575 : to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1111 : static GEN
1112 578850 : Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1113 :
1114 : GEN
1115 3437 : Fp_to_mod(GEN z, GEN p)
1116 : {
1117 3437 : retmkintmod(modii(z, p), icopy(p));
1118 : }
1119 :
1120 : GEN
1121 997248 : FpX_to_mod_raw(GEN z, GEN p)
1122 : {
1123 997248 : long i, l = lg(z);
1124 : GEN x;
1125 997248 : if (l == 2)
1126 : {
1127 522424 : x = cgetg(3,t_POL); x[1] = z[1];
1128 522424 : gel(x,2) = mkintmod(gen_0,p); return x;
1129 : }
1130 474824 : x = cgetg(l,t_POL); x[1] = z[1];
1131 3918509 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1132 474824 : return normalizepol_lg(x,l);
1133 : }
1134 :
1135 : /* z in Z[X], return z * Mod(1,p), normalized*/
1136 : GEN
1137 1669085 : FpX_to_mod(GEN z, GEN p)
1138 : {
1139 1669085 : long i, l = lg(z);
1140 : GEN x;
1141 1669085 : if (l == 2)
1142 : {
1143 3171 : x = cgetg(3,t_POL); x[1] = z[1];
1144 3171 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1145 : }
1146 1665914 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1147 18370815 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1148 1665914 : return normalizepol_lg(x,l);
1149 : }
1150 :
1151 : GEN
1152 84021 : FpXC_to_mod(GEN z, GEN p)
1153 : {
1154 84021 : long i,l = lg(z);
1155 84021 : GEN x = cgetg(l, t_COL);
1156 84021 : p = icopy(p);
1157 474166 : for (i=1; i<l; i++)
1158 390145 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1159 84021 : return x;
1160 : }
1161 :
1162 : static GEN
1163 0 : FpXC_to_mod_raw(GEN x, GEN p)
1164 0 : { pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1165 :
1166 : GEN
1167 0 : FpXM_to_mod(GEN z, GEN p)
1168 : {
1169 0 : long i,l = lg(z);
1170 0 : GEN x = cgetg(l, t_MAT);
1171 0 : p = icopy(p);
1172 0 : for (i=1; i<l; i++)
1173 0 : gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1174 0 : return x;
1175 : }
1176 :
1177 : /* z in Z^n, return z * Mod(1,p), normalized*/
1178 : GEN
1179 35984 : FpV_to_mod(GEN z, GEN p)
1180 : {
1181 35984 : long i,l = lg(z);
1182 35984 : GEN x = cgetg(l, t_VEC);
1183 35984 : if (l == 1) return x;
1184 35984 : p = icopy(p);
1185 72367 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1186 35984 : return x;
1187 : }
1188 : /* z in Z^n, return z * Mod(1,p), normalized*/
1189 : GEN
1190 105 : FpC_to_mod(GEN z, GEN p)
1191 : {
1192 105 : long i, l = lg(z);
1193 105 : GEN x = cgetg(l, t_COL);
1194 105 : if (l == 1) return x;
1195 91 : p = icopy(p);
1196 805 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1197 91 : return x;
1198 : }
1199 : /* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1200 : GEN
1201 226 : FpM_to_mod(GEN z, GEN p)
1202 : {
1203 226 : long i, j, m, l = lg(z);
1204 226 : GEN x = cgetg(l,t_MAT), y, zi;
1205 226 : if (l == 1) return x;
1206 205 : m = lgcols(z);
1207 205 : p = icopy(p);
1208 2456 : for (i=1; i<l; i++)
1209 : {
1210 2251 : gel(x,i) = cgetg(m,t_COL);
1211 2251 : y = gel(x,i); zi= gel(z,i);
1212 66799 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1213 : }
1214 205 : return x;
1215 : }
1216 : GEN
1217 28 : Flc_to_mod(GEN z, ulong pp)
1218 : {
1219 28 : long i, l = lg(z);
1220 28 : GEN p, x = cgetg(l, t_COL);
1221 28 : if (l == 1) return x;
1222 28 : p = utoipos(pp);
1223 112 : for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1224 28 : return x;
1225 : }
1226 : GEN
1227 59694 : Flm_to_mod(GEN z, ulong pp)
1228 : {
1229 59694 : long i, j, m, l = lg(z);
1230 59694 : GEN p, x = cgetg(l,t_MAT), y, zi;
1231 59694 : if (l == 1) return x;
1232 59666 : m = lgcols(z);
1233 59666 : p = utoipos(pp);
1234 201328 : for (i=1; i<l; i++)
1235 : {
1236 141662 : gel(x,i) = cgetg(m,t_COL);
1237 141662 : y = gel(x,i); zi= gel(z,i);
1238 720428 : for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1239 : }
1240 59666 : return x;
1241 : }
1242 :
1243 : GEN
1244 574 : FpVV_to_mod(GEN z, GEN p)
1245 : {
1246 574 : long i, j, m, l = lg(z);
1247 574 : GEN x = cgetg(l,t_VEC), y, zi;
1248 574 : if (l == 1) return x;
1249 574 : m = lgcols(z);
1250 574 : p = icopy(p);
1251 1246 : for (i=1; i<l; i++)
1252 : {
1253 672 : gel(x,i) = cgetg(m,t_VEC);
1254 672 : y = gel(x,i); zi= gel(z,i);
1255 2016 : for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1256 : }
1257 574 : return x;
1258 : }
1259 :
1260 : /* z in Z^n, return z * Mod(1,p), normalized*/
1261 : GEN
1262 7 : FpXQC_to_mod(GEN z, GEN T, GEN p)
1263 : {
1264 7 : long i,l = lg(z);
1265 7 : GEN x = cgetg(l, t_COL);
1266 7 : if (l == 1) return x;
1267 7 : p = icopy(p);
1268 7 : T = FpX_to_mod_raw(T, p);
1269 21 : for (i=1; i<l; i++)
1270 14 : gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1271 7 : return x;
1272 : }
1273 :
1274 : static GEN
1275 582533 : Fq_to_mod_raw(GEN x, GEN T, GEN p)
1276 : {
1277 582533 : GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1278 582533 : return mkpolmod(z, T);
1279 : }
1280 :
1281 : /* z in Z^n, return z * Mod(1,p), normalized*/
1282 : GEN
1283 28 : FqC_to_mod(GEN z, GEN T, GEN p)
1284 : {
1285 : GEN x;
1286 28 : long i,l = lg(z);
1287 28 : if (!T) return FpC_to_mod(z, p);
1288 28 : x = cgetg(l, t_COL);
1289 28 : if (l == 1) return x;
1290 28 : p = icopy(p);
1291 28 : T = FpX_to_mod_raw(T, p);
1292 504 : for (i=1; i<l; i++)
1293 476 : gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1294 28 : return x;
1295 : }
1296 :
1297 : /* z in Z^n, return z * Mod(1,p), normalized*/
1298 : static GEN
1299 107975 : FqC_to_mod_raw(GEN x, GEN T, GEN p)
1300 690032 : { pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1301 :
1302 : /* z in Z^n, return z * Mod(1,p), normalized*/
1303 : GEN
1304 26145 : FqM_to_mod(GEN z, GEN T, GEN p)
1305 : {
1306 : GEN x;
1307 26145 : long i,l = lg(z);
1308 26145 : if (!T) return FpM_to_mod(z, p);
1309 26145 : x = cgetg(l, t_MAT);
1310 26145 : if (l == 1) return x;
1311 26117 : p = icopy(p);
1312 26117 : T = FpX_to_mod_raw(T, p);
1313 134092 : for (i=1; i<l; i++)
1314 107975 : gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1315 26117 : return x;
1316 : }
1317 :
1318 : /********************************************************************/
1319 : /* */
1320 : /* Blackbox linear algebra */
1321 : /* */
1322 : /********************************************************************/
1323 :
1324 : /* A sparse column (zCs) v is a t_COL with two components C and E which are
1325 : * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1326 : * (e_j) is the canonical basis.
1327 : * A sparse matrix (zMs) is a t_VEC of zCs */
1328 :
1329 : /* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1330 : * integer representing an element of Fp. This is important since p can be
1331 : * large and p+E[i] would not fit in a C long. */
1332 :
1333 : /* RgCs and RgMs are similar, except that the type of the component is
1334 : * unspecified. Functions handling RgCs/RgMs must be independent of the type
1335 : * of E. */
1336 :
1337 : /* Most functions take an argument nbrow which is the number of lines of the
1338 : * column/matrix, which cannot be derived from the data. */
1339 :
1340 : GEN
1341 0 : zCs_to_ZC(GEN R, long nbrow)
1342 : {
1343 0 : GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1344 0 : long j, l = lg(C);
1345 0 : for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1346 0 : return c;
1347 : }
1348 :
1349 : GEN
1350 0 : zMs_to_ZM(GEN x, long nbrow)
1351 0 : { pari_APPLY_type(t_MAT, zCs_to_ZC(gel(x, i), nbrow)) }
1352 :
1353 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1354 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1355 : GEN
1356 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1357 : {
1358 147 : pari_sp ltop = avma;
1359 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1360 147 : if (ZV_equal0(B)) return zerocol(n);
1361 147 : while (++col <= n)
1362 : {
1363 147 : pari_sp btop = avma, av;
1364 : long i, lQ;
1365 147 : GEN V, Q, M, W = B;
1366 147 : GEN b = cgetg(m+2, t_POL);
1367 147 : b[1] = evalsigne(1)|evalvarn(0);
1368 147 : gel(b, 2) = gel(W, col);
1369 46219 : for (i = 3; i<m+2; i++)
1370 46072 : gel(b, i) = cgeti(lgefint(p));
1371 147 : av = avma;
1372 46219 : for (i = 3; i<m+2; i++)
1373 : {
1374 46072 : W = f(E, W);
1375 46072 : affii(gel(W, col),gel(b, i));
1376 46072 : if (gc_needed(av,1))
1377 : {
1378 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1379 72 : W = gerepileupto(av, W);
1380 : }
1381 : }
1382 147 : b = FpX_renormalize(b, m+2);
1383 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1384 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1385 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1386 147 : W = B; lQ =lg(Q);
1387 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1388 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1389 147 : av = avma;
1390 22742 : for (i = lQ-3; i > 1; i--)
1391 : {
1392 22595 : W = f(E, W);
1393 22595 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1394 22595 : if (gc_needed(av,1))
1395 : {
1396 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1397 110 : gerepileall(av, 2, &V, &W);
1398 : }
1399 : }
1400 147 : V = FpC_red(V, p);
1401 147 : W = FpC_sub(f(E,V), B, p);
1402 147 : if (ZV_equal0(W)) return gerepilecopy(ltop, V);
1403 0 : av = avma;
1404 0 : for (i = 1; i <= n; ++i)
1405 : {
1406 0 : V = W;
1407 0 : W = f(E, V);
1408 0 : if (ZV_equal0(W))
1409 0 : return gerepilecopy(ltop, shallowtrans(V));
1410 0 : gerepileall(av, 2, &V, &W);
1411 : }
1412 0 : set_avma(btop);
1413 : }
1414 0 : return NULL;
1415 : }
1416 :
1417 : GEN
1418 0 : zMs_ZC_mul(GEN M, GEN B)
1419 : {
1420 : long i, j;
1421 0 : long n = lg(B)-1;
1422 0 : GEN V = zerocol(n);
1423 0 : for (i = 1; i <= n; ++i)
1424 0 : if (signe(gel(B, i)))
1425 : {
1426 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1427 0 : long l = lg(C);
1428 0 : for (j = 1; j < l; ++j)
1429 : {
1430 0 : long k = C[j];
1431 0 : switch(E[j])
1432 : {
1433 0 : case 1:
1434 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1435 0 : break;
1436 0 : case -1:
1437 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1438 0 : break;
1439 0 : default:
1440 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]));
1441 0 : break;
1442 : }
1443 : }
1444 : }
1445 0 : return V;
1446 : }
1447 :
1448 : GEN
1449 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1450 :
1451 : GEN
1452 68961 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1453 : {
1454 68961 : long i, j, lM = lg(M);
1455 68961 : GEN V = cgetg(lM,t_VEC);
1456 28331175 : for (i = 1; i < lM; ++i)
1457 : {
1458 28262214 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1459 28262214 : pari_sp av = avma;
1460 28262214 : long lC = lg(C);
1461 28262214 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1462 28262214 : z = mulis(gel(B, C[1]), E[1]);
1463 228720437 : for (j = 2; j < lC; ++j)
1464 : {
1465 200458223 : GEN b = gel(B, C[j]);
1466 200458223 : switch(E[j])
1467 : {
1468 131566519 : case 1: z = addii(z, b); break;
1469 57521178 : case -1: z = subii(z, b); break;
1470 11370526 : default: z = addii(z, mulis(b, E[j])); break;
1471 : }
1472 : }
1473 28262214 : gel(V,i) = gerepileuptoint(av, p? Fp_red(z, p): z);
1474 : }
1475 68961 : return V;
1476 : }
1477 : GEN
1478 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1479 :
1480 : GEN
1481 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1482 : {
1483 0 : pari_sp av = avma, av2;
1484 0 : GEN xi, xb, pi = gen_1, P;
1485 : long i;
1486 0 : if (!C) {
1487 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1488 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1489 : }
1490 0 : P = utoipos(p);
1491 0 : av2 = avma;
1492 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1493 0 : xb = Flm_to_ZM(xi);
1494 0 : for (i = 2; i <= e; i++)
1495 : {
1496 0 : pi = muliu(pi, p); /* = p^(i-1) */
1497 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1498 0 : if (gc_needed(av,2))
1499 : {
1500 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1501 0 : gerepileall(av2,3, &pi,&b,&xb);
1502 : }
1503 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1504 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1505 : }
1506 0 : return gerepileupto(av, xb);
1507 : }
1508 :
1509 : struct wrapper_modp_s {
1510 : GEN (*f)(void*, GEN);
1511 : void *E;
1512 : GEN p;
1513 : };
1514 :
1515 : /* compute f(x) mod p */
1516 : static GEN
1517 0 : wrap_relcomb_modp(void *E, GEN x)
1518 : {
1519 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1520 0 : return FpC_red(W->f(W->E, x), W->p);
1521 : }
1522 : static GEN
1523 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1524 :
1525 : static GEN
1526 68814 : wrap_relker(void*E, GEN x)
1527 : {
1528 68814 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1529 68814 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1530 : }
1531 :
1532 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1533 : GEN
1534 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1535 : {
1536 : struct wrapper_modp_s W;
1537 0 : pari_sp av = avma;
1538 0 : GEN xb, xi, pi = gen_1;
1539 : long i;
1540 0 : W.E = E;
1541 0 : W.f = f;
1542 0 : W.p = p;
1543 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1544 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1545 0 : xb = xi;
1546 0 : for (i = 2; i <= e; i++)
1547 : {
1548 0 : pi = mulii(pi, p); /* = p^(i-1) */
1549 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1550 0 : if (gc_needed(av,2))
1551 : {
1552 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1553 0 : gerepileall(av,3, &pi,&B,&xb);
1554 : }
1555 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1556 0 : if (!xi) return NULL;
1557 0 : if (typ(xi) == t_VEC) return gerepileupto(av, xi);
1558 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1559 : }
1560 0 : return gerepileupto(av, xb);
1561 : }
1562 :
1563 : static GEN
1564 23036 : vecprow(GEN A, GEN prow)
1565 : {
1566 23036 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1567 : }
1568 :
1569 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1570 : * or the index of a column which is linearly dependent from the others as a
1571 : * t_VECSMALL with a single component. */
1572 : GEN
1573 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1574 : {
1575 0 : pari_sp av = avma;
1576 : GEN pcol, prow;
1577 0 : long nbi=lg(M)-1, lR;
1578 : long i, n;
1579 : GEN Mp, Ap, Rp;
1580 : pari_timer ti;
1581 0 : if (DEBUGLEVEL) timer_start(&ti);
1582 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1583 0 : if (!pcol) return gc_NULL(av);
1584 0 : if (DEBUGLEVEL)
1585 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1586 0 : n = lg(pcol)-1;
1587 0 : Mp = cgetg(n+1, t_MAT);
1588 0 : for(i=1; i<=n; i++)
1589 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1590 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1591 0 : if (DEBUGLEVEL) timer_start(&ti);
1592 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1593 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1594 0 : if (!Rp) return gc_NULL(av);
1595 0 : lR = lg(Rp)-1;
1596 0 : if (typ(Rp) == t_COL)
1597 : {
1598 0 : GEN R = zerocol(nbi+1);
1599 0 : for (i=1; i<=lR; i++)
1600 0 : gel(R,pcol[i]) = gel(Rp,i);
1601 0 : return gerepilecopy(av, R);
1602 : }
1603 0 : for (i = 1; i <= lR; ++i)
1604 0 : if (signe(gel(Rp, i)))
1605 0 : return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
1606 : return NULL; /* LCOV_EXCL_LINE */
1607 : }
1608 :
1609 : GEN
1610 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1611 : {
1612 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1613 : }
1614 :
1615 : GEN
1616 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1617 : {
1618 : GEN res;
1619 0 : pari_CATCH(e_INV)
1620 : {
1621 0 : GEN E = pari_err_last();
1622 0 : GEN x = err_get_compo(E,2);
1623 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1624 0 : if (DEBUGLEVEL)
1625 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1626 0 : res = NULL;
1627 : } pari_TRY {
1628 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1629 0 : } pari_ENDCATCH
1630 0 : return res;
1631 : }
1632 :
1633 : static GEN
1634 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1635 : {
1636 147 : long i, j, oldf = 0, f = 0;
1637 147 : long lrow = lg(prow), lM = lg(M);
1638 147 : GEN W = const_vecsmall(lM-1,1);
1639 147 : GEN R = cgetg(lrow, t_VEC);
1640 228354 : for (i=1; i<lrow; i++)
1641 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1642 : do
1643 : {
1644 442 : oldf = f;
1645 374995 : for(i=1; i<lM; i++)
1646 374553 : if (W[i])
1647 : {
1648 131789 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1649 131789 : long c=0, cj=0, lF = lg(F);
1650 1094147 : for(j=1; j<lF; j++)
1651 962358 : if (!gel(R,F[j])) { c++; cj=j; }
1652 131789 : if (c>=2) continue;
1653 112311 : if (c==1)
1654 : {
1655 32898 : pari_sp av = avma;
1656 32898 : GEN s = gen_0;
1657 272977 : for(j=1; j<lF; j++)
1658 240079 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1659 : /* s /= -E[cj] mod p */
1660 32898 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1661 32898 : gel(R,F[cj]) = gerepileuptoint(av, s);
1662 32898 : f++;
1663 : }
1664 112311 : W[i]=0;
1665 : }
1666 442 : } while(oldf!=f);
1667 228354 : for (i=1; i<lrow; i++)
1668 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1669 147 : return R;
1670 : }
1671 :
1672 : /* Return a linear form Y such that YM is essentially 0 */
1673 : GEN
1674 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1675 : {
1676 147 : pari_sp av = avma, av2;
1677 : GEN Mp, pcol, prow;
1678 : long i, n;
1679 : pari_timer ti;
1680 : struct wrapper_modp_s W;
1681 147 : if (DEBUGLEVEL) timer_start(&ti);
1682 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1683 147 : if (!pcol) return gc_NULL(av);
1684 147 : if (DEBUGLEVEL)
1685 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1686 147 : n = lg(pcol)-1;
1687 147 : Mp = cgetg(n+1, t_MAT);
1688 23183 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1689 147 : W.E = (void*)Mp;
1690 147 : W.p = p;
1691 147 : av2 = avma;
1692 0 : for(;; set_avma(av2))
1693 0 : {
1694 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1695 147 : if (DEBUGLEVEL) timer_start(&ti);
1696 147 : pari_CATCH(e_INV)
1697 : {
1698 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1699 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1700 0 : if (DEBUGLEVEL)
1701 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1702 0 : Rp = NULL;
1703 : } pari_TRY {
1704 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1705 147 : } pari_ENDCATCH
1706 147 : if (!Rp) continue;
1707 147 : if (typ(Rp)==t_COL)
1708 : {
1709 147 : Rp = FpC_sub(Rp,B,p);
1710 147 : if (ZV_equal0(Rp)) continue;
1711 : }
1712 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1713 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1714 147 : return gerepilecopy(av, R);
1715 : }
1716 : }
1717 :
1718 : GEN
1719 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1720 : {
1721 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1722 : }
|
