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