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