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 31527474 : FpC_red(GEN x, GEN p)
28 276874269 : { 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 7061870 : FpC_center(GEN x, GEN p, GEN pov2)
37 43973121 : { 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 5655485 : FpM_red(GEN x, GEN p)
46 25921550 : { pari_APPLY_same(FpC_red(gel(x,i), p)) }
47 :
48 : GEN
49 2809877 : FpM_center(GEN x, GEN p, GEN pov2)
50 8440200 : { 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 23186 : 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 41416 : random_Flv(long n, ulong p)
127 : {
128 41416 : GEN y = cgetg(n+1, t_VECSMALL);
129 : long i;
130 262158 : for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
131 41417 : return y;
132 : }
133 :
134 : /********************************************************************/
135 : /** **/
136 : /** ADD, SUB **/
137 : /** **/
138 : /********************************************************************/
139 : GEN
140 283601 : FpC_add(GEN x, GEN y, GEN p)
141 3551393 : { pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
142 :
143 : GEN
144 0 : FpV_add(GEN x, GEN y, GEN p)
145 0 : { pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
146 :
147 : GEN
148 0 : FpM_add(GEN x, GEN y, GEN p)
149 0 : { pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
150 :
151 : GEN
152 1812496 : Flv_add(GEN x, GEN y, ulong p)
153 : {
154 1812496 : if (p==2)
155 7096847 : pari_APPLY_ulong(x[i]^y[i])
156 : else
157 41999014 : pari_APPLY_ulong(Fl_add(x[i], y[i], p))
158 : }
159 :
160 : void
161 509595 : Flv_add_inplace(GEN x, GEN y, ulong p)
162 : {
163 509595 : long i, l = lg(x);
164 509595 : if (p==2)
165 1570127 : 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 509595 : }
169 :
170 : ulong
171 3871 : Flv_sum(GEN x, ulong p)
172 : {
173 3871 : long i, l = lg(x);
174 3871 : ulong s = 0;
175 3871 : if (p==2)
176 17689 : for (i = 1; i < l; i++) s ^= x[i];
177 : else
178 0 : for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
179 3871 : return s;
180 : }
181 :
182 : GEN
183 965634 : FpC_sub(GEN x, GEN y, GEN p)
184 18959601 : { 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 236788299 : Flv_sub(GEN x, GEN y, ulong p)
196 1057815286 : { 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 51196 : Flm_Fl_add(GEN x, ulong y, ulong p)
207 : {
208 51196 : long l = lg(x), i, j;
209 51196 : GEN z = cgetg(l,t_MAT);
210 :
211 51196 : if (l==1) return z;
212 51196 : if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
213 241315 : for (i=1; i<l; i++)
214 : {
215 190119 : GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
216 190119 : gel(z,i) = zi;
217 1331276 : for (j=1; j<l; j++) zi[j] = xi[j];
218 190119 : zi[i] = Fl_add(zi[i], y, p);
219 : }
220 51196 : return z;
221 : }
222 : GEN
223 3724 : Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
224 :
225 : GEN
226 29322 : Flm_add(GEN x, GEN y, ulong p)
227 695951 : { pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
228 :
229 : GEN
230 26078373 : Flm_sub(GEN x, GEN y, ulong p)
231 262839986 : { pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
232 :
233 : /********************************************************************/
234 : /** **/
235 : /** MULTIPLICATION **/
236 : /** **/
237 : /********************************************************************/
238 : GEN
239 973022 : FpC_Fp_mul(GEN x, GEN y, GEN p)
240 11762513 : { pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
241 :
242 : GEN
243 1754386 : Flv_Fl_mul(GEN x, ulong y, ulong p)
244 41056983 : { 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 766408 : FpM_Fp_mul(GEN X, GEN c, GEN p) {
258 766408 : long i, j, h, l = lg(X);
259 766408 : GEN A = cgetg(l, t_MAT);
260 766409 : if (l == 1) return A;
261 766409 : h = lgcols(X);
262 4083046 : for (j=1; j<l; j++)
263 : {
264 3316768 : GEN a = cgetg(h, t_COL), x = gel(X, j);
265 24459804 : for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
266 3316635 : gel(A,j) = a;
267 : }
268 766278 : return A;
269 : }
270 :
271 : /* x *= y */
272 : void
273 7910303 : Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
274 : {
275 : long i;
276 7910303 : if (HIGHWORD(y | p))
277 8606023 : for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
278 : else
279 32281641 : for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
280 7910302 : }
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 17529503 : Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
303 : {
304 17529503 : long i, j, m, l = lg(y);
305 17529503 : GEN z = cgetg(l, t_MAT);
306 17526967 : if (l == 1) return z;
307 17526967 : m = lgcols(y);
308 159580597 : for(j=1; j<l; j++) {
309 141942098 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
310 381637921 : for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
311 : }
312 17638499 : return z;
313 : }
314 :
315 : /* return x * y */
316 : static GEN
317 8647709 : Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
318 : {
319 8647709 : long i, j, m, l = lg(y);
320 8647709 : GEN z = cgetg(l, t_MAT);
321 8647706 : if (l == 1) return z;
322 8647706 : m = lgcols(y);
323 53947546 : for(j=1; j<l; j++) {
324 45299844 : GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
325 163523483 : for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
326 : }
327 8647702 : return z;
328 : }
329 :
330 : /* return x * y */
331 : GEN
332 26094686 : Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
333 : {
334 26094686 : if (HIGHWORD(p))
335 17529588 : return Flm_Fl_mul_pre_i(y, x, p, pi);
336 : else
337 8565098 : return Flm_Fl_mul_OK(y, x, p);
338 : }
339 :
340 : /* return x * y */
341 : GEN
342 82649 : Flm_Fl_mul(GEN y, ulong x, ulong p)
343 : {
344 82649 : if (HIGHWORD(p))
345 74 : return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
346 : else
347 82575 : return Flm_Fl_mul_OK(y, x, p);
348 : }
349 :
350 : GEN
351 4102742 : Flv_neg(GEN x, ulong p)
352 31302875 : { pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
353 :
354 : void
355 6585 : Flv_neg_inplace(GEN v, ulong p)
356 : {
357 : long i;
358 192319 : for (i = 1; i < lg(v); ++i)
359 185734 : v[i] = Fl_neg(v[i], p);
360 6585 : }
361 :
362 : GEN
363 324988 : Flm_neg(GEN x, ulong p)
364 4396985 : { 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 18945987 : FpMrow_FpC_mul_i(GEN x, GEN y, long lx, long i, GEN p)
369 : {
370 18945987 : pari_sp av = avma;
371 18945987 : GEN c = mulii(gcoeff(x,i,1), gel(y,1));
372 : long k;
373 230157469 : for (k = 2; k < lx; k++)
374 : {
375 211213484 : GEN t = mulii(gcoeff(x,i,k), gel(y,k));
376 211208725 : if (signe(t)) c = addii(c, t);
377 : }
378 18943985 : 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 771340447 : Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
422 : {
423 771340447 : ulong c = ucoeff(x,i,1) * uel(y,1);
424 : long k;
425 12024864557 : for (k = 2; k < lx; k++) {
426 11253524110 : c += ucoeff(x,i,k) * uel(y,k);
427 11253524110 : if (c & HIGHBIT) c %= p;
428 : }
429 771340447 : return c % p;
430 : }
431 :
432 : static ulong
433 548758756 : 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 548758756 : long k = 1;
437 : LOCAL_OVERFLOW;
438 : LOCAL_HIREMAINDER;
439 548758756 : l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
440 8975696752 : while (++k < lx) {
441 8426937996 : l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
442 8426937996 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
443 : }
444 548758756 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
445 53809381 : else return remlll_pre(v1, h1, l1, p, pi);
446 : }
447 :
448 : static GEN
449 339422 : Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
450 : {
451 : long i,j;
452 339422 : GEN z = NULL;
453 :
454 3753581 : for (j=1; j<lx; j++)
455 : {
456 3414159 : if (!y[j]) continue;
457 1121700 : if (!z) z = Flv_copy(gel(x,j));
458 29830274 : else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
459 : }
460 339422 : if (!z) z = zero_zv(l-1);
461 339422 : return z;
462 : }
463 :
464 : static GEN
465 1762539 : FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
466 : {
467 1762539 : GEN z = cgetg(l,t_COL);
468 : long i;
469 16816788 : for (i = 1; i < l; i++) gel(z,i) = FpMrow_FpC_mul_i(x, y, lx, i, p);
470 1762366 : return z;
471 : }
472 :
473 : static void
474 86833457 : __Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
475 : {
476 : long i;
477 858182493 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
478 86846829 : }
479 : static GEN
480 83611174 : Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
481 : {
482 83611174 : GEN z = cgetg(l,t_VECSMALL);
483 83610347 : __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
484 83611074 : return z;
485 : }
486 :
487 : static void
488 98247603 : __Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
489 : {
490 : long i;
491 649640184 : for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
492 98610958 : }
493 : static GEN
494 93885063 : Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
495 : {
496 93885063 : GEN z = cgetg(l,t_VECSMALL);
497 93733869 : __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
498 93944214 : return z;
499 : }
500 :
501 : GEN
502 2062861 : FpM_mul(GEN x, GEN y, GEN p)
503 : {
504 2062861 : long lx=lg(x), ly=lg(y);
505 : GEN z;
506 2062861 : pari_sp av = avma;
507 2062861 : if (ly==1) return cgetg(1,t_MAT);
508 2062861 : if (lx==1) return zeromat(0, ly-1);
509 2062861 : if (lgefint(p) == 3)
510 : {
511 1886039 : ulong pp = uel(p,2);
512 1886039 : if (pp == 2)
513 : {
514 724114 : x = ZM_to_F2m(x);
515 724164 : y = ZM_to_F2m(y);
516 724170 : z = F2m_to_ZM(F2m_mul(x,y));
517 : }
518 : else
519 : {
520 1161925 : x = ZM_to_Flm(x, pp);
521 1161948 : y = ZM_to_Flm(y, pp);
522 1161950 : z = Flm_to_ZM(Flm_mul(x,y, pp));
523 : }
524 : } else
525 176822 : z = FpM_red(ZM_mul(x, y), p);
526 2062944 : return gc_upto(av, z);
527 : }
528 :
529 : static GEN
530 28651564 : Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
531 : {
532 : long j;
533 28651564 : GEN z = cgetg(ly, t_MAT);
534 28650358 : if (SMALL_ULONG(p))
535 103553707 : for (j=1; j<ly; j++)
536 83101241 : gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
537 : else
538 101881443 : for (j=1; j<ly; j++)
539 93684599 : gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
540 28649310 : 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 2782358 : for (j = 1; j <= min_e; j++) {
555 2716054 : C = gel(M, j + dx) + dy;
556 123912365 : for (i = 1; i <= min_d; i++)
557 121196311 : uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
558 121195522 : ucoeff(B, mb + i, nb + j), p);
559 2781403 : for (; i <= da; i++)
560 64560 : uel(C, i) = ucoeff(A, ma + i, na + j);
561 2716843 : for (; i <= db; i++)
562 0 : uel(C, i) = ucoeff(B, mb + i, nb + j);
563 2716843 : 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 1415648 : for (j = 1; j <= n; j++)
596 1382481 : gel(M, j) = cgetg(m + 1, t_VECSMALL);
597 33167 : 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 58014 : 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 58014 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
609 58014 : GEN M = cgetg(n + 1, t_MAT), C;
610 :
611 2527718 : for (j = 1; j <= min_e; j++) {
612 2469703 : gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
613 116106937 : for (i = 1; i <= min_d; i++)
614 113637235 : uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
615 113637440 : ucoeff(B, mb + i, nb + j), p);
616 2569943 : for (; i <= da; i++)
617 100446 : uel(C, i) = ucoeff(A, ma + i, na + j);
618 2544541 : for (; i <= db; i++)
619 75044 : uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
620 2469497 : for (; i <= m; i++)
621 0 : uel(C, i) = 0;
622 : }
623 58015 : 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 59434 : 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 59434 : for (; j <= n; j++)
638 1419 : gel(M, j) = zero_Flv(m);
639 58015 : 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 28658731 : Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
723 : {
724 28658731 : ulong e = expu(p);
725 : #ifdef LONG_IS_64BIT /* Beware to update ZM_mul_i if this changes */
726 24066309 : long Flm_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
727 : #else
728 4592902 : long Flm_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
729 : #endif
730 28659211 : if (l <= Flm_sw_bound || lx <= Flm_sw_bound || ly <= Flm_sw_bound)
731 28650923 : 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 13950181 : Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
738 : {
739 13950181 : long lx=lg(x), ly=lg(y);
740 13950181 : if (ly==1) return cgetg(1,t_MAT);
741 13950181 : if (lx==1) return zero_Flm(0, ly-1);
742 13950181 : return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
743 : }
744 :
745 : GEN
746 14556865 : Flm_mul(GEN x, GEN y, ulong p)
747 : {
748 14556865 : long lx=lg(x), ly=lg(y);
749 14556865 : if (ly==1) return cgetg(1,t_MAT);
750 14556578 : if (lx==1) return zero_Flm(0, ly-1);
751 14556578 : 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 847 : _Flm_one(void *E)
776 847 : { 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 847 : Flm_powers(GEN x, ulong n, ulong p)
787 : {
788 : struct _Flm d;
789 847 : d.p = p;
790 847 : d.n = lg(x)-1;
791 847 : 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 5875 : FpC_FpV_mul(GEN x, GEN y, GEN p)
840 : {
841 5875 : long i,j, lx=lg(x), ly=lg(y);
842 : GEN z;
843 5875 : if (ly==1) return cgetg(1,t_MAT);
844 5875 : z = cgetg(ly,t_MAT);
845 64914 : for (j=1; j < ly; j++)
846 : {
847 59039 : GEN zj = cgetg(lx,t_COL);
848 1383646 : for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
849 59039 : gel(z, j) = zj;
850 : }
851 5875 : return z;
852 : }
853 :
854 : /* Multiply a line vector by a column and return a scalar (t_INT) */
855 : GEN
856 8684238 : FpV_dotproduct(GEN x, GEN y, GEN p)
857 : {
858 8684238 : long i, lx = lg(x);
859 : pari_sp av;
860 : GEN c;
861 8684238 : if (lx == 1) return gen_0;
862 8681725 : av = avma; c = mulii(gel(x,1),gel(y,1));
863 36971090 : for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
864 8681706 : 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 9624165 : Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
880 : {
881 9624165 : ulong c = uel(x,0) * uel(y,0);
882 : long k;
883 72860022 : for (k = 1; k < lx; k++) {
884 63235857 : c += uel(x,k) * uel(y,k);
885 63235857 : if (c & HIGHBIT) c %= p;
886 : }
887 9624165 : return c % p;
888 : }
889 :
890 : INLINE ulong
891 741860 : Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
892 : {
893 : ulong l0, l1, v1, h0, h1;
894 741860 : long i = 0;
895 : LOCAL_OVERFLOW;
896 : LOCAL_HIREMAINDER;
897 741860 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
898 15430105 : while (++i < lx) {
899 14688245 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
900 14688245 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
901 : }
902 741860 : if (v1 == 0) return remll_pre(h1, l1, p, pi);
903 464137 : else return remlll_pre(v1, h1, l1, p, pi);
904 : }
905 :
906 : ulong
907 730736 : Flv_dotproduct(GEN x, GEN y, ulong p)
908 : {
909 730736 : long lx = lg(x)-1;
910 730736 : if (lx == 0) return 0;
911 730736 : if (SMALL_ULONG(p))
912 730736 : 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 66256 : Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
919 : {
920 66256 : long lx = lg(x)-1;
921 66256 : if (lx == 0) return 0;
922 66256 : if (SMALL_ULONG(p))
923 60044 : return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
924 : else
925 6212 : return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
926 : }
927 :
928 : ulong
929 9638598 : Flx_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
930 : {
931 9638598 : long lx = minss(lgpol(x), lgpol(y));
932 9638651 : if (lx == 0) return 0;
933 9569026 : if (pi)
934 735648 : return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
935 : else
936 8833378 : 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 1762539 : FpM_FpC_mul(GEN x, GEN y, GEN p)
944 : {
945 1762539 : long lx = lg(x);
946 1762539 : return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
947 : }
948 : GEN
949 1086777 : Flm_Flc_mul(GEN x, GEN y, ulong p)
950 : {
951 1086777 : long l, lx = lg(x);
952 1086777 : if (lx==1) return cgetg(1,t_VECSMALL);
953 1086777 : l = lgcols(x);
954 1086777 : if (p==2)
955 339422 : return Flm_Flc_mul_i_2(x, y, lx, l);
956 747355 : else if (SMALL_ULONG(p))
957 510015 : 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 7729457 : Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
981 : {
982 7729457 : long l, lx = lg(x);
983 : GEN z;
984 7729457 : if (lx==1) return pol0_Flx(sv);
985 7729457 : l = lgcols(x);
986 7728002 : z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
987 7727488 : if (SMALL_ULONG(p))
988 3217869 : __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
989 : else
990 4509619 : __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
991 7741347 : 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 1234390 : FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
997 : {
998 1234390 : long i, l, lx = lg(x);
999 : GEN z;
1000 1234390 : if (lx==1) return pol_0(v);
1001 1234390 : l = lgcols(x);
1002 1234390 : z = new_chunk(l+1);
1003 1688751 : for (i=l-1; i; i--)
1004 : {
1005 1631765 : pari_sp av = avma;
1006 1631765 : GEN t = FpMrow_FpC_mul_i(x,y,lx,i,p);
1007 1631763 : if (signe(t))
1008 : {
1009 1177403 : if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
1010 1177404 : gel(z,i+1) = t;
1011 1177404 : break;
1012 : }
1013 454360 : set_avma(av);
1014 : }
1015 1234390 : if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
1016 1177404 : z[0] = evaltyp(t_POL) | _evallg(i+2);
1017 1177404 : z[1] = evalsigne(1) | evalvarn(v);
1018 3437149 : for (; i; i--) gel(z,i+1) = FpMrow_FpC_mul_i(x,y,lx,i,p);
1019 1177404 : return z;
1020 : }
1021 :
1022 : /********************************************************************/
1023 : /** **/
1024 : /** TRANSPOSITION **/
1025 : /** **/
1026 : /********************************************************************/
1027 :
1028 : /* == zm_transpose */
1029 : GEN
1030 735509 : Flm_transpose(GEN x)
1031 : {
1032 735509 : long i, dx, lx = lg(x);
1033 : GEN y;
1034 735509 : if (lx == 1) return cgetg(1,t_MAT);
1035 735376 : dx = lgcols(x); y = cgetg(dx,t_MAT);
1036 9446689 : for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1037 735379 : 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 1440810 : matid_Flm(long n)
1079 : {
1080 1440810 : GEN y = cgetg(n+1,t_MAT);
1081 : long i;
1082 1440802 : if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1083 9212553 : for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1084 1440799 : 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 58663280 : ZV_to_Flv(GEN x, ulong p)
1103 727656922 : { pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1104 :
1105 : GEN
1106 13779234 : ZM_to_Flm(GEN x, ulong p)
1107 71214461 : { 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 18978300 : 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 4007 : Fp_to_mod(GEN z, GEN p)
1125 : {
1126 4007 : retmkintmod(modii(z, p), icopy(p));
1127 : }
1128 :
1129 : GEN
1130 997269 : FpX_to_mod_raw(GEN z, GEN p)
1131 : {
1132 997269 : long i, l = lg(z);
1133 : GEN x;
1134 997269 : 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 474845 : x = cgetg(l,t_POL); x[1] = z[1];
1140 3918691 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1141 474845 : return normalizepol_lg(x,l);
1142 : }
1143 :
1144 : /* z in Z[X], return z * Mod(1,p), normalized*/
1145 : GEN
1146 1260826 : FpX_to_mod(GEN z, GEN p)
1147 : {
1148 1260826 : long i, l = lg(z);
1149 : GEN x;
1150 1260826 : if (l == 2)
1151 : {
1152 1197 : x = cgetg(3,t_POL); x[1] = z[1];
1153 1197 : gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1154 : }
1155 1259629 : x = cgetg(l,t_POL); x[1] = z[1]; p = icopy(p);
1156 16690387 : for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1157 1259629 : return normalizepol_lg(x,l);
1158 : }
1159 :
1160 : GEN
1161 84021 : FpXC_to_mod(GEN z, GEN p)
1162 : {
1163 84021 : long i,l = lg(z);
1164 84021 : GEN x = cgetg(l, t_COL);
1165 84021 : p = icopy(p);
1166 474159 : for (i=1; i<l; i++)
1167 390138 : gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1168 84021 : 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 36495 : FpV_to_mod(GEN z, GEN p)
1189 : {
1190 36495 : long i,l = lg(z);
1191 36495 : GEN x = cgetg(l, t_VEC);
1192 36495 : if (l == 1) return x;
1193 36495 : p = icopy(p);
1194 73585 : for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1195 36495 : 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 : /* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1363 : * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1364 : GEN
1365 147 : gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1366 : {
1367 147 : pari_sp ltop = avma;
1368 147 : long col = 0, n = lg(B)-1, m = 2*n+1;
1369 147 : if (ZV_equal0(B)) return zerocol(n);
1370 147 : while (++col <= n)
1371 : {
1372 147 : pari_sp btop = avma, av;
1373 : long i, lQ;
1374 147 : GEN V, Q, M, W = B;
1375 147 : GEN b = cgetg(m+2, t_POL);
1376 147 : b[1] = evalsigne(1)|evalvarn(0);
1377 147 : gel(b, 2) = gel(W, col);
1378 46225 : for (i = 3; i<m+2; i++)
1379 46078 : gel(b, i) = cgeti(lgefint(p));
1380 147 : av = avma;
1381 46225 : for (i = 3; i<m+2; i++)
1382 : {
1383 46078 : W = f(E, W);
1384 46078 : affii(gel(W, col),gel(b, i));
1385 46078 : if (gc_needed(av,1))
1386 : {
1387 72 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1388 72 : W = gc_upto(av, W);
1389 : }
1390 : }
1391 147 : b = FpX_renormalize(b, m+2);
1392 147 : if (lgpol(b)==0) {set_avma(btop); continue; }
1393 147 : M = FpX_halfgcd(b, pol_xn(m, 0), p);
1394 147 : Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1395 147 : W = B; lQ =lg(Q);
1396 147 : if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1397 147 : V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1398 147 : av = avma;
1399 22745 : for (i = lQ-3; i > 1; i--)
1400 : {
1401 22598 : W = f(E, W);
1402 22598 : V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1403 22598 : if (gc_needed(av,1))
1404 : {
1405 110 : if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1406 110 : (void)gc_all(av, 2, &V, &W);
1407 : }
1408 : }
1409 147 : V = FpC_red(V, p);
1410 147 : W = FpC_sub(f(E,V), B, p);
1411 147 : if (ZV_equal0(W)) return gc_GEN(ltop, V);
1412 0 : av = avma;
1413 0 : for (i = 1; i <= n; ++i)
1414 : {
1415 0 : V = W;
1416 0 : W = f(E, V);
1417 0 : if (ZV_equal0(W))
1418 0 : return gc_GEN(ltop, shallowtrans(V));
1419 0 : (void)gc_all(av, 2, &V, &W);
1420 : }
1421 0 : set_avma(btop);
1422 : }
1423 0 : return NULL;
1424 : }
1425 :
1426 : GEN
1427 0 : zMs_ZC_mul(GEN M, GEN B)
1428 : {
1429 : long i, j;
1430 0 : long n = lg(B)-1;
1431 0 : GEN V = zerocol(n);
1432 0 : for (i = 1; i <= n; ++i)
1433 0 : if (signe(gel(B, i)))
1434 : {
1435 0 : GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1436 0 : long l = lg(C);
1437 0 : for (j = 1; j < l; ++j)
1438 : {
1439 0 : long k = C[j];
1440 0 : switch(E[j])
1441 : {
1442 0 : case 1:
1443 0 : gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1444 0 : break;
1445 0 : case -1:
1446 0 : gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1447 0 : break;
1448 0 : default:
1449 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]));
1450 0 : break;
1451 : }
1452 : }
1453 : }
1454 0 : return V;
1455 : }
1456 :
1457 : GEN
1458 0 : FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1459 :
1460 : GEN
1461 68970 : FpV_FpMs_mul(GEN B, GEN M, GEN p)
1462 : {
1463 68970 : long i, j, lM = lg(M);
1464 68970 : GEN V = cgetg(lM,t_VEC);
1465 28333281 : for (i = 1; i < lM; ++i)
1466 : {
1467 28264311 : GEN z, R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1468 28264311 : pari_sp av = avma;
1469 28264311 : long lC = lg(C);
1470 28264311 : if (lC == 1) { gel(V,i) = gen_0; continue; }
1471 28264311 : z = mulis(gel(B, C[1]), E[1]);
1472 228739358 : for (j = 2; j < lC; ++j)
1473 : {
1474 200475047 : GEN b = gel(B, C[j]);
1475 200475047 : switch(E[j])
1476 : {
1477 131566519 : case 1: z = addii(z, b); break;
1478 57528522 : case -1: z = subii(z, b); break;
1479 11380006 : default: z = addii(z, mulis(b, E[j])); break;
1480 : }
1481 : }
1482 28264311 : gel(V,i) = gc_INT(av, p? Fp_red(z, p): z);
1483 : }
1484 68970 : return V;
1485 : }
1486 : GEN
1487 0 : ZV_zMs_mul(GEN B, GEN M) { return FpV_FpMs_mul(B, M, NULL); }
1488 :
1489 : GEN
1490 0 : ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1491 : {
1492 0 : pari_sp av = avma, av2;
1493 0 : GEN xi, xb, pi = gen_1, P;
1494 : long i;
1495 0 : if (!C) {
1496 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
1497 0 : if (!C) pari_err_INV("ZlM_gauss", a);
1498 : }
1499 0 : P = utoipos(p);
1500 0 : av2 = avma;
1501 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1502 0 : xb = Flm_to_ZM(xi);
1503 0 : for (i = 2; i <= e; i++)
1504 : {
1505 0 : pi = muliu(pi, p); /* = p^(i-1) */
1506 0 : b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1507 0 : if (gc_needed(av,2))
1508 : {
1509 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1510 0 : (void)gc_all(av2,3, &pi,&b,&xb);
1511 : }
1512 0 : xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1513 0 : xb = ZM_add(xb, nm_Z_mul(xi, pi));
1514 : }
1515 0 : return gc_upto(av, xb);
1516 : }
1517 :
1518 : struct wrapper_modp_s {
1519 : GEN (*f)(void*, GEN);
1520 : void *E;
1521 : GEN p;
1522 : };
1523 :
1524 : /* compute f(x) mod p */
1525 : static GEN
1526 0 : wrap_relcomb_modp(void *E, GEN x)
1527 : {
1528 0 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1529 0 : return FpC_red(W->f(W->E, x), W->p);
1530 : }
1531 : static GEN
1532 0 : wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1533 :
1534 : static GEN
1535 68823 : wrap_relker(void*E, GEN x)
1536 : {
1537 68823 : struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1538 68823 : return FpV_FpMs_mul(x, (GEN) W->E, W->p);
1539 : }
1540 :
1541 : /* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1542 : GEN
1543 0 : gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1544 : {
1545 : struct wrapper_modp_s W;
1546 0 : pari_sp av = avma;
1547 0 : GEN xb, xi, pi = gen_1;
1548 : long i;
1549 0 : W.E = E;
1550 0 : W.f = f;
1551 0 : W.p = p;
1552 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1553 0 : if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1554 0 : xb = xi;
1555 0 : for (i = 2; i <= e; i++)
1556 : {
1557 0 : pi = mulii(pi, p); /* = p^(i-1) */
1558 0 : B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1559 0 : if (gc_needed(av,2))
1560 : {
1561 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1562 0 : (void)gc_all(av,3, &pi,&B,&xb);
1563 : }
1564 0 : xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1565 0 : if (!xi) return NULL;
1566 0 : if (typ(xi) == t_VEC) return gc_upto(av, xi);
1567 0 : xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1568 : }
1569 0 : return gc_upto(av, xb);
1570 : }
1571 :
1572 : static GEN
1573 23039 : vecprow(GEN A, GEN prow)
1574 : {
1575 23039 : return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1576 : }
1577 :
1578 : /* Solve the equation MX = A. Return either a solution as a t_COL,
1579 : * or the index of a column which is linearly dependent from the others as a
1580 : * t_VECSMALL with a single component. */
1581 : GEN
1582 0 : ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1583 : {
1584 0 : pari_sp av = avma;
1585 : GEN pcol, prow;
1586 0 : long nbi=lg(M)-1, lR;
1587 : long i, n;
1588 : GEN Mp, Ap, Rp;
1589 : pari_timer ti;
1590 0 : if (DEBUGLEVEL) timer_start(&ti);
1591 0 : RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1592 0 : if (!pcol) return gc_NULL(av);
1593 0 : if (DEBUGLEVEL)
1594 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1595 0 : n = lg(pcol)-1;
1596 0 : Mp = cgetg(n+1, t_MAT);
1597 0 : for(i=1; i<=n; i++)
1598 0 : gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1599 0 : Ap = zCs_to_ZC(vecprow(A, prow), n);
1600 0 : if (DEBUGLEVEL) timer_start(&ti);
1601 0 : Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1602 0 : if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1603 0 : if (!Rp) return gc_NULL(av);
1604 0 : lR = lg(Rp)-1;
1605 0 : if (typ(Rp) == t_COL)
1606 : {
1607 0 : GEN R = zerocol(nbi+1);
1608 0 : for (i=1; i<=lR; i++)
1609 0 : gel(R,pcol[i]) = gel(Rp,i);
1610 0 : return gc_GEN(av, R);
1611 : }
1612 0 : for (i = 1; i <= lR; ++i)
1613 0 : if (signe(gel(Rp, i)))
1614 0 : return gc_leaf(av, mkvecsmall(pcol[i]));
1615 : return NULL; /* LCOV_EXCL_LINE */
1616 : }
1617 :
1618 : GEN
1619 0 : FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1620 : {
1621 0 : return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1622 : }
1623 :
1624 : GEN
1625 0 : FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1626 : {
1627 : GEN res;
1628 0 : pari_CATCH(e_INV)
1629 : {
1630 0 : GEN E = pari_err_last();
1631 0 : GEN x = err_get_compo(E,2);
1632 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1633 0 : if (DEBUGLEVEL)
1634 0 : pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1635 0 : res = NULL;
1636 : } pari_TRY {
1637 0 : res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1638 0 : } pari_ENDCATCH
1639 0 : return res;
1640 : }
1641 :
1642 : static GEN
1643 147 : FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1644 : {
1645 147 : long i, j, oldf = 0, f = 0;
1646 147 : long lrow = lg(prow), lM = lg(M);
1647 147 : GEN W = const_vecsmall(lM-1,1);
1648 147 : GEN R = cgetg(lrow, t_VEC);
1649 228354 : for (i=1; i<lrow; i++)
1650 228207 : gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1651 : do
1652 : {
1653 442 : oldf = f;
1654 374995 : for(i=1; i<lM; i++)
1655 374553 : if (W[i])
1656 : {
1657 131804 : GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1658 131804 : long c=0, cj=0, lF = lg(F);
1659 1094258 : for(j=1; j<lF; j++)
1660 962454 : if (!gel(R,F[j])) { c++; cj=j; }
1661 131804 : if (c>=2) continue;
1662 112308 : if (c==1)
1663 : {
1664 32895 : pari_sp av = avma;
1665 32895 : GEN s = gen_0;
1666 272914 : for(j=1; j<lF; j++)
1667 240019 : if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1668 : /* s /= -E[cj] mod p */
1669 32895 : s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1670 32895 : gel(R,F[cj]) = gc_INT(av, s);
1671 32895 : f++;
1672 : }
1673 112308 : W[i]=0;
1674 : }
1675 442 : } while(oldf!=f);
1676 228354 : for (i=1; i<lrow; i++)
1677 228207 : if (!gel(R,i)) gel(R,i) = gen_0;
1678 147 : return R;
1679 : }
1680 :
1681 : /* Return a linear form Y such that YM is essentially 0 */
1682 : GEN
1683 147 : FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1684 : {
1685 147 : pari_sp av = avma, av2;
1686 : GEN Mp, pcol, prow;
1687 : long i, n;
1688 : pari_timer ti;
1689 : struct wrapper_modp_s W;
1690 147 : if (DEBUGLEVEL) timer_start(&ti);
1691 147 : RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1692 147 : if (!pcol) return gc_NULL(av);
1693 147 : if (DEBUGLEVEL)
1694 0 : timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1695 147 : n = lg(pcol)-1;
1696 147 : Mp = cgetg(n+1, t_MAT);
1697 23186 : for (i=1; i<=n; i++) gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1698 147 : W.E = (void*)Mp;
1699 147 : W.p = p;
1700 147 : av2 = avma;
1701 0 : for(;; set_avma(av2))
1702 0 : {
1703 147 : GEN R, Rp, B = random_FpV(n, p), MB = FpV_FpMs_mul(B, Mp, p);
1704 147 : if (DEBUGLEVEL) timer_start(&ti);
1705 147 : pari_CATCH(e_INV)
1706 : {
1707 0 : GEN E = pari_err_last(), x = err_get_compo(E,2);
1708 0 : if (typ(x) != t_INTMOD) pari_err(0,E);
1709 0 : if (DEBUGLEVEL)
1710 0 : pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1711 0 : Rp = NULL;
1712 : } pari_TRY {
1713 147 : Rp = gen_FpM_Wiedemann((void*)&W, wrap_relker, MB, p);
1714 147 : } pari_ENDCATCH
1715 147 : if (!Rp) continue;
1716 147 : if (typ(Rp)==t_COL)
1717 : {
1718 147 : Rp = FpC_sub(Rp,B,p);
1719 147 : if (ZV_equal0(Rp)) continue;
1720 : }
1721 147 : R = FpMs_structelim_back(M, Rp, prow, p);
1722 147 : if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1723 147 : return gc_GEN(av, R);
1724 : }
1725 : }
1726 :
1727 : GEN
1728 42 : FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1729 : {
1730 42 : return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1731 : }
|
