Line data Source code
1 : /* Copyright (C) 2000, 2012 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 : /********************************************************************/
16 : /** **/
17 : /** LINEAR ALGEBRA **/
18 : /** (first part) **/
19 : /** **/
20 : /********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_mat
25 :
26 : /*******************************************************************/
27 : /* */
28 : /* GEREPILE */
29 : /* */
30 : /*******************************************************************/
31 :
32 : static void
33 0 : gerepile_mat(pari_sp av, pari_sp tetpil, GEN x, long k, long m, long n, long t)
34 : {
35 0 : pari_sp A, bot = pari_mainstack->bot;
36 : long u, i;
37 : size_t dec;
38 :
39 0 : (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
40 :
41 0 : for (u=t+1; u<=m; u++)
42 : {
43 0 : A = (pari_sp)coeff(x,u,k);
44 0 : if (A < av && A >= bot) coeff(x,u,k) += dec;
45 : }
46 0 : for (i=k+1; i<=n; i++)
47 0 : for (u=1; u<=m; u++)
48 : {
49 0 : A = (pari_sp)coeff(x,u,i);
50 0 : if (A < av && A >= bot) coeff(x,u,i) += dec;
51 : }
52 0 : }
53 :
54 : static void
55 0 : gen_gerepile_gauss_ker(GEN x, long k, long t, pari_sp av, void *E, GEN (*copy)(void*, GEN))
56 : {
57 0 : pari_sp tetpil = avma;
58 0 : long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
59 :
60 0 : if (DEBUGMEM > 1) pari_warn(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
61 0 : for (u=t+1; u<=m; u++) gcoeff(x,u,k) = copy(E,gcoeff(x,u,k));
62 0 : for (i=k+1; i<=n; i++)
63 0 : for (u=1; u<=m; u++) gcoeff(x,u,i) = copy(E,gcoeff(x,u,i));
64 0 : gerepile_mat(av,tetpil,x,k,m,n,t);
65 0 : }
66 :
67 : /* special gerepile for huge matrices */
68 :
69 : #define COPY(x) {\
70 : GEN _t = (x); if (!is_universal_constant(_t)) x = gcopy(_t); \
71 : }
72 :
73 : INLINE GEN
74 0 : _copy(void *E, GEN x)
75 : {
76 0 : (void) E; COPY(x);
77 0 : return x;
78 : }
79 :
80 : static void
81 0 : gerepile_gauss_ker(GEN x, long k, long t, pari_sp av)
82 : {
83 0 : gen_gerepile_gauss_ker(x, k, t, av, NULL, &_copy);
84 0 : }
85 :
86 : static void
87 0 : gerepile_gauss(GEN x,long k,long t,pari_sp av, long j, GEN c)
88 : {
89 0 : pari_sp tetpil = avma, A, bot;
90 0 : long u,i, n = lg(x)-1, m = n? nbrows(x): 0;
91 : size_t dec;
92 :
93 0 : if (DEBUGMEM > 1) pari_warn(warnmem,"RgM_pivots. k=%ld, n=%ld",k,n);
94 0 : for (u=t+1; u<=m; u++)
95 0 : if (u==j || !c[u]) COPY(gcoeff(x,u,k));
96 0 : for (u=1; u<=m; u++)
97 0 : if (u==j || !c[u])
98 0 : for (i=k+1; i<=n; i++) COPY(gcoeff(x,u,i));
99 :
100 0 : (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
101 0 : bot = pari_mainstack->bot;
102 0 : for (u=t+1; u<=m; u++)
103 0 : if (u==j || !c[u])
104 : {
105 0 : A=(pari_sp)coeff(x,u,k);
106 0 : if (A<av && A>=bot) coeff(x,u,k)+=dec;
107 : }
108 0 : for (u=1; u<=m; u++)
109 0 : if (u==j || !c[u])
110 0 : for (i=k+1; i<=n; i++)
111 : {
112 0 : A=(pari_sp)coeff(x,u,i);
113 0 : if (A<av && A>=bot) coeff(x,u,i)+=dec;
114 : }
115 0 : }
116 :
117 : /*******************************************************************/
118 : /* */
119 : /* GENERIC */
120 : /* */
121 : /*******************************************************************/
122 : GEN
123 1271 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
124 : {
125 1271 : pari_sp av0 = avma, av, tetpil;
126 : GEN y, c, d;
127 : long i, j, k, r, t, n, m;
128 :
129 1271 : n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
130 1271 : m=nbrows(x); r=0;
131 1271 : x = RgM_shallowcopy(x);
132 1271 : c = zero_zv(m);
133 1271 : d=new_chunk(n+1);
134 1271 : av=avma;
135 4558 : for (k=1; k<=n; k++)
136 : {
137 9501 : for (j=1; j<=m; j++)
138 8047 : if (!c[j])
139 : {
140 5583 : gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
141 5583 : if (!ff->equal0(gcoeff(x,j,k))) break;
142 : }
143 3322 : if (j>m)
144 : {
145 1454 : if (deplin)
146 : {
147 35 : GEN c = cgetg(n+1, t_COL), g0 = ff->s(E,0), g1=ff->s(E,1);
148 98 : for (i=1; i<k; i++) gel(c,i) = ff->red(E, gcoeff(x,d[i],k));
149 63 : gel(c,k) = g1; for (i=k+1; i<=n; i++) gel(c,i) = g0;
150 35 : return gerepileupto(av0, c);
151 : }
152 1419 : r++; d[k]=0;
153 3313 : for(j=1; j<k; j++)
154 1894 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
155 : }
156 : else
157 : {
158 1868 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
159 1868 : c[j] = k; d[k] = j;
160 1868 : gcoeff(x,j,k) = ff->s(E,-1);
161 4554 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
162 9916 : for (t=1; t<=m; t++)
163 : {
164 8048 : if (t==j) continue;
165 :
166 6180 : piv = ff->red(E,gcoeff(x,t,k));
167 6180 : if (ff->equal0(piv)) continue;
168 :
169 2240 : gcoeff(x,t,k) = ff->s(E,0);
170 5527 : for (i=k+1; i<=n; i++)
171 3287 : gcoeff(x,t,i) = ff->red(E, ff->add(E, gcoeff(x,t,i),
172 3287 : ff->mul(E,piv,gcoeff(x,j,i))));
173 2240 : if (gc_needed(av,1))
174 0 : gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
175 : }
176 : }
177 : }
178 1236 : if (deplin) return gc_NULL(av0);
179 :
180 1208 : tetpil=avma; y=cgetg(r+1,t_MAT);
181 2627 : for (j=k=1; j<=r; j++,k++)
182 : {
183 1419 : GEN C = cgetg(n+1,t_COL);
184 1419 : GEN g0 = ff->s(E,0), g1 = ff->s(E,1);
185 2640 : gel(y,j) = C; while (d[k]) k++;
186 3313 : for (i=1; i<k; i++)
187 1894 : if (d[i])
188 : {
189 1512 : GEN p1=gcoeff(x,d[i],k);
190 1512 : gel(C,i) = ff->red(E,p1); gunclone(p1);
191 : }
192 : else
193 382 : gel(C,i) = g0;
194 2096 : gel(C,k) = g1; for (i=k+1; i<=n; i++) gel(C,i) = g0;
195 : }
196 1208 : return gerepile(av0,tetpil,y);
197 : }
198 :
199 : GEN
200 1119 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
201 : {
202 : pari_sp av;
203 : GEN c, d;
204 1119 : long i, j, k, r, t, m, n = lg(x)-1;
205 :
206 1119 : if (!n) { *rr = 0; return NULL; }
207 :
208 1119 : m=nbrows(x); r=0;
209 1119 : d = cgetg(n+1, t_VECSMALL);
210 1119 : x = RgM_shallowcopy(x);
211 1119 : c = zero_zv(m);
212 1119 : av=avma;
213 3816 : for (k=1; k<=n; k++)
214 : {
215 7828 : for (j=1; j<=m; j++)
216 7220 : if (!c[j])
217 : {
218 5367 : gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
219 5367 : if (!ff->equal0(gcoeff(x,j,k))) break;
220 : }
221 2697 : if (j>m) { r++; d[k]=0; }
222 : else
223 : {
224 2089 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
225 2089 : GEN g0 = ff->s(E,0);
226 2089 : c[j] = k; d[k] = j;
227 4012 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
228 11858 : for (t=1; t<=m; t++)
229 : {
230 9769 : if (c[t]) continue; /* already a pivot on that line */
231 :
232 6344 : piv = ff->red(E,gcoeff(x,t,k));
233 6344 : if (ff->equal0(piv)) continue;
234 3021 : gcoeff(x,t,k) = g0;
235 4712 : for (i=k+1; i<=n; i++)
236 1691 : gcoeff(x,t,i) = ff->red(E, ff->add(E,gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i))));
237 3021 : if (gc_needed(av,1))
238 0 : gerepile_gauss(x,k,t,av,j,c);
239 : }
240 6101 : for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
241 : }
242 : }
243 1119 : *rr = r; return gc_const((pari_sp)d, d);
244 : }
245 :
246 : GEN
247 294 : gen_det(GEN a, void *E, const struct bb_field *ff)
248 : {
249 294 : pari_sp av = avma;
250 294 : long i,j,k, s = 1, nbco = lg(a)-1;
251 294 : GEN x = ff->s(E,1);
252 294 : if (!nbco) return x;
253 287 : a = RgM_shallowcopy(a);
254 1064 : for (i=1; i<nbco; i++)
255 : {
256 : GEN q;
257 1029 : for(k=i; k<=nbco; k++)
258 : {
259 994 : gcoeff(a,k,i) = ff->red(E,gcoeff(a,k,i));
260 994 : if (!ff->equal0(gcoeff(a,k,i))) break;
261 : }
262 812 : if (k > nbco) return gerepileupto(av, gcoeff(a,i,i));
263 777 : if (k != i)
264 : { /* exchange the lines s.t. k = i */
265 413 : for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
266 105 : s = -s;
267 : }
268 777 : q = gcoeff(a,i,i);
269 777 : x = ff->red(E,ff->mul(E,x,q));
270 777 : q = ff->inv(E,q);
271 2324 : for (k=i+1; k<=nbco; k++)
272 : {
273 1547 : GEN m = ff->red(E,gcoeff(a,i,k));
274 1547 : if (ff->equal0(m)) continue;
275 1092 : m = ff->neg(E, ff->red(E,ff->mul(E,m, q)));
276 3528 : for (j=i+1; j<=nbco; j++)
277 2436 : gcoeff(a,j,k) = ff->red(E, ff->add(E, gcoeff(a,j,k),
278 2436 : ff->mul(E, m, gcoeff(a,j,i))));
279 : }
280 777 : if (gc_needed(av,2))
281 : {
282 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
283 0 : gerepileall(av,2, &a,&x);
284 : }
285 : }
286 252 : if (s < 0) x = ff->neg(E,x);
287 252 : return gerepileupto(av, ff->red(E,ff->mul(E, x, gcoeff(a,nbco,nbco))));
288 : }
289 :
290 : INLINE void
291 57335 : _gen_addmul(GEN b, long k, long i, GEN m, void *E, const struct bb_field *ff)
292 : {
293 57335 : gel(b,i) = ff->red(E,gel(b,i));
294 57335 : gel(b,k) = ff->add(E,gel(b,k), ff->mul(E,m, gel(b,i)));
295 57335 : }
296 :
297 : static GEN
298 21485 : _gen_get_col(GEN a, GEN b, long li, void *E, const struct bb_field *ff)
299 : {
300 21485 : GEN u = cgetg(li+1,t_COL);
301 21485 : pari_sp av = avma;
302 : long i, j;
303 :
304 21485 : gel(u,li) = gerepileupto(av, ff->red(E,ff->mul(E,gel(b,li), gcoeff(a,li,li))));
305 94317 : for (i=li-1; i>0; i--)
306 : {
307 72832 : pari_sp av = avma;
308 72832 : GEN m = gel(b,i);
309 296314 : for (j=i+1; j<=li; j++) m = ff->add(E,m, ff->neg(E,ff->mul(E,gcoeff(a,i,j), gel(u,j))));
310 72832 : m = ff->red(E, m);
311 72832 : gel(u,i) = gerepileupto(av, ff->red(E,ff->mul(E,m, gcoeff(a,i,i))));
312 : }
313 21485 : return u;
314 : }
315 :
316 : GEN
317 5816 : gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff)
318 : {
319 : long i, j, k, li, bco, aco;
320 5816 : GEN u, g0 = ff->s(E,0);
321 5816 : pari_sp av = avma;
322 5816 : a = RgM_shallowcopy(a);
323 5816 : b = RgM_shallowcopy(b);
324 5816 : aco = lg(a)-1; bco = lg(b)-1; li = nbrows(a);
325 20204 : for (i=1; i<=aco; i++)
326 : {
327 : GEN invpiv;
328 22393 : for (k = i; k <= li; k++)
329 : {
330 22351 : GEN piv = ff->red(E,gcoeff(a,k,i));
331 22351 : if (!ff->equal0(piv)) { gcoeff(a,k,i) = ff->inv(E,piv); break; }
332 2189 : gcoeff(a,k,i) = g0;
333 : }
334 : /* found a pivot on line k */
335 20204 : if (k > li) return NULL;
336 20162 : if (k != i)
337 : { /* swap lines so that k = i */
338 9276 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
339 12419 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
340 : }
341 20162 : if (i == aco) break;
342 :
343 14388 : invpiv = gcoeff(a,i,i); /* 1/piv mod p */
344 51602 : for (k=i+1; k<=li; k++)
345 : {
346 37214 : GEN m = ff->red(E,gcoeff(a,k,i)); gcoeff(a,k,i) = g0;
347 37214 : if (ff->equal0(m)) continue;
348 :
349 7799 : m = ff->red(E,ff->neg(E,ff->mul(E,m, invpiv)));
350 29376 : for (j=i+1; j<=aco; j++) _gen_addmul(gel(a,j),k,i,m,E,ff);
351 43557 : for (j=1 ; j<=bco; j++) _gen_addmul(gel(b,j),k,i,m,E,ff);
352 : }
353 14388 : if (gc_needed(av,1))
354 : {
355 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_Gauss. i=%ld",i);
356 0 : gerepileall(av,2, &a,&b);
357 : }
358 : }
359 :
360 5774 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
361 5774 : u = cgetg(bco+1,t_MAT);
362 27259 : for (j=1; j<=bco; j++) gel(u,j) = _gen_get_col(a, gel(b,j), aco, E, ff);
363 5774 : return u;
364 : }
365 :
366 : /* compatible t_MAT * t_COL, lgA = lg(A) = lg(B) > 1, l = lgcols(A) */
367 : static GEN
368 350531 : gen_matcolmul_i(GEN A, GEN B, ulong lgA, ulong l,
369 : void *E, const struct bb_field *ff)
370 : {
371 350531 : GEN C = cgetg(l, t_COL);
372 : ulong i;
373 2231912 : for (i = 1; i < l; i++) {
374 1881381 : pari_sp av = avma;
375 1881381 : GEN e = ff->mul(E, gcoeff(A, i, 1), gel(B, 1));
376 : ulong k;
377 5531603 : for(k = 2; k < lgA; k++)
378 3650222 : e = ff->add(E, e, ff->mul(E, gcoeff(A, i, k), gel(B, k)));
379 1881381 : gel(C, i) = gerepileupto(av, ff->red(E, e));
380 : }
381 350531 : return C;
382 : }
383 :
384 : GEN
385 48648 : gen_matcolmul(GEN A, GEN B, void *E, const struct bb_field *ff)
386 : {
387 48648 : ulong lgA = lg(A);
388 48648 : if (lgA != (ulong)lg(B))
389 0 : pari_err_OP("operation 'gen_matcolmul'", A, B);
390 48648 : if (lgA == 1)
391 0 : return cgetg(1, t_COL);
392 48648 : return gen_matcolmul_i(A, B, lgA, lgcols(A), E, ff);
393 : }
394 :
395 : static GEN
396 75983 : gen_matmul_classical(GEN A, GEN B, long l, long la, long lb,
397 : void *E, const struct bb_field *ff)
398 : {
399 : long j;
400 75983 : GEN C = cgetg(lb, t_MAT);
401 377866 : for(j = 1; j < lb; j++)
402 301883 : gel(C, j) = gen_matcolmul_i(A, gel(B, j), la, l, E, ff);
403 75983 : return C;
404 : }
405 :
406 : /* Strassen-Winograd algorithm */
407 :
408 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
409 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
410 : static GEN
411 0 : add_slices(long m, long n,
412 : GEN A, long ma, long da, long na, long ea,
413 : GEN B, long mb, long db, long nb, long eb,
414 : void *E, const struct bb_field *ff)
415 : {
416 0 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
417 0 : GEN M = cgetg(n + 1, t_MAT), C;
418 :
419 0 : for (j = 1; j <= min_e; j++) {
420 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
421 0 : for (i = 1; i <= min_d; i++)
422 0 : gel(C, i) = ff->add(E, gcoeff(A, ma + i, na + j),
423 0 : gcoeff(B, mb + i, nb + j));
424 0 : for (; i <= da; i++)
425 0 : gel(C, i) = gcoeff(A, ma + i, na + j);
426 0 : for (; i <= db; i++)
427 0 : gel(C, i) = gcoeff(B, mb + i, nb + j);
428 0 : for (; i <= m; i++)
429 0 : gel(C, i) = ff->s(E, 0);
430 : }
431 0 : for (; j <= ea; j++) {
432 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
433 0 : for (i = 1; i <= da; i++)
434 0 : gel(C, i) = gcoeff(A, ma + i, na + j);
435 0 : for (; i <= m; i++)
436 0 : gel(C, i) = ff->s(E, 0);
437 : }
438 0 : for (; j <= eb; j++) {
439 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
440 0 : for (i = 1; i <= db; i++)
441 0 : gel(C, i) = gcoeff(B, mb + i, nb + j);
442 0 : for (; i <= m; i++)
443 0 : gel(C, i) = ff->s(E, 0);
444 : }
445 0 : for (; j <= n; j++) {
446 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
447 0 : for (i = 1; i <= m; i++)
448 0 : gel(C, i) = ff->s(E, 0);
449 : }
450 0 : return M;
451 : }
452 :
453 : /* Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
454 : * as an (m x n)-matrix, padding the input with zeroes as necessary. */
455 : static GEN
456 0 : subtract_slices(long m, long n,
457 : GEN A, long ma, long da, long na, long ea,
458 : GEN B, long mb, long db, long nb, long eb,
459 : void *E, const struct bb_field *ff)
460 : {
461 0 : long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
462 0 : GEN M = cgetg(n + 1, t_MAT), C;
463 :
464 0 : for (j = 1; j <= min_e; j++) {
465 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
466 0 : for (i = 1; i <= min_d; i++)
467 0 : gel(C, i) = ff->add(E, gcoeff(A, ma + i, na + j),
468 0 : ff->neg(E, gcoeff(B, mb + i, nb + j)));
469 0 : for (; i <= da; i++)
470 0 : gel(C, i) = gcoeff(A, ma + i, na + j);
471 0 : for (; i <= db; i++)
472 0 : gel(C, i) = ff->neg(E, gcoeff(B, mb + i, nb + j));
473 0 : for (; i <= m; i++)
474 0 : gel(C, i) = ff->s(E, 0);
475 : }
476 0 : for (; j <= ea; j++) {
477 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
478 0 : for (i = 1; i <= da; i++)
479 0 : gel(C, i) = gcoeff(A, ma + i, na + j);
480 0 : for (; i <= m; i++)
481 0 : gel(C, i) = ff->s(E, 0);
482 : }
483 0 : for (; j <= eb; j++) {
484 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
485 0 : for (i = 1; i <= db; i++)
486 0 : gel(C, i) = ff->neg(E, gcoeff(B, mb + i, nb + j));
487 0 : for (; i <= m; i++)
488 0 : gel(C, i) = ff->s(E, 0);
489 : }
490 0 : for (; j <= n; j++) {
491 0 : gel(M, j) = C = cgetg(m + 1, t_COL);
492 0 : for (i = 1; i <= m; i++)
493 0 : gel(C, i) = ff->s(E, 0);
494 : }
495 0 : return M;
496 : }
497 :
498 : static GEN gen_matmul_i(GEN A, GEN B, long l, long la, long lb,
499 : void *E, const struct bb_field *ff);
500 :
501 : static GEN
502 0 : gen_matmul_sw(GEN A, GEN B, long m, long n, long p,
503 : void *E, const struct bb_field *ff)
504 : {
505 0 : pari_sp av = avma;
506 0 : long m1 = (m + 1)/2, m2 = m/2,
507 0 : n1 = (n + 1)/2, n2 = n/2,
508 0 : p1 = (p + 1)/2, p2 = p/2;
509 : GEN A11, A12, A22, B11, B21, B22,
510 : S1, S2, S3, S4, T1, T2, T3, T4,
511 : M1, M2, M3, M4, M5, M6, M7,
512 : V1, V2, V3, C11, C12, C21, C22, C;
513 :
514 0 : T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, E, ff);
515 0 : S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, E, ff);
516 0 : M2 = gen_matmul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, E, ff);
517 0 : if (gc_needed(av, 1))
518 0 : gerepileall(av, 2, &T2, &M2); /* destroy S1 */
519 0 : T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, E, ff);
520 0 : if (gc_needed(av, 1))
521 0 : gerepileall(av, 2, &M2, &T3); /* destroy T2 */
522 0 : S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, E, ff);
523 0 : T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, E, ff);
524 0 : M3 = gen_matmul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, E, ff);
525 0 : if (gc_needed(av, 1))
526 0 : gerepileall(av, 4, &M2, &T3, &S2, &M3); /* destroy T1 */
527 0 : S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, E, ff);
528 0 : if (gc_needed(av, 1))
529 0 : gerepileall(av, 4, &M2, &T3, &M3, &S3); /* destroy S2 */
530 0 : A11 = matslice(A, 1, m1, 1, n1);
531 0 : B11 = matslice(B, 1, n1, 1, p1);
532 0 : M1 = gen_matmul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, E, ff);
533 0 : if (gc_needed(av, 1))
534 0 : gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1); /* destroy A11, B11 */
535 0 : A12 = matslice(A, 1, m1, n1 + 1, n);
536 0 : B21 = matslice(B, n1 + 1, n, 1, p1);
537 0 : M4 = gen_matmul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, E, ff);
538 0 : if (gc_needed(av, 1))
539 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4); /* destroy A12, B21 */
540 0 : C11 = add_slices(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, E, ff);
541 0 : if (gc_needed(av, 1))
542 0 : gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &C11); /* destroy M4 */
543 0 : M5 = gen_matmul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, E, ff);
544 0 : S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, E, ff);
545 0 : if (gc_needed(av, 1))
546 0 : gerepileall(av, 7, &M2, &T3, &M3, &M1, &C11, &M5, &S4); /* destroy S3 */
547 0 : T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, E, ff);
548 0 : if (gc_needed(av, 1))
549 0 : gerepileall(av, 7, &M2, &M3, &M1, &C11, &M5, &S4, &T4); /* destroy T3 */
550 0 : V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, E, ff);
551 0 : if (gc_needed(av, 1))
552 0 : gerepileall(av, 6, &M2, &M3, &S4, &T4, &C11, &V1); /* destroy M1, M5 */
553 0 : B22 = matslice(B, n1 + 1, n, p1 + 1, p);
554 0 : M6 = gen_matmul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, E, ff);
555 0 : if (gc_needed(av, 1))
556 0 : gerepileall(av, 6, &M2, &M3, &T4, &C11, &V1, &M6); /* destroy S4, B22 */
557 0 : A22 = matslice(A, m1 + 1, m, n1 + 1, n);
558 0 : M7 = gen_matmul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, E, ff);
559 0 : if (gc_needed(av, 1))
560 0 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M6, &M7); /* destroy A22, T4 */
561 0 : V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, E, ff);
562 0 : C12 = add_slices(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, E, ff);
563 0 : if (gc_needed(av, 1))
564 0 : gerepileall(av, 6, &M2, &M3, &C11, &V1, &M7, &C12); /* destroy V3, M6 */
565 0 : V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, E, ff);
566 0 : if (gc_needed(av, 1))
567 0 : gerepileall(av, 5, &M3, &C11, &M7, &C12, &V2); /* destroy V1, M2 */
568 0 : C21 = add_slices(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, E, ff);
569 0 : if (gc_needed(av, 1))
570 0 : gerepileall(av, 5, &M3, &C11, &C12, &V2, &C21); /* destroy M7 */
571 0 : C22 = add_slices(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, E, ff);
572 0 : if (gc_needed(av, 1))
573 0 : gerepileall(av, 4, &C11, &C12, &C21, &C22); /* destroy V2, M3 */
574 0 : C = mkmat2(mkcol2(C11, C21), mkcol2(C12, C22));
575 0 : return gerepileupto(av, matconcat(C));
576 : }
577 :
578 : /* Strassen-Winograd used for dim >= gen_matmul_sw_bound */
579 : static const long gen_matmul_sw_bound = 24;
580 :
581 : static GEN
582 75983 : gen_matmul_i(GEN A, GEN B, long l, long la, long lb,
583 : void *E, const struct bb_field *ff)
584 : {
585 75983 : if (l <= gen_matmul_sw_bound
586 7 : || la <= gen_matmul_sw_bound
587 0 : || lb <= gen_matmul_sw_bound)
588 75983 : return gen_matmul_classical(A, B, l, la, lb, E, ff);
589 : else
590 0 : return gen_matmul_sw(A, B, l - 1, la - 1, lb - 1, E, ff);
591 : }
592 :
593 : GEN
594 75983 : gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
595 : {
596 75983 : ulong lgA, lgB = lg(B);
597 75983 : if (lgB == 1)
598 0 : return cgetg(1, t_MAT);
599 75983 : lgA = lg(A);
600 75983 : if (lgA != (ulong)lgcols(B))
601 0 : pari_err_OP("operation 'gen_matmul'", A, B);
602 75983 : if (lgA == 1)
603 0 : return zeromat(0, lgB - 1);
604 75983 : return gen_matmul_i(A, B, lgcols(A), lgA, lgB, E, ff);
605 : }
606 :
607 : static GEN
608 14704 : gen_colneg(GEN x, void *E, const struct bb_field *ff)
609 59481 : { pari_APPLY_same(ff->neg(E, gel(x,i))); }
610 :
611 : static GEN
612 2450 : gen_matneg(GEN x, void *E, const struct bb_field *ff)
613 17084 : { pari_APPLY_same(gen_colneg(gel(x,i), E, ff)); }
614 :
615 : static GEN
616 310391 : gen_colscalmul(GEN x, GEN b, void *E, const struct bb_field *ff)
617 734209 : { pari_APPLY_same(ff->red(E, ff->mul(E, gel(x,i), b))); }
618 :
619 : static GEN
620 47140 : gen_matscalmul(GEN x, GEN b, void *E, const struct bb_field *ff)
621 357531 : { pari_APPLY_same(gen_colscalmul(gel(x,i), b, E, ff)); }
622 :
623 : static GEN
624 574403 : gen_colsub(GEN x, GEN y, void *E, const struct bb_field *ff)
625 2141596 : { pari_APPLY_same(ff->add(E, gel(x,i), ff->neg(E, gel(y,i)))); }
626 :
627 : static GEN
628 63502 : gen_matsub(GEN x, GEN y, void *E, const struct bb_field *ff)
629 637905 : { pari_APPLY_same(gen_colsub(gel(x,i), gel(y,i), E, ff)); }
630 :
631 : static GEN
632 13108 : gen_zerocol(long n, void* data, const struct bb_field *R)
633 13108 : { return const_col(n, R->s(data, 0)); }
634 :
635 : static GEN
636 13108 : gen_zeromat(long m, long n, void* data, const struct bb_field *R)
637 : {
638 13108 : GEN M = const_vec(n, gen_zerocol(m, data, R));
639 13108 : settyp(M, t_MAT); return M;
640 : }
641 :
642 : static GEN
643 154 : gen_colei(long n, long i, void *E, const struct bb_field *S)
644 : {
645 154 : GEN y = cgetg(n+1,t_COL), _0, _1;
646 : long j;
647 154 : if (n < 0) pari_err_DOMAIN("gen_colei", "dimension","<",gen_0,stoi(n));
648 154 : _0 = S->s(E,0);
649 154 : _1 = S->s(E,1);
650 2422 : for (j=1; j<=n; j++)
651 2268 : gel(y, j) = i==j ? _1: _0;
652 154 : return y;
653 : }
654 :
655 : /* assume dim A >= 1, A invertible + upper triangular */
656 : static GEN
657 77 : gen_matinv_upper_ind(GEN A, long index, void *E, const struct bb_field *ff)
658 : {
659 77 : long n = lg(A) - 1, i, j;
660 77 : GEN u = cgetg(n + 1, t_COL);
661 147 : for (i = n; i > index; i--)
662 70 : gel(u, i) = ff->s(E, 0);
663 77 : gel(u, i) = ff->inv(E, gcoeff(A, i, i));
664 147 : for (i--; i > 0; i--) {
665 70 : pari_sp av = avma;
666 70 : GEN m = ff->neg(E, ff->mul(E, gcoeff(A, i, i + 1), gel(u, i + 1)));
667 112 : for (j = i + 2; j <= n; j++)
668 42 : m = ff->add(E, m, ff->neg(E, ff->mul(E, gcoeff(A, i, j), gel(u, j))));
669 70 : gel(u, i) = gerepileupto(av, ff->red(E, ff->mul(E, m, ff->inv(E, gcoeff(A, i, i)))));
670 : }
671 77 : return u;
672 : }
673 :
674 : static GEN
675 28 : gen_matinv_upper(GEN A, void *E, const struct bb_field *ff)
676 : {
677 : long i, l;
678 28 : GEN B = cgetg_copy(A, &l);
679 105 : for (i = 1; i < l; i++)
680 77 : gel(B,i) = gen_matinv_upper_ind(A, i, E, ff);
681 28 : return B;
682 : }
683 :
684 : /* find z such that A z = y. Return NULL if no solution */
685 : GEN
686 0 : gen_matcolinvimage(GEN A, GEN y, void *E, const struct bb_field *ff)
687 : {
688 0 : pari_sp av = avma;
689 0 : long i, l = lg(A);
690 : GEN M, x, t;
691 :
692 0 : M = gen_ker(shallowconcat(A, y), 0, E, ff);
693 0 : i = lg(M) - 1;
694 0 : if (!i) return gc_NULL(av);
695 :
696 0 : x = gel(M, i);
697 0 : t = gel(x, l);
698 0 : if (ff->equal0(t)) return gc_NULL(av);
699 :
700 0 : t = ff->neg(E, ff->inv(E, t));
701 0 : setlg(x, l);
702 0 : for (i = 1; i < l; i++)
703 0 : gel(x, i) = ff->red(E, ff->mul(E, t, gel(x, i)));
704 0 : return gc_GEN(av, x);
705 : }
706 :
707 : /* find Z such that A Z = B. Return NULL if no solution */
708 : GEN
709 77 : gen_matinvimage(GEN A, GEN B, void *E, const struct bb_field *ff)
710 : {
711 77 : pari_sp av = avma;
712 : GEN d, x, X, Y;
713 : long i, j, nY, nA, nB;
714 77 : x = gen_ker(shallowconcat(gen_matneg(A, E, ff), B), 0, E, ff);
715 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
716 : * We must find T such that Y T = Id_nB then X T = Z. This exists
717 : * iff Y has at least nB columns and full rank. */
718 77 : nY = lg(x) - 1;
719 77 : nB = lg(B) - 1;
720 77 : if (nY < nB) return gc_NULL(av);
721 77 : nA = lg(A) - 1;
722 77 : Y = rowslice(x, nA + 1, nA + nB); /* nB rows */
723 77 : d = cgetg(nB + 1, t_VECSMALL);
724 182 : for (i = nB, j = nY; i >= 1; i--, j--) {
725 224 : for (; j >= 1; j--)
726 175 : if (!ff->equal0(gcoeff(Y, i, j))) { d[i] = j; break; }
727 154 : if (!j) return gc_NULL(av);
728 : }
729 : /* reduce to the case Y square, upper triangular with 1s on diagonal */
730 28 : Y = vecpermute(Y, d);
731 28 : x = vecpermute(x, d);
732 28 : X = rowslice(x, 1, nA);
733 28 : return gerepileupto(av, gen_matmul(X, gen_matinv_upper(Y, E, ff), E, ff));
734 : }
735 :
736 : static GEN
737 371491 : image_from_pivot(GEN x, GEN d, long r)
738 : {
739 : GEN y;
740 : long j, k;
741 :
742 371491 : if (!d) return gcopy(x);
743 : /* d left on stack for efficiency */
744 365764 : r = lg(x)-1 - r; /* = dim Im(x) */
745 365764 : y = cgetg(r+1,t_MAT);
746 2098572 : for (j=k=1; j<=r; k++)
747 1732808 : if (d[k]) gel(y,j++) = gcopy(gel(x,k));
748 365764 : return y;
749 : }
750 :
751 : /* r = dim Ker x, n = nbrows(x) */
752 : static GEN
753 267807 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
754 : {
755 : pari_sp av;
756 : GEN y, c;
757 267807 : long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
758 :
759 267807 : if (rx == n && r == 0) return gcopy(x);
760 197393 : y = cgetg(n+1, t_MAT);
761 197394 : av = avma; c = zero_zv(n);
762 : /* c = lines containing pivots (could get it from RgM_pivots, but cheap)
763 : * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
764 836160 : for (k = j = 1; j<=rx; j++)
765 638766 : if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
766 1198222 : for (j=1; j<=n; j++)
767 1000828 : if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
768 197394 : set_avma(av);
769 :
770 197394 : rx -= r;
771 836088 : for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
772 559522 : for ( ; j<=n; j++) gel(y,j) = ei(n, y[j]);
773 197393 : return y;
774 : }
775 :
776 : /* n = dim x, r = dim Ker(x), d from RgM_pivots */
777 : static GEN
778 1940948 : indexrank0(long n, long r, GEN d)
779 : {
780 1940948 : GEN p1, p2, res = cgetg(3,t_VEC);
781 : long i, j;
782 :
783 1940944 : r = n - r; /* now r = dim Im(x) */
784 1940944 : p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
785 1940944 : p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
786 1940945 : if (d)
787 : {
788 7831893 : for (i=0,j=1; j<=n; j++)
789 5894455 : if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
790 1937438 : vecsmall_sort(p1);
791 : }
792 1940955 : return res;
793 : }
794 :
795 : /*******************************************************************/
796 : /* */
797 : /* Echelon form and CUP decomposition */
798 : /* */
799 : /*******************************************************************/
800 :
801 : /* By Peter Bruin, based on
802 : C.-P. Jeannerod, C. Pernet and A. Storjohann, Rank-profile revealing
803 : Gaussian elimination and the CUP matrix decomposition. J. Symbolic
804 : Comput. 56 (2013), 46-68.
805 :
806 : Decompose an m x n-matrix A of rank r as C*U*P, with
807 : - C: m x r-matrix in column echelon form (not necessarily reduced)
808 : with all pivots equal to 1
809 : - U: upper-triangular r x n-matrix
810 : - P: permutation matrix
811 : The pivots of C and the known zeroes in C and U are not necessarily
812 : filled in; instead, we also return the vector R of pivot rows.
813 : Instead of the matrix P, we return the permutation p of [1..n]
814 : (t_VECSMALL) such that P[i,j] = 1 if and only if j = p[i].
815 : */
816 :
817 : /* complement of a strictly increasing subsequence of (1, 2, ..., n) */
818 : static GEN
819 12199 : indexcompl(GEN v, long n)
820 : {
821 12199 : long i, j, k, m = lg(v) - 1;
822 12199 : GEN w = cgetg(n - m + 1, t_VECSMALL);
823 127236 : for (i = j = k = 1; i <= n; i++)
824 115037 : if (j <= m && v[j] == i) j++; else w[k++] = i;
825 12199 : return w;
826 : }
827 :
828 : static GEN
829 4035 : gen_solve_upper_1(GEN U, GEN B, void *E, const struct bb_field *ff)
830 4035 : { return gen_matscalmul(B, ff->inv(E, gcoeff(U, 1, 1)), E, ff); }
831 :
832 : static GEN
833 2256 : gen_rsolve_upper_2(GEN U, GEN B, void *E, const struct bb_field *ff)
834 : {
835 2256 : GEN a = gcoeff(U, 1, 1), b = gcoeff(U, 1, 2), d = gcoeff(U, 2, 2);
836 2256 : GEN D = ff->red(E, ff->mul(E, a, d)), Dinv = ff->inv(E, D);
837 2256 : GEN ainv = ff->red(E, ff->mul(E, d, Dinv));
838 2256 : GEN dinv = ff->red(E, ff->mul(E, a, Dinv));
839 2256 : GEN B1 = rowslice(B, 1, 1);
840 2256 : GEN B2 = rowslice(B, 2, 2);
841 2256 : GEN X2 = gen_matscalmul(B2, dinv, E, ff);
842 2256 : GEN X1 = gen_matscalmul(gen_matsub(B1, gen_matscalmul(X2, b, E, ff), E, ff),
843 : ainv, E, ff);
844 2256 : return vconcat(X1, X2);
845 : }
846 :
847 : /* solve U*X = B, U upper triangular and invertible */
848 : static GEN
849 5840 : gen_rsolve_upper(GEN U, GEN B, void *E, const struct bb_field *ff,
850 : GEN (*mul)(void *E, GEN a, GEN))
851 : {
852 5840 : long n = lg(U) - 1, n1;
853 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
854 5840 : pari_sp av = avma;
855 :
856 5840 : if (n == 0) return B;
857 5840 : if (n == 1) return gen_solve_upper_1(U, B, E, ff);
858 4914 : if (n == 2) return gen_rsolve_upper_2(U, B, E, ff);
859 2658 : n1 = (n + 1)/2;
860 2658 : U2 = vecslice(U, n1 + 1, n);
861 2658 : U11 = matslice(U, 1,n1, 1,n1);
862 2658 : U12 = rowslice(U2, 1, n1);
863 2658 : U22 = rowslice(U2, n1 + 1, n);
864 2658 : B1 = rowslice(B, 1, n1);
865 2658 : B2 = rowslice(B, n1 + 1, n);
866 2658 : X2 = gen_rsolve_upper(U22, B2, E, ff, mul);
867 2658 : B1 = gen_matsub(B1, mul(E, U12, X2), E, ff);
868 2658 : if (gc_needed(av, 1)) gerepileall(av, 3, &B1, &U11, &X2);
869 2658 : X1 = gen_rsolve_upper(U11, B1, E, ff, mul);
870 2658 : X = vconcat(X1, X2);
871 2658 : if (gc_needed(av, 1)) X = gc_GEN(av, X);
872 2658 : return X;
873 : }
874 :
875 : static GEN
876 5894 : gen_lsolve_upper_2(GEN U, GEN B, void *E, const struct bb_field *ff)
877 : {
878 5894 : GEN a = gcoeff(U, 1, 1), b = gcoeff(U, 1, 2), d = gcoeff(U, 2, 2);
879 5894 : GEN D = ff->red(E, ff->mul(E, a, d)), Dinv = ff->inv(E, D);
880 5894 : GEN ainv = ff->red(E, ff->mul(E, d, Dinv)), dinv = ff->red(E, ff->mul(E, a, Dinv));
881 5894 : GEN B1 = vecslice(B, 1, 1);
882 5894 : GEN B2 = vecslice(B, 2, 2);
883 5894 : GEN X1 = gen_matscalmul(B1, ainv, E, ff);
884 5894 : GEN X2 = gen_matscalmul(gen_matsub(B2, gen_matscalmul(X1, b, E, ff), E, ff), dinv, E, ff);
885 5894 : return shallowconcat(X1, X2);
886 : }
887 :
888 : /* solve X*U = B, U upper triangular and invertible */
889 : static GEN
890 13882 : gen_lsolve_upper(GEN U, GEN B, void *E, const struct bb_field *ff,
891 : GEN (*mul)(void *E, GEN a, GEN))
892 : {
893 13882 : long n = lg(U) - 1, n1;
894 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
895 13882 : pari_sp av = avma;
896 :
897 13882 : if (n == 0) return B;
898 13882 : if (n == 1) return gen_solve_upper_1(U, B, E, ff);
899 10773 : if (n == 2) return gen_lsolve_upper_2(U, B, E, ff);
900 4879 : n1 = (n + 1)/2;
901 4879 : U2 = vecslice(U, n1 + 1, n);
902 4879 : U11 = matslice(U, 1,n1, 1,n1);
903 4879 : U12 = rowslice(U2, 1, n1);
904 4879 : U22 = rowslice(U2, n1 + 1, n);
905 4879 : B1 = vecslice(B, 1, n1);
906 4879 : B2 = vecslice(B, n1 + 1, n);
907 4879 : X1 = gen_lsolve_upper(U11, B1, E, ff, mul);
908 4879 : B2 = gen_matsub(B2, mul(E, X1, U12), E, ff);
909 4879 : if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
910 4879 : X2 = gen_lsolve_upper(U22, B2, E, ff, mul);
911 4879 : X = shallowconcat(X1, X2);
912 4879 : if (gc_needed(av, 1)) X = gc_GEN(av, X);
913 4879 : return X;
914 : }
915 :
916 : static GEN
917 12590 : gen_rsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
918 : {
919 12590 : GEN X1 = rowslice(A, 1, 1);
920 12590 : GEN X2 = gen_matsub(rowslice(A, 2, 2), gen_matscalmul(X1, gcoeff(L, 2, 1), E, ff), E, ff);
921 12590 : return vconcat(X1, X2);
922 : }
923 :
924 : /* solve L*X = A, L lower triangular with ones on the diagonal
925 : * (at least as many rows as columns) */
926 : static GEN
927 30422 : gen_rsolve_lower_unit(GEN L, GEN A, void *E, const struct bb_field *ff,
928 : GEN (*mul)(void *E, GEN a, GEN))
929 : {
930 30422 : long m = lg(L) - 1, m1, n;
931 : GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
932 30422 : pari_sp av = avma;
933 :
934 30422 : if (m == 0) return zeromat(0, lg(A) - 1);
935 30422 : if (m == 1) return rowslice(A, 1, 1);
936 24201 : if (m == 2) return gen_rsolve_lower_unit_2(L, A, E, ff);
937 11611 : m1 = (m + 1)/2;
938 11611 : n = nbrows(L);
939 11611 : L1 = vecslice(L, 1, m1);
940 11611 : L11 = rowslice(L1, 1, m1);
941 11611 : L21 = rowslice(L1, m1 + 1, n);
942 11611 : A1 = rowslice(A, 1, m1);
943 11611 : X1 = gen_rsolve_lower_unit(L11, A1, E, ff, mul);
944 11611 : A2 = rowslice(A, m1 + 1, n);
945 11611 : A2 = gen_matsub(A2, mul(E, L21, X1), E, ff);
946 11611 : if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
947 11611 : L22 = matslice(L, m1+1,n, m1+1,m);
948 11611 : X2 = gen_rsolve_lower_unit(L22, A2, E, ff, mul);
949 11611 : X = vconcat(X1, X2);
950 11611 : if (gc_needed(av, 1)) X = gc_GEN(av, X);
951 11611 : return X;
952 : }
953 :
954 : static GEN
955 6065 : gen_lsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
956 : {
957 6065 : GEN X2 = vecslice(A, 2, 2);
958 6065 : GEN X1 = gen_matsub(vecslice(A, 1, 1),
959 6065 : gen_matscalmul(X2, gcoeff(L, 2, 1), E, ff), E, ff);
960 6065 : return shallowconcat(X1, X2);
961 : }
962 :
963 : /* solve L*X = A, L lower triangular with ones on the diagonal
964 : * (at least as many rows as columns) */
965 : static GEN
966 16025 : gen_lsolve_lower_unit(GEN L, GEN A, void *E, const struct bb_field *ff,
967 : GEN (*mul)(void *E, GEN a, GEN))
968 : {
969 16025 : long m = lg(L) - 1, m1;
970 : GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
971 16025 : pari_sp av = avma;
972 :
973 16025 : if (m <= 1) return A;
974 12856 : if (m == 2) return gen_lsolve_lower_unit_2(L, A, E, ff);
975 6791 : m1 = (m + 1)/2;
976 6791 : L2 = vecslice(L, m1 + 1, m);
977 6791 : L22 = rowslice(L2, m1 + 1, m);
978 6791 : A2 = vecslice(A, m1 + 1, m);
979 6791 : X2 = gen_lsolve_lower_unit(L22, A2, E, ff, mul);
980 6791 : if (gc_needed(av, 1)) X2 = gc_GEN(av, X2);
981 6791 : L1 = vecslice(L, 1, m1);
982 6791 : L21 = rowslice(L1, m1 + 1, m);
983 6791 : A1 = vecslice(A, 1, m1);
984 6791 : A1 = gen_matsub(A1, mul(E, X2, L21), E, ff);
985 6791 : L11 = rowslice(L1, 1, m1);
986 6791 : if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
987 6791 : X1 = gen_lsolve_lower_unit(L11, A1, E, ff, mul);
988 6791 : X = shallowconcat(X1, X2);
989 6791 : if (gc_needed(av, 1)) X = gc_GEN(av, X);
990 6791 : return X;
991 : }
992 :
993 : /* destroy A */
994 : static long
995 16007 : gen_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, void *E, const struct bb_field *ff)
996 : {
997 16007 : long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
998 : pari_sp av;
999 : GEN u, v;
1000 :
1001 16007 : if (P) *P = identity_perm(n);
1002 16007 : *R = cgetg(m + 1, t_VECSMALL);
1003 16007 : av = avma;
1004 45918 : for (j = 1, pr = 0; j <= n; j++)
1005 : {
1006 104382 : for (pr++, pc = 0; pr <= m; pr++)
1007 : {
1008 544100 : for (k = j; k <= n; k++)
1009 : {
1010 451722 : v = ff->red(E, gcoeff(A, pr, k));
1011 451722 : gcoeff(A, pr, k) = v;
1012 451722 : if (!pc && !ff->equal0(v)) pc = k;
1013 : }
1014 92378 : if (pc) break;
1015 : }
1016 41915 : if (!pc) break;
1017 29911 : (*R)[j] = pr;
1018 29911 : if (pc != j)
1019 : {
1020 4258 : swap(gel(A, j), gel(A, pc));
1021 4258 : if (P) lswap((*P)[j], (*P)[pc]);
1022 : }
1023 29911 : u = ff->inv(E, gcoeff(A, pr, j));
1024 154973 : for (i = pr + 1; i <= m; i++)
1025 : {
1026 125062 : v = ff->red(E, ff->mul(E, gcoeff(A, i, j), u));
1027 125062 : gcoeff(A, i, j) = v;
1028 125062 : v = ff->neg(E, v);
1029 413234 : for (k = j + 1; k <= n; k++)
1030 288172 : gcoeff(A, i, k) = ff->add(E, gcoeff(A, i, k),
1031 288172 : ff->red(E, ff->mul(E, gcoeff(A, pr, k), v)));
1032 : }
1033 29911 : if (gc_needed(av, 2)) A = gc_GEN(av, A);
1034 : }
1035 16007 : setlg(*R, j);
1036 16007 : *C = vecslice(A, 1, j - 1);
1037 16007 : if (U) *U = rowpermute(A, *R);
1038 16007 : return j - 1;
1039 : }
1040 :
1041 : static const long gen_CUP_LIMIT = 5;
1042 :
1043 : static long
1044 10598 : gen_CUP(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, void *E, const struct bb_field *ff,
1045 : GEN (*mul)(void *E, GEN a, GEN))
1046 : {
1047 10598 : long m = nbrows(A), m1, n = lg(A) - 1, i, r1, r2, r;
1048 : GEN R1, C1, U1, P1, R2, C2, U2, P2;
1049 : GEN A1, A2, B2, C21, U11, U12, T21, T22;
1050 10598 : pari_sp av = avma;
1051 :
1052 10598 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1053 : /* destroy A; not called at the outermost recursion level */
1054 5985 : return gen_CUP_basecase(A, R, C, U, P, E, ff);
1055 4613 : m1 = (minss(m, n) + 1)/2;
1056 4613 : A1 = rowslice(A, 1, m1);
1057 4613 : A2 = rowslice(A, m1 + 1, m);
1058 4613 : r1 = gen_CUP(A1, &R1, &C1, &U1, &P1, E, ff, mul);
1059 4613 : if (r1 == 0)
1060 : {
1061 489 : r2 = gen_CUP(A2, &R2, &C2, &U2, &P2, E, ff, mul);
1062 489 : *R = cgetg(r2 + 1, t_VECSMALL);
1063 798 : for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
1064 489 : *C = vconcat(gen_zeromat(m1, r2, E, ff), C2);
1065 489 : *U = U2;
1066 489 : *P = P2;
1067 489 : r = r2;
1068 : }
1069 : else
1070 : {
1071 4124 : U11 = vecslice(U1, 1, r1);
1072 4124 : U12 = vecslice(U1, r1 + 1, n);
1073 4124 : T21 = vecslicepermute(A2, P1, 1, r1);
1074 4124 : T22 = vecslicepermute(A2, P1, r1 + 1, n);
1075 4124 : C21 = gen_lsolve_upper(U11, T21, E, ff, mul);
1076 4124 : if (gc_needed(av, 1))
1077 0 : gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
1078 4124 : B2 = gen_matsub(T22, mul(E, C21, U12), E, ff);
1079 4124 : r2 = gen_CUP(B2, &R2, &C2, &U2, &P2, E, ff, mul);
1080 4124 : r = r1 + r2;
1081 4124 : *R = cgetg(r + 1, t_VECSMALL);
1082 19021 : for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
1083 19879 : for ( ; i <= r; i++) (*R)[i] = R2[i - r1] + m1;
1084 4124 : *C = shallowconcat(vconcat(C1, C21),
1085 : vconcat(gen_zeromat(m1, r2, E, ff), C2));
1086 4124 : *U = shallowconcat(vconcat(U11, gen_zeromat(r2, r1, E, ff)),
1087 : vconcat(vecpermute(U12, P2), U2));
1088 :
1089 4124 : *P = cgetg(n + 1, t_VECSMALL);
1090 19021 : for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
1091 49559 : for ( ; i <= n; i++) (*P)[i] = P1[P2[i - r1] + r1];
1092 : }
1093 4613 : if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
1094 4613 : return r;
1095 : }
1096 :
1097 : /* column echelon form */
1098 : static long
1099 17685 : gen_echelon(GEN A, GEN *R, GEN *C, void *E, const struct bb_field *ff,
1100 : GEN (*mul)(void*, GEN, GEN))
1101 : {
1102 17685 : long j, j1, j2, m = nbrows(A), n = lg(A) - 1, n1, r, r1, r2;
1103 : GEN A1, A2, R1, R1c, C1, R2, C2;
1104 : GEN A12, A22, B2, C11, C21, M12;
1105 17685 : pari_sp av = avma;
1106 :
1107 17685 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1108 10022 : return gen_CUP_basecase(shallowcopy(A), R, C, NULL, NULL, E, ff);
1109 :
1110 7663 : n1 = (n + 1)/2;
1111 7663 : A1 = vecslice(A, 1, n1);
1112 7663 : A2 = vecslice(A, n1 + 1, n);
1113 7663 : r1 = gen_echelon(A1, &R1, &C1, E, ff, mul);
1114 7663 : if (!r1) return gen_echelon(A2, R, C, E, ff, mul);
1115 6781 : if (r1 == m) { *R = R1; *C = C1; return r1; }
1116 6634 : R1c = indexcompl(R1, m);
1117 6634 : C11 = rowpermute(C1, R1);
1118 6634 : C21 = rowpermute(C1, R1c);
1119 6634 : A12 = rowpermute(A2, R1);
1120 6634 : A22 = rowpermute(A2, R1c);
1121 6634 : M12 = gen_rsolve_lower_unit(C11, A12, E, ff, mul);
1122 6634 : B2 = gen_matsub(A22, mul(E, C21, M12), E, ff);
1123 6634 : r2 = gen_echelon(B2, &R2, &C2, E, ff, mul);
1124 6634 : if (!r2) { *R = R1; *C = C1; r = r1; }
1125 : else
1126 : {
1127 4350 : R2 = perm_mul(R1c, R2);
1128 4350 : C2 = rowpermute(vconcat(gen_zeromat(r1, r2, E, ff), C2),
1129 : perm_inv(vecsmall_concat(R1, R1c)));
1130 4350 : r = r1 + r2;
1131 4350 : *R = cgetg(r + 1, t_VECSMALL);
1132 4350 : *C = cgetg(r + 1, t_MAT);
1133 33176 : for (j = j1 = j2 = 1; j <= r; j++)
1134 28826 : if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
1135 : {
1136 16362 : gel(*C, j) = gel(C1, j1);
1137 16362 : (*R)[j] = R1[j1++];
1138 : }
1139 : else
1140 : {
1141 12464 : gel(*C, j) = gel(C2, j2);
1142 12464 : (*R)[j] = R2[j2++];
1143 : }
1144 : }
1145 6634 : if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
1146 6634 : return r;
1147 : }
1148 :
1149 : static GEN
1150 610 : gen_pivots_CUP(GEN x, long *rr, void *E, const struct bb_field *ff,
1151 : GEN (*mul)(void*, GEN, GEN))
1152 : {
1153 : pari_sp av;
1154 610 : long i, n = lg(x) - 1, r;
1155 610 : GEN R, C, U, P, d = zero_zv(n);
1156 610 : av = avma;
1157 610 : r = gen_CUP(x, &R, &C, &U, &P, E, ff, mul);
1158 5157 : for(i = 1; i <= r; i++)
1159 4547 : d[P[i]] = R[i];
1160 610 : set_avma(av);
1161 610 : *rr = n - r;
1162 610 : return d;
1163 : }
1164 :
1165 : static GEN
1166 140 : gen_det_CUP(GEN a, void *E, const struct bb_field *ff,
1167 : GEN (*mul)(void*, GEN, GEN))
1168 : {
1169 140 : pari_sp av = avma;
1170 : GEN R, C, U, P, d;
1171 140 : long i, n = lg(a) - 1, r;
1172 140 : r = gen_CUP(a, &R, &C, &U, &P, E, ff, mul);
1173 140 : if (r < n)
1174 0 : d = ff->s(E, 0);
1175 : else {
1176 140 : d = ff->s(E, perm_sign(P) == 1 ? 1: - 1);
1177 2730 : for (i = 1; i <= n; i++)
1178 2590 : d = ff->red(E, ff->mul(E, d, gcoeff(U, i, i)));
1179 : }
1180 140 : return gerepileupto(av, d);
1181 : }
1182 :
1183 : static long
1184 35 : gen_matrank(GEN x, void *E, const struct bb_field *ff,
1185 : GEN (*mul)(void*, GEN, GEN))
1186 : {
1187 35 : pari_sp av = avma;
1188 : long r;
1189 35 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1190 : {
1191 : GEN R, C;
1192 28 : return gc_long(av, gen_echelon(x, &R, &C, E, ff, mul));
1193 : }
1194 7 : (void) gen_Gauss_pivot(x, &r, E, ff);
1195 7 : return gc_long(av, lg(x)-1 - r);
1196 : }
1197 :
1198 : static GEN
1199 63 : gen_invimage_CUP(GEN A, GEN B, void *E, const struct bb_field *ff,
1200 : GEN (*mul)(void*, GEN, GEN))
1201 : {
1202 63 : pari_sp av = avma;
1203 : GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
1204 63 : long r = gen_CUP(A, &R, &C, &U, &P, E, ff, mul);
1205 63 : Rc = indexcompl(R, nbrows(B));
1206 63 : C1 = rowpermute(C, R);
1207 63 : C2 = rowpermute(C, Rc);
1208 63 : B1 = rowpermute(B, R);
1209 63 : B2 = rowpermute(B, Rc);
1210 63 : Z = gen_rsolve_lower_unit(C1, B1, E, ff, mul);
1211 63 : if (!gequal(mul(E, C2, Z), B2))
1212 42 : return NULL;
1213 21 : Y = vconcat(gen_rsolve_upper(vecslice(U, 1, r), Z, E, ff, mul),
1214 21 : gen_zeromat(lg(A) - 1 - r, lg(B) - 1, E, ff));
1215 21 : X = rowpermute(Y, perm_inv(P));
1216 21 : return gc_GEN(av, X);
1217 : }
1218 :
1219 : static GEN
1220 2373 : gen_ker_echelon(GEN x, void *E, const struct bb_field *ff,
1221 : GEN (*mul)(void*, GEN, GEN))
1222 : {
1223 2373 : pari_sp av = avma;
1224 : GEN R, Rc, C, C1, C2, S, K;
1225 2373 : long n = lg(x) - 1, r;
1226 2373 : r = gen_echelon(shallowtrans(x), &R, &C, E, ff, mul);
1227 2373 : Rc = indexcompl(R, n);
1228 2373 : C1 = rowpermute(C, R);
1229 2373 : C2 = rowpermute(C, Rc);
1230 2373 : S = gen_lsolve_lower_unit(C1, C2, E, ff, mul);
1231 2373 : K = vecpermute(shallowconcat(gen_matneg(S, E, ff), gen_matid(n - r, E, ff)),
1232 : perm_inv(vecsmall_concat(R, Rc)));
1233 2373 : K = shallowtrans(K);
1234 2373 : return gc_GEN(av, K);
1235 : }
1236 :
1237 : static GEN
1238 105 : gen_deplin_echelon(GEN x, void *E, const struct bb_field *ff,
1239 : GEN (*mul)(void*, GEN, GEN))
1240 : {
1241 105 : pari_sp av = avma;
1242 : GEN R, Rc, C, C1, C2, s, v;
1243 105 : long i, n = lg(x) - 1, r;
1244 105 : r = gen_echelon(shallowtrans(x), &R, &C, E, ff, mul);
1245 105 : if (r == n) return gc_NULL(av);
1246 70 : Rc = indexcompl(R, n);
1247 70 : i = Rc[1];
1248 70 : C1 = rowpermute(C, R);
1249 70 : C2 = rowslice(C, i, i);
1250 70 : s = row(gen_lsolve_lower_unit(C1, C2, E, ff, mul), 1);
1251 70 : settyp(s, t_COL);
1252 70 : v = vecpermute(shallowconcat(gen_colneg(s, E, ff), gen_colei(n - r, 1, E, ff)),
1253 : perm_inv(vecsmall_concat(R, Rc)));
1254 70 : return gc_GEN(av, v);
1255 : }
1256 :
1257 : static GEN
1258 559 : gen_gauss_CUP(GEN a, GEN b, void *E, const struct bb_field *ff,
1259 : GEN (*mul)(void*, GEN, GEN))
1260 : {
1261 : GEN R, C, U, P, Y;
1262 559 : long n = lg(a) - 1, r;
1263 559 : if (nbrows(a) < n || (r = gen_CUP(a, &R, &C, &U, &P, E, ff, mul)) < n)
1264 56 : return NULL;
1265 503 : Y = gen_rsolve_lower_unit(rowpermute(C, R), rowpermute(b, R), E, ff, mul);
1266 503 : return rowpermute(gen_rsolve_upper(U, Y, E, ff, mul), perm_inv(P));
1267 : }
1268 :
1269 : static GEN
1270 3036 : gen_gauss(GEN a, GEN b, void *E, const struct bb_field *ff,
1271 : GEN (*mul)(void*, GEN, GEN))
1272 : {
1273 3036 : if (lg(a) - 1 >= gen_CUP_LIMIT)
1274 559 : return gen_gauss_CUP(a, b, E, ff, mul);
1275 2477 : return gen_Gauss(a, b, E, ff);
1276 : }
1277 :
1278 : static GEN
1279 3672 : gen_ker_i(GEN x, long deplin, void *E, const struct bb_field *ff,
1280 : GEN (*mul)(void*, GEN, GEN)) {
1281 3672 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1282 2478 : return deplin? gen_deplin_echelon(x, E, ff, mul): gen_ker_echelon(x, E, ff, mul);
1283 1194 : return gen_ker(x, deplin, E, ff);
1284 : }
1285 :
1286 : static GEN
1287 140 : gen_invimage(GEN A, GEN B, void *E, const struct bb_field *ff,
1288 : GEN (*mul)(void*, GEN, GEN))
1289 : {
1290 140 : long nA = lg(A)-1, nB = lg(B)-1;
1291 :
1292 140 : if (!nB) return cgetg(1, t_MAT);
1293 140 : if (nA + nB >= gen_CUP_LIMIT && nbrows(B) >= gen_CUP_LIMIT)
1294 63 : return gen_invimage_CUP(A, B, E, ff, mul);
1295 77 : return gen_matinvimage(A, B, E, ff);
1296 : }
1297 :
1298 : /* find z such that A z = y. Return NULL if no solution */
1299 : static GEN
1300 71 : gen_matcolinvimage_i(GEN A, GEN y, void *E, const struct bb_field *ff,
1301 : GEN (*mul)(void*, GEN, GEN))
1302 : {
1303 71 : pari_sp av = avma;
1304 71 : long i, l = lg(A);
1305 : GEN M, x, t;
1306 :
1307 71 : M = gen_ker_i(shallowconcat(A, y), 0, E, ff, mul);
1308 71 : i = lg(M) - 1;
1309 71 : if (!i) return gc_NULL(av);
1310 :
1311 71 : x = gel(M, i);
1312 71 : t = gel(x, l);
1313 71 : if (ff->equal0(t)) return gc_NULL(av);
1314 :
1315 50 : t = ff->neg(E, ff->inv(E, t));
1316 50 : setlg(x, l);
1317 178 : for (i = 1; i < l; i++)
1318 128 : gel(x, i) = ff->red(E, ff->mul(E, t, gel(x, i)));
1319 50 : return gc_GEN(av, x);
1320 : }
1321 :
1322 : static GEN
1323 420 : gen_det_i(GEN a, void *E, const struct bb_field *ff,
1324 : GEN (*mul)(void*, GEN, GEN))
1325 : {
1326 420 : if (lg(a) - 1 >= gen_CUP_LIMIT)
1327 140 : return gen_det_CUP(a, E, ff, mul);
1328 : else
1329 280 : return gen_det(a, E, ff);
1330 : }
1331 :
1332 : static GEN
1333 1722 : gen_pivots(GEN x, long *rr, void *E, const struct bb_field *ff,
1334 : GEN (*mul)(void*, GEN, GEN))
1335 : {
1336 1722 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1337 610 : return gen_pivots_CUP(x, rr, E, ff, mul);
1338 1112 : return gen_Gauss_pivot(x, rr, E, ff);
1339 : }
1340 :
1341 : /* r = dim Ker x, n = nbrows(x) */
1342 : static GEN
1343 21 : gen_get_suppl(GEN x, GEN d, long n, long r, void *E, const struct bb_field *ff)
1344 : {
1345 : GEN y, c;
1346 21 : long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
1347 :
1348 21 : if (rx == n && r == 0) return gcopy(x);
1349 21 : c = zero_zv(n);
1350 21 : y = cgetg(n+1, t_MAT);
1351 : /* c = lines containing pivots (could get it from RgM_pivots, but cheap)
1352 : * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
1353 119 : for (k = j = 1; j<=rx; j++)
1354 98 : if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gcopy(gel(x,j)); }
1355 203 : for (j=1; j<=n; j++)
1356 182 : if (!c[j]) gel(y,k++) = gen_colei(n, j, E, ff);
1357 21 : return y;
1358 : }
1359 :
1360 : static GEN
1361 21 : gen_suppl(GEN x, void *E, const struct bb_field *ff,
1362 : GEN (*mul)(void*, GEN, GEN))
1363 : {
1364 : GEN d;
1365 21 : long n = nbrows(x), r;
1366 :
1367 21 : if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
1368 21 : d = gen_pivots(x, &r, E, ff, mul);
1369 21 : return gen_get_suppl(x, d, n, r, E, ff);
1370 : }
1371 :
1372 : /*******************************************************************/
1373 : /* */
1374 : /* MATRIX MULTIPLICATION MODULO P */
1375 : /* */
1376 : /*******************************************************************/
1377 :
1378 : GEN
1379 21 : F2xqM_F2xqC_mul(GEN A, GEN B, GEN T) {
1380 : void *E;
1381 21 : const struct bb_field *ff = get_F2xq_field(&E, T);
1382 21 : return gen_matcolmul(A, B, E, ff);
1383 : }
1384 :
1385 : GEN
1386 35 : FlxqM_FlxqC_mul(GEN A, GEN B, GEN T, ulong p) {
1387 : void *E;
1388 35 : const struct bb_field *ff = get_Flxq_field(&E, T, p);
1389 35 : return gen_matcolmul(A, B, E, ff);
1390 : }
1391 :
1392 : GEN
1393 63 : FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
1394 : void *E;
1395 63 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1396 63 : return gen_matcolmul(A, B, E, ff);
1397 : }
1398 :
1399 : GEN
1400 1449 : F2xqM_mul(GEN A, GEN B, GEN T) {
1401 : void *E;
1402 1449 : const struct bb_field *ff = get_F2xq_field(&E, T);
1403 1449 : return gen_matmul(A, B, E, ff);
1404 : }
1405 :
1406 : GEN
1407 148988 : FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
1408 : void *E;
1409 : const struct bb_field *ff;
1410 148988 : long n = lg(A) - 1;
1411 :
1412 148988 : if (n == 0)
1413 0 : return cgetg(1, t_MAT);
1414 148988 : if (n > 1)
1415 81581 : return FlxqM_mul_Kronecker(A, B, T, p);
1416 67407 : ff = get_Flxq_field(&E, T, p);
1417 67407 : return gen_matmul(A, B, E, ff);
1418 : }
1419 :
1420 : GEN
1421 86016 : FqM_mul(GEN A, GEN B, GEN T, GEN p) {
1422 : void *E;
1423 86016 : long n = lg(A) - 1;
1424 : const struct bb_field *ff;
1425 86016 : if (n == 0)
1426 0 : return cgetg(1, t_MAT);
1427 86016 : if (n > 1)
1428 81851 : return FqM_mul_Kronecker(A, B, T, p);
1429 4165 : ff = get_Fq_field(&E, T, p);
1430 4165 : return gen_matmul(A, B, E, ff);
1431 : }
1432 :
1433 : /*******************************************************************/
1434 : /* */
1435 : /* LINEAR ALGEBRA MODULO P */
1436 : /* */
1437 : /*******************************************************************/
1438 :
1439 : static GEN
1440 0 : _F2xqM_mul(void *E, GEN A, GEN B)
1441 0 : { return F2xqM_mul(A, B, (GEN) E); }
1442 :
1443 : struct _Flxq {
1444 : GEN aut;
1445 : GEN T;
1446 : ulong p;
1447 : };
1448 :
1449 : static GEN
1450 7924 : _FlxqM_mul(void *E, GEN A, GEN B)
1451 : {
1452 7924 : struct _Flxq *D = (struct _Flxq*)E;
1453 7924 : return FlxqM_mul(A, B, D->T, D->p);
1454 : }
1455 :
1456 : static GEN
1457 22487 : _FpM_mul(void *E, GEN A, GEN B)
1458 22487 : { return FpM_mul(A, B, (GEN) E); }
1459 :
1460 : struct _Fq_field
1461 : {
1462 : GEN T, p;
1463 : };
1464 :
1465 : static GEN
1466 6349 : _FqM_mul(void *E, GEN A, GEN B)
1467 : {
1468 6349 : struct _Fq_field *D = (struct _Fq_field*) E;
1469 6349 : return FqM_mul(A, B, D->T, D->p);
1470 : }
1471 :
1472 : static GEN
1473 1289600 : FpM_init(GEN a, GEN p, ulong *pp)
1474 : {
1475 1289600 : if (lgefint(p) == 3)
1476 : {
1477 1285316 : *pp = uel(p,2);
1478 1285316 : return (*pp==2)? ZM_to_F2m(a): ZM_to_Flm(a, *pp);
1479 : }
1480 4284 : *pp = 0; return a;
1481 : }
1482 : static GEN
1483 1308678 : FpM_init3(GEN a, GEN p, ulong *pp)
1484 : {
1485 1308678 : if (lgefint(p) == 3)
1486 : {
1487 1306107 : *pp = uel(p,2);
1488 1306107 : switch(*pp)
1489 : {
1490 705315 : case 2: return ZM_to_F2m(a);
1491 157134 : case 3: return ZM_to_F3m(a);
1492 443658 : default:return ZM_to_Flm(a, *pp);
1493 : }
1494 : }
1495 2571 : *pp = 0; return a;
1496 : }
1497 : GEN
1498 4599 : RgM_Fp_init(GEN a, GEN p, ulong *pp)
1499 : {
1500 4599 : if (lgefint(p) == 3)
1501 : {
1502 4319 : *pp = uel(p,2);
1503 4319 : return (*pp==2)? RgM_to_F2m(a): RgM_to_Flm(a, *pp);
1504 : }
1505 280 : *pp = 0; return RgM_to_FpM(a,p);
1506 : }
1507 : static GEN
1508 658 : RgM_Fp_init3(GEN a, GEN p, ulong *pp)
1509 : {
1510 658 : if (lgefint(p) == 3)
1511 : {
1512 588 : *pp = uel(p,2);
1513 588 : switch(*pp)
1514 : {
1515 35 : case 2: return RgM_to_F2m(a);
1516 77 : case 3: return RgM_to_F3m(a);
1517 476 : default:return RgM_to_Flm(a, *pp);
1518 : }
1519 : }
1520 70 : *pp = 0; return RgM_to_FpM(a,p);
1521 : }
1522 :
1523 : static GEN
1524 315 : FpM_det_gen(GEN a, GEN p)
1525 : {
1526 : void *E;
1527 315 : const struct bb_field *S = get_Fp_field(&E,p);
1528 315 : return gen_det_i(a, E, S, _FpM_mul);
1529 : }
1530 : GEN
1531 4683 : FpM_det(GEN a, GEN p)
1532 : {
1533 4683 : pari_sp av = avma;
1534 : ulong pp, d;
1535 4683 : a = FpM_init(a, p, &pp);
1536 4683 : switch(pp)
1537 : {
1538 315 : case 0: return FpM_det_gen(a, p);
1539 1750 : case 2: d = F2m_det_sp(a); break;
1540 2618 : default:d = Flm_det_sp(a,pp); break;
1541 : }
1542 4368 : return gc_utoi(av, d);
1543 : }
1544 :
1545 : GEN
1546 7 : F2xqM_det(GEN a, GEN T)
1547 : {
1548 : void *E;
1549 7 : const struct bb_field *S = get_F2xq_field(&E, T);
1550 7 : return gen_det_i(a, E, S, _F2xqM_mul);
1551 : }
1552 :
1553 : GEN
1554 28 : FlxqM_det(GEN a, GEN T, ulong p) {
1555 : void *E;
1556 28 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1557 28 : return gen_det_i(a, E, S, _FlxqM_mul);
1558 : }
1559 :
1560 : GEN
1561 70 : FqM_det(GEN x, GEN T, GEN p)
1562 : {
1563 : void *E;
1564 70 : const struct bb_field *S = get_Fq_field(&E,T,p);
1565 70 : return gen_det_i(x, E, S, _FqM_mul);
1566 : }
1567 :
1568 : static GEN
1569 1214 : FpM_gauss_pivot_gen(GEN x, GEN p, long *rr)
1570 : {
1571 : void *E;
1572 1214 : const struct bb_field *S = get_Fp_field(&E,p);
1573 1214 : return gen_pivots(x, rr, E, S, _FpM_mul);
1574 : }
1575 :
1576 : static GEN
1577 640609 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
1578 : {
1579 : ulong pp;
1580 640609 : if (lg(x)==1) { *rr = 0; return NULL; }
1581 635533 : x = FpM_init(x, p, &pp);
1582 635540 : switch(pp)
1583 : {
1584 1214 : case 0: return FpM_gauss_pivot_gen(x, p, rr);
1585 352685 : case 2: return F2m_gauss_pivot(x, rr);
1586 281641 : default:return Flm_pivots(x, pp, rr, 1);
1587 : }
1588 : }
1589 :
1590 : static GEN
1591 21 : F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
1592 : {
1593 : void *E;
1594 21 : const struct bb_field *S = get_F2xq_field(&E,T);
1595 21 : return gen_pivots(x, rr, E, S, _F2xqM_mul);
1596 : }
1597 :
1598 : static GEN
1599 361 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr) {
1600 : void *E;
1601 361 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1602 361 : return gen_pivots(x, rr, E, S, _FlxqM_mul);
1603 : }
1604 :
1605 : static GEN
1606 105 : FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
1607 : {
1608 : void *E;
1609 105 : const struct bb_field *S = get_Fq_field(&E,T,p);
1610 105 : return gen_pivots(x, rr, E, S, _FqM_mul);
1611 : }
1612 : static GEN
1613 438 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
1614 : {
1615 438 : if (lg(x)==1) { *rr = 0; return NULL; }
1616 438 : if (!T) return FpM_gauss_pivot(x, p, rr);
1617 438 : if (lgefint(p) == 3)
1618 : {
1619 333 : pari_sp av = avma;
1620 333 : ulong pp = uel(p,2);
1621 333 : GEN Tp = ZXT_to_FlxT(T, pp);
1622 333 : GEN d = FlxqM_gauss_pivot(ZXM_to_FlxM(x, pp, get_Flx_var(Tp)), Tp, pp, rr);
1623 333 : return d ? gerepileuptoleaf(av, d): d;
1624 : }
1625 105 : return FqM_gauss_pivot_gen(x, T, p, rr);
1626 : }
1627 :
1628 : GEN
1629 330532 : FpM_image(GEN x, GEN p)
1630 : {
1631 : long r;
1632 330532 : GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
1633 330534 : return image_from_pivot(x,d,r);
1634 : }
1635 :
1636 : GEN
1637 40859 : Flm_image(GEN x, ulong p)
1638 : {
1639 : long r;
1640 40859 : GEN d = Flm_pivots(x, p, &r, 0); /* d left on stack for efficiency */
1641 40859 : return image_from_pivot(x,d,r);
1642 : }
1643 :
1644 : GEN
1645 7 : F2m_image(GEN x)
1646 : {
1647 : long r;
1648 7 : GEN d = F2m_gauss_pivot(F2m_copy(x),&r); /* d left on stack for efficiency */
1649 7 : return image_from_pivot(x,d,r);
1650 : }
1651 :
1652 : GEN
1653 7 : F2xqM_image(GEN x, GEN T)
1654 : {
1655 : long r;
1656 7 : GEN d = F2xqM_gauss_pivot(x,T,&r); /* d left on stack for efficiency */
1657 7 : return image_from_pivot(x,d,r);
1658 : }
1659 :
1660 : GEN
1661 21 : FlxqM_image(GEN x, GEN T, ulong p)
1662 : {
1663 : long r;
1664 21 : GEN d = FlxqM_gauss_pivot(x, T, p, &r); /* d left on stack for efficiency */
1665 21 : return image_from_pivot(x,d,r);
1666 : }
1667 :
1668 : GEN
1669 49 : FqM_image(GEN x, GEN T, GEN p)
1670 : {
1671 : long r;
1672 49 : GEN d = FqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
1673 49 : return image_from_pivot(x,d,r);
1674 : }
1675 :
1676 : long
1677 56 : FpM_rank(GEN x, GEN p)
1678 : {
1679 56 : pari_sp av = avma;
1680 : long r;
1681 56 : (void)FpM_gauss_pivot(x,p,&r);
1682 56 : return gc_long(av, lg(x)-1 - r);
1683 : }
1684 :
1685 : long
1686 7 : F2xqM_rank(GEN x, GEN T)
1687 : {
1688 7 : pari_sp av = avma;
1689 : long r;
1690 7 : (void)F2xqM_gauss_pivot(x,T,&r);
1691 7 : return gc_long(av, lg(x)-1 - r);
1692 : }
1693 :
1694 : long
1695 35 : FlxqM_rank(GEN x, GEN T, ulong p)
1696 : {
1697 : void *E;
1698 35 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1699 35 : return gen_matrank(x, E, S, _FlxqM_mul);
1700 : }
1701 :
1702 : long
1703 70 : FqM_rank(GEN x, GEN T, GEN p)
1704 : {
1705 70 : pari_sp av = avma;
1706 : long r;
1707 70 : (void)FqM_gauss_pivot(x,T,p,&r);
1708 70 : return gc_long(av, lg(x)-1 - r);
1709 : }
1710 :
1711 : static GEN
1712 35 : FpM_invimage_gen(GEN A, GEN B, GEN p)
1713 : {
1714 : void *E;
1715 35 : const struct bb_field *ff = get_Fp_field(&E, p);
1716 35 : return gen_invimage(A, B, E, ff, _FpM_mul);
1717 : }
1718 :
1719 : GEN
1720 0 : FpM_invimage(GEN A, GEN B, GEN p)
1721 : {
1722 0 : pari_sp av = avma;
1723 : ulong pp;
1724 : GEN y;
1725 :
1726 0 : A = FpM_init(A, p, &pp);
1727 0 : switch(pp)
1728 : {
1729 0 : case 0: return FpM_invimage_gen(A, B, p);
1730 0 : case 2:
1731 0 : y = F2m_invimage(A, ZM_to_F2m(B));
1732 0 : if (!y) return gc_NULL(av);
1733 0 : y = F2m_to_ZM(y);
1734 0 : return gerepileupto(av, y);
1735 0 : default:
1736 0 : y = Flm_invimage_i(A, ZM_to_Flm(B, pp), pp);
1737 0 : if (!y) return gc_NULL(av);
1738 0 : y = Flm_to_ZM(y);
1739 0 : return gerepileupto(av, y);
1740 : }
1741 : }
1742 :
1743 : GEN
1744 21 : F2xqM_invimage(GEN A, GEN B, GEN T) {
1745 : void *E;
1746 21 : const struct bb_field *ff = get_F2xq_field(&E, T);
1747 21 : return gen_invimage(A, B, E, ff, _F2xqM_mul);
1748 : }
1749 :
1750 : GEN
1751 42 : FlxqM_invimage(GEN A, GEN B, GEN T, ulong p) {
1752 : void *E;
1753 42 : const struct bb_field *ff = get_Flxq_field(&E, T, p);
1754 42 : return gen_invimage(A, B, E, ff, _FlxqM_mul);
1755 : }
1756 :
1757 : GEN
1758 42 : FqM_invimage(GEN A, GEN B, GEN T, GEN p) {
1759 : void *E;
1760 42 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1761 42 : return gen_invimage(A, B, E, ff, _FqM_mul);
1762 : }
1763 :
1764 : static GEN
1765 8 : FpM_FpC_invimage_gen(GEN A, GEN y, GEN p)
1766 : {
1767 : void *E;
1768 8 : const struct bb_field *ff = get_Fp_field(&E, p);
1769 8 : return gen_matcolinvimage_i(A, y, E, ff, _FpM_mul);
1770 : }
1771 :
1772 : GEN
1773 297545 : FpM_FpC_invimage(GEN A, GEN x, GEN p)
1774 : {
1775 297545 : pari_sp av = avma;
1776 : ulong pp;
1777 : GEN y;
1778 :
1779 297545 : A = FpM_init(A, p, &pp);
1780 297556 : switch(pp)
1781 : {
1782 8 : case 0: return FpM_FpC_invimage_gen(A, x, p);
1783 193256 : case 2:
1784 193256 : y = F2m_F2c_invimage(A, ZV_to_F2v(x));
1785 193253 : if (!y) return y;
1786 193253 : y = F2c_to_ZC(y);
1787 193254 : return gerepileupto(av, y);
1788 104292 : default:
1789 104292 : y = Flm_Flc_invimage(A, ZV_to_Flv(x, pp), pp);
1790 104292 : if (!y) return y;
1791 104292 : y = Flc_to_ZC(y);
1792 104292 : return gerepileupto(av, y);
1793 : }
1794 : }
1795 :
1796 : GEN
1797 21 : F2xqM_F2xqC_invimage(GEN A, GEN B, GEN T) {
1798 : void *E;
1799 21 : const struct bb_field *ff = get_F2xq_field(&E, T);
1800 21 : return gen_matcolinvimage_i(A, B, E, ff, _F2xqM_mul);
1801 : }
1802 :
1803 : GEN
1804 21 : FlxqM_FlxqC_invimage(GEN A, GEN B, GEN T, ulong p) {
1805 : void *E;
1806 21 : const struct bb_field *ff = get_Flxq_field(&E, T, p);
1807 21 : return gen_matcolinvimage_i(A, B, E, ff, _FlxqM_mul);
1808 : }
1809 :
1810 : GEN
1811 21 : FqM_FqC_invimage(GEN A, GEN B, GEN T, GEN p) {
1812 : void *E;
1813 21 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1814 21 : return gen_matcolinvimage_i(A, B, E, ff, _FqM_mul);
1815 : }
1816 :
1817 : static GEN
1818 2642 : FpM_ker_gen(GEN x, GEN p, long deplin)
1819 : {
1820 : void *E;
1821 2642 : const struct bb_field *S = get_Fp_field(&E,p);
1822 2642 : return gen_ker_i(x, deplin, E, S, _FpM_mul);
1823 : }
1824 : static GEN
1825 1308679 : FpM_ker_i(GEN x, GEN p, long deplin)
1826 : {
1827 1308679 : pari_sp av = avma;
1828 : ulong pp;
1829 : GEN y;
1830 :
1831 1308679 : if (lg(x)==1) return cgetg(1,t_MAT);
1832 1308679 : x = FpM_init3(x, p, &pp);
1833 1308714 : switch(pp)
1834 : {
1835 2572 : case 0: return FpM_ker_gen(x,p,deplin);
1836 705343 : case 2:
1837 705343 : y = F2m_ker_sp(x, deplin);
1838 705330 : if (!y) return gc_NULL(av);
1839 705337 : y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
1840 705340 : return gerepileupto(av, y);
1841 157134 : case 3:
1842 157134 : y = F3m_ker_sp(x, deplin);
1843 157134 : if (!y) return gc_NULL(av);
1844 157134 : y = deplin? F3c_to_ZC(y): F3m_to_ZM(y);
1845 157134 : return gerepileupto(av, y);
1846 443665 : default:
1847 443665 : y = Flm_ker_sp(x, pp, deplin);
1848 443665 : if (!y) return gc_NULL(av);
1849 443665 : y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
1850 443665 : return gerepileupto(av, y);
1851 : }
1852 : }
1853 :
1854 : GEN
1855 849698 : FpM_ker(GEN x, GEN p) { return FpM_ker_i(x,p,0); }
1856 :
1857 : static GEN
1858 21 : F2xqM_ker_i(GEN x, GEN T, long deplin)
1859 : {
1860 : const struct bb_field *ff;
1861 : void *E;
1862 :
1863 21 : if (lg(x)==1) return cgetg(1,t_MAT);
1864 21 : ff = get_F2xq_field(&E,T);
1865 21 : return gen_ker_i(x,deplin, E, ff, _F2xqM_mul);
1866 : }
1867 :
1868 : GEN
1869 7 : F2xqM_ker(GEN x, GEN T)
1870 : {
1871 7 : return F2xqM_ker_i(x, T, 0);
1872 : }
1873 :
1874 : static GEN
1875 812 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin) {
1876 : void *E;
1877 812 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1878 812 : return gen_ker_i(x, deplin, E, S, _FlxqM_mul);
1879 : }
1880 :
1881 : GEN
1882 28 : FlxqM_ker(GEN x, GEN T, ulong p)
1883 : {
1884 28 : return FlxqM_ker_i(x, T, p, 0);
1885 : }
1886 :
1887 : static GEN
1888 126 : FqM_ker_gen(GEN x, GEN T, GEN p, long deplin)
1889 : {
1890 : void *E;
1891 126 : const struct bb_field *S = get_Fq_field(&E,T,p);
1892 126 : return gen_ker_i(x,deplin,E,S,_FqM_mul);
1893 : }
1894 : static GEN
1895 3535 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
1896 : {
1897 3535 : if (!T) return FpM_ker_i(x,p,deplin);
1898 875 : if (lg(x)==1) return cgetg(1,t_MAT);
1899 :
1900 875 : if (lgefint(p)==3)
1901 : {
1902 749 : pari_sp av = avma;
1903 749 : ulong l = p[2];
1904 749 : GEN Tl = ZXT_to_FlxT(T,l);
1905 749 : GEN Ml = ZXM_to_FlxM(x, l, get_Flx_var(Tl));
1906 749 : GEN K = FlxqM_ker_i(Ml, Tl, l, deplin);
1907 749 : if (!deplin) K = FlxM_to_ZXM(K);
1908 28 : else if (!K) return gc_NULL(av);
1909 21 : else K = FlxC_to_ZXC(K);
1910 742 : return gerepileupto(av, K);
1911 : }
1912 126 : return FqM_ker_gen(x, T, p, deplin);
1913 : }
1914 :
1915 : GEN
1916 3451 : FqM_ker(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,0); }
1917 :
1918 : GEN
1919 456337 : FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x,p,1); }
1920 :
1921 : GEN
1922 14 : F2xqM_deplin(GEN x, GEN T)
1923 : {
1924 14 : return F2xqM_ker_i(x, T, 1);
1925 : }
1926 :
1927 : GEN
1928 35 : FlxqM_deplin(GEN x, GEN T, ulong p)
1929 : {
1930 35 : return FlxqM_ker_i(x, T, p, 1);
1931 : }
1932 :
1933 : GEN
1934 84 : FqM_deplin(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,1); }
1935 :
1936 : static GEN
1937 2749 : FpM_gauss_gen(GEN a, GEN b, GEN p)
1938 : {
1939 : void *E;
1940 2749 : const struct bb_field *S = get_Fp_field(&E,p);
1941 2749 : return gen_gauss(a,b, E, S, _FpM_mul);
1942 : }
1943 : /* a an FpM, lg(a)>1; b an FpM or NULL (replace by identity) */
1944 : static GEN
1945 351867 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
1946 : {
1947 351867 : long n = nbrows(a);
1948 351867 : a = FpM_init(a,p,pp);
1949 351865 : switch(*pp)
1950 : {
1951 2749 : case 0:
1952 2749 : if (!b) b = matid(n);
1953 2749 : return FpM_gauss_gen(a,b,p);
1954 228720 : case 2:
1955 228720 : if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
1956 228721 : return F2m_gauss_sp(a,b);
1957 120396 : default:
1958 120396 : if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
1959 120396 : return Flm_gauss_sp(a,b, NULL, *pp);
1960 : }
1961 : }
1962 : GEN
1963 35 : FpM_gauss(GEN a, GEN b, GEN p)
1964 : {
1965 35 : pari_sp av = avma;
1966 : ulong pp;
1967 : GEN u;
1968 35 : if (lg(a) == 1 || lg(b)==1) return cgetg(1, t_MAT);
1969 35 : u = FpM_gauss_i(a, b, p, &pp);
1970 35 : if (!u) return gc_NULL(av);
1971 28 : switch(pp)
1972 : {
1973 28 : case 0: return gc_GEN(av, u);
1974 0 : case 2: u = F2m_to_ZM(u); break;
1975 0 : default: u = Flm_to_ZM(u); break;
1976 : }
1977 0 : return gerepileupto(av, u);
1978 : }
1979 :
1980 : static GEN
1981 63 : F2xqM_gauss_gen(GEN a, GEN b, GEN T)
1982 : {
1983 : void *E;
1984 63 : const struct bb_field *S = get_F2xq_field(&E, T);
1985 63 : return gen_gauss(a, b, E, S, _F2xqM_mul);
1986 : }
1987 :
1988 : GEN
1989 14 : F2xqM_gauss(GEN a, GEN b, GEN T)
1990 : {
1991 14 : pari_sp av = avma;
1992 14 : long n = lg(a)-1;
1993 : GEN u;
1994 14 : if (!n || lg(b)==1) retgc_const(av, cgetg(1, t_MAT));
1995 14 : u = F2xqM_gauss_gen(a, b, T);
1996 14 : if (!u) return gc_NULL(av);
1997 14 : return gc_GEN(av, u);
1998 : }
1999 :
2000 : static GEN
2001 91 : FlxqM_gauss_i(GEN a, GEN b, GEN T, ulong p) {
2002 : void *E;
2003 91 : const struct bb_field *S = get_Flxq_field(&E, T, p);
2004 91 : return gen_gauss(a, b, E, S, _FlxqM_mul);
2005 : }
2006 :
2007 : GEN
2008 21 : FlxqM_gauss(GEN a, GEN b, GEN T, ulong p)
2009 : {
2010 21 : pari_sp av = avma;
2011 21 : long n = lg(a)-1;
2012 : GEN u;
2013 21 : if (!n || lg(b)==1) retgc_const(av, cgetg(1, t_MAT));
2014 21 : u = FlxqM_gauss_i(a, b, T, p);
2015 21 : if (!u) return gc_NULL(av);
2016 14 : return gc_GEN(av, u);
2017 : }
2018 :
2019 : static GEN
2020 133 : FqM_gauss_gen(GEN a, GEN b, GEN T, GEN p)
2021 : {
2022 : void *E;
2023 133 : const struct bb_field *S = get_Fq_field(&E,T,p);
2024 133 : return gen_gauss(a,b,E,S,_FqM_mul);
2025 : }
2026 : GEN
2027 21 : FqM_gauss(GEN a, GEN b, GEN T, GEN p)
2028 : {
2029 21 : pari_sp av = avma;
2030 : GEN u;
2031 : long n;
2032 21 : if (!T) return FpM_gauss(a,b,p);
2033 21 : n = lg(a)-1; if (!n || lg(b)==1) return cgetg(1, t_MAT);
2034 21 : u = FqM_gauss_gen(a,b,T,p);
2035 21 : if (!u) return gc_NULL(av);
2036 14 : return gc_GEN(av, u);
2037 : }
2038 :
2039 : GEN
2040 14 : FpM_FpC_gauss(GEN a, GEN b, GEN p)
2041 : {
2042 14 : pari_sp av = avma;
2043 : ulong pp;
2044 : GEN u;
2045 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2046 14 : u = FpM_gauss_i(a, mkmat(b), p, &pp);
2047 14 : if (!u) return gc_NULL(av);
2048 14 : switch(pp)
2049 : {
2050 14 : case 0: return gc_GEN(av, gel(u,1));
2051 0 : case 2: u = F2c_to_ZC(gel(u,1)); break;
2052 0 : default: u = Flc_to_ZC(gel(u,1)); break;
2053 : }
2054 0 : return gerepileupto(av, u);
2055 : }
2056 :
2057 : GEN
2058 14 : F2xqM_F2xqC_gauss(GEN a, GEN b, GEN T)
2059 : {
2060 14 : pari_sp av = avma;
2061 : GEN u;
2062 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2063 14 : u = F2xqM_gauss_gen(a, mkmat(b), T);
2064 14 : if (!u) return gc_NULL(av);
2065 7 : return gc_GEN(av, gel(u,1));
2066 : }
2067 :
2068 : GEN
2069 14 : FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p)
2070 : {
2071 14 : pari_sp av = avma;
2072 : GEN u;
2073 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2074 14 : u = FlxqM_gauss_i(a, mkmat(b), T, p);
2075 14 : if (!u) return gc_NULL(av);
2076 7 : return gc_GEN(av, gel(u,1));
2077 : }
2078 :
2079 : GEN
2080 14 : FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p)
2081 : {
2082 14 : pari_sp av = avma;
2083 : GEN u;
2084 14 : if (!T) return FpM_FpC_gauss(a,b,p);
2085 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2086 14 : u = FqM_gauss_gen(a,mkmat(b),T,p);
2087 14 : if (!u) return gc_NULL(av);
2088 7 : return gc_GEN(av, gel(u,1));
2089 : }
2090 :
2091 : GEN
2092 351818 : FpM_inv(GEN a, GEN p)
2093 : {
2094 351818 : pari_sp av = avma;
2095 : ulong pp;
2096 : GEN u;
2097 351818 : if (lg(a) == 1) return cgetg(1, t_MAT);
2098 351818 : u = FpM_gauss_i(a, NULL, p, &pp);
2099 351819 : if (!u) return gc_NULL(av);
2100 351805 : switch(pp)
2101 : {
2102 2693 : case 0: return gc_GEN(av, u);
2103 228716 : case 2: u = F2m_to_ZM(u); break;
2104 120396 : default: u = Flm_to_ZM(u); break;
2105 : }
2106 349108 : return gerepileupto(av, u);
2107 : }
2108 :
2109 : GEN
2110 35 : F2xqM_inv(GEN a, GEN T)
2111 : {
2112 35 : pari_sp av = avma;
2113 : GEN u;
2114 35 : if (lg(a) == 1) retgc_const(av, cgetg(1, t_MAT));
2115 35 : u = F2xqM_gauss_gen(a, matid_F2xqM(nbrows(a),T), T);
2116 35 : if (!u) return gc_NULL(av);
2117 28 : return gc_GEN(av, u);
2118 : }
2119 :
2120 : GEN
2121 56 : FlxqM_inv(GEN a, GEN T, ulong p)
2122 : {
2123 56 : pari_sp av = avma;
2124 : GEN u;
2125 56 : if (lg(a) == 1) retgc_const(av, cgetg(1, t_MAT));
2126 56 : u = FlxqM_gauss_i(a, matid_FlxqM(nbrows(a),T,p), T,p);
2127 56 : if (!u) return gc_NULL(av);
2128 42 : return gc_GEN(av, u);
2129 : }
2130 :
2131 : GEN
2132 98 : FqM_inv(GEN a, GEN T, GEN p)
2133 : {
2134 98 : pari_sp av = avma;
2135 : GEN u;
2136 98 : if (!T) return FpM_inv(a,p);
2137 98 : if (lg(a) == 1) return cgetg(1, t_MAT);
2138 98 : u = FqM_gauss_gen(a,matid(nbrows(a)),T,p);
2139 98 : if (!u) return gc_NULL(av);
2140 70 : return gc_GEN(av, u);
2141 : }
2142 :
2143 : GEN
2144 263894 : FpM_intersect_i(GEN x, GEN y, GEN p)
2145 : {
2146 263894 : long j, lx = lg(x);
2147 : GEN z;
2148 :
2149 263894 : if (lx == 1 || lg(y) == 1) return cgetg(1,t_MAT);
2150 263894 : if (lgefint(p) == 3)
2151 : {
2152 263894 : ulong pp = p[2];
2153 263894 : return Flm_to_ZM(Flm_intersect_i(ZM_to_Flm(x,pp), ZM_to_Flm(y,pp), pp));
2154 : }
2155 0 : z = FpM_ker(shallowconcat(x,y), p);
2156 0 : for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
2157 0 : return FpM_mul(x,z,p);
2158 : }
2159 : GEN
2160 0 : FpM_intersect(GEN x, GEN y, GEN p)
2161 : {
2162 0 : pari_sp av = avma;
2163 : GEN z;
2164 0 : if (lgefint(p) == 3)
2165 : {
2166 0 : ulong pp = p[2];
2167 0 : z = Flm_to_ZM(Flm_image(Flm_intersect_i(ZM_to_Flm(x,pp), ZM_to_Flm(y,pp), pp), pp));
2168 : }
2169 : else
2170 0 : z = FpM_image(FpM_intersect_i(x,y,p), p);
2171 0 : return gerepileupto(av, z);
2172 : }
2173 :
2174 : /* HACK: avoid overwriting d from RgM_pivots after set_avma(av) in suppl
2175 : * or indexrank-type functions */
2176 : static void
2177 267804 : init_suppl(GEN x)
2178 : {
2179 267804 : if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
2180 267804 : (void)new_chunk(lgcols(x) * 2);
2181 267803 : }
2182 : static void
2183 1925365 : init_pivot_list(GEN x) { (void)new_chunk(3 + 2*lg(x)); /* HACK */ }
2184 :
2185 : GEN
2186 267303 : FpM_suppl(GEN x, GEN p)
2187 : {
2188 : GEN d;
2189 : long r;
2190 267303 : init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
2191 267306 : return get_suppl(x,d,nbrows(x),r,&col_ei);
2192 : }
2193 :
2194 : GEN
2195 14 : F2m_suppl(GEN x)
2196 : {
2197 : GEN d;
2198 : long r;
2199 14 : init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
2200 14 : return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
2201 : }
2202 :
2203 : GEN
2204 105 : Flm_suppl(GEN x, ulong p)
2205 : {
2206 : GEN d;
2207 : long r;
2208 105 : init_suppl(x); d = Flm_pivots(x, p, &r, 0);
2209 105 : return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
2210 : }
2211 :
2212 : GEN
2213 7 : F2xqM_suppl(GEN x, GEN T)
2214 : {
2215 : void *E;
2216 7 : const struct bb_field *S = get_F2xq_field(&E, T);
2217 7 : return gen_suppl(x, E, S, _F2xqM_mul);
2218 : }
2219 :
2220 : GEN
2221 14 : FlxqM_suppl(GEN x, GEN T, ulong p)
2222 : {
2223 : void *E;
2224 14 : const struct bb_field *S = get_Flxq_field(&E, T, p);
2225 14 : return gen_suppl(x, E, S, _FlxqM_mul);
2226 : }
2227 :
2228 : GEN
2229 1481 : FqM_suppl(GEN x, GEN T, GEN p)
2230 : {
2231 1481 : pari_sp av = avma;
2232 : GEN d;
2233 : long r;
2234 :
2235 1481 : if (!T) return FpM_suppl(x,p);
2236 312 : init_suppl(x);
2237 312 : d = FqM_gauss_pivot(x,T,p,&r);
2238 312 : set_avma(av); return get_suppl(x,d,nbrows(x),r,&col_ei);
2239 : }
2240 :
2241 : GEN
2242 42726 : FpM_indexrank(GEN x, GEN p) {
2243 42726 : pari_sp av = avma;
2244 : long r;
2245 : GEN d;
2246 42726 : init_pivot_list(x);
2247 42726 : d = FpM_gauss_pivot(x,p,&r);
2248 42726 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2249 : }
2250 :
2251 : GEN
2252 58309 : Flm_indexrank(GEN x, ulong p) {
2253 58309 : pari_sp av = avma;
2254 : long r;
2255 : GEN d;
2256 58309 : init_pivot_list(x);
2257 58309 : d = Flm_pivots(x, p, &r, 0);
2258 58310 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2259 : }
2260 :
2261 : GEN
2262 53 : F2m_indexrank(GEN x) {
2263 53 : pari_sp av = avma;
2264 : long r;
2265 : GEN d;
2266 53 : init_pivot_list(x);
2267 53 : d = F2m_gauss_pivot(F2m_copy(x),&r);
2268 53 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2269 : }
2270 :
2271 : GEN
2272 7 : F2xqM_indexrank(GEN x, GEN T) {
2273 7 : pari_sp av = avma;
2274 : long r;
2275 : GEN d;
2276 7 : init_pivot_list(x);
2277 7 : d = F2xqM_gauss_pivot(x, T, &r);
2278 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2279 : }
2280 :
2281 : GEN
2282 7 : FlxqM_indexrank(GEN x, GEN T, ulong p) {
2283 7 : pari_sp av = avma;
2284 : long r;
2285 : GEN d;
2286 7 : init_pivot_list(x);
2287 7 : d = FlxqM_gauss_pivot(x, T, p, &r);
2288 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2289 : }
2290 :
2291 : GEN
2292 7 : FqM_indexrank(GEN x, GEN T, GEN p) {
2293 7 : pari_sp av = avma;
2294 : long r;
2295 : GEN d;
2296 7 : init_pivot_list(x);
2297 7 : d = FqM_gauss_pivot(x, T, p, &r);
2298 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2299 : }
2300 :
2301 : /*******************************************************************/
2302 : /* */
2303 : /* Solve A*X=B (Gauss pivot) */
2304 : /* */
2305 : /*******************************************************************/
2306 : /* x a column, x0 same column in the original input matrix (for reference),
2307 : * c list of pivots so far */
2308 : static long
2309 2593082 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
2310 : {
2311 2593082 : GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
2312 2593082 : long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
2313 2593082 : if (c)
2314 : {
2315 585456 : for (i=1; i<lx; i++)
2316 365768 : if (!c[i])
2317 : {
2318 149179 : long e = gexpo(gel(x,i));
2319 149179 : if (e > ex) { ex = e; k = i; }
2320 : }
2321 : }
2322 : else
2323 : {
2324 8419163 : for (i=ix; i<lx; i++)
2325 : {
2326 6045749 : long e = gexpo(gel(x,i));
2327 6045769 : if (e > ex) { ex = e; k = i; }
2328 : }
2329 : }
2330 2593102 : if (!k) return lx;
2331 2477941 : p = gel(x,k);
2332 2477941 : r = gel(x0,k); if (isrationalzero(r)) r = x0;
2333 2477946 : return cx_approx0(p, r)? lx: k;
2334 : }
2335 : static long
2336 201820 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
2337 : {
2338 201820 : GEN x = gel(X, ix);
2339 201820 : long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
2340 201820 : if (c)
2341 : {
2342 504 : for (i=1; i<lx; i++)
2343 378 : if (!c[i] && !gequal0(gel(x,i)))
2344 : {
2345 245 : long e = gvaluation(gel(x,i), p);
2346 245 : if (e < ex) { ex = e; k = i; }
2347 : }
2348 : }
2349 : else
2350 : {
2351 1721352 : for (i=ix; i<lx; i++)
2352 1519658 : if (!gequal0(gel(x,i)))
2353 : {
2354 1147068 : long e = gvaluation(gel(x,i), p);
2355 1147068 : if (e < ex) { ex = e; k = i; }
2356 : }
2357 : }
2358 201820 : return k? k: lx;
2359 : }
2360 : static long
2361 3752 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
2362 : {
2363 3752 : GEN x = gel(X, ix);
2364 3752 : long i, lx = lg(x);
2365 : (void)x0;
2366 3752 : if (c)
2367 : {
2368 9891 : for (i=1; i<lx; i++)
2369 9002 : if (!c[i] && !gequal0(gel(x,i))) return i;
2370 : }
2371 : else
2372 : {
2373 2002 : for (i=ix; i<lx; i++)
2374 1988 : if (!gequal0(gel(x,i))) return i;
2375 : }
2376 903 : return lx;
2377 : }
2378 :
2379 : /* Set pivot seeking function appropriate for the domain of x with RgM_type t
2380 : * (first non zero pivot, maximal pivot...)
2381 : * x0 is a reference point used when guessing whether x[i,j] ~ 0
2382 : * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
2383 : * but use original x when deciding whether a prospective pivot is nonzero */
2384 : static void
2385 1908117 : set_pivot_fun(pivot_fun *fun, GEN *data, long t, GEN x0, GEN p)
2386 : {
2387 1908117 : switch(t)
2388 : {
2389 1799251 : case t_REAL:
2390 1799251 : case t_COMPLEX: *data = x0; *fun = gauss_get_pivot_max; break;
2391 26998 : case t_PADIC: *data = p; *fun = gauss_get_pivot_padic; break;
2392 81868 : default: *data = NULL; *fun = gauss_get_pivot_NZ;
2393 : }
2394 1908117 : }
2395 : static void
2396 26820 : set_pivot_fun_all(pivot_fun *fun, GEN *data, GEN x)
2397 : {
2398 : GEN p, pol;
2399 26820 : long pa, t = RgM_type(x, &p,&pol,&pa);
2400 26820 : set_pivot_fun(fun, data, t, x, p);
2401 26820 : }
2402 :
2403 : static GEN
2404 1265736 : get_col(GEN a, GEN b, GEN p, long li)
2405 : {
2406 1265736 : GEN u = cgetg(li+1,t_COL);
2407 : long i, j;
2408 :
2409 1265739 : gel(u,li) = gdiv(gel(b,li), p);
2410 5151741 : for (i=li-1; i>0; i--)
2411 : {
2412 3886014 : pari_sp av = avma;
2413 3886014 : GEN m = gel(b,i);
2414 17097574 : for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
2415 3885923 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
2416 : }
2417 1265727 : return u;
2418 : }
2419 :
2420 : /* bk -= m * bi */
2421 : static void
2422 18302316 : _submul(GEN b, long k, long i, GEN m)
2423 : {
2424 18302316 : gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
2425 18302173 : }
2426 : static int
2427 2376810 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
2428 : {
2429 2376810 : *iscol = *b ? (typ(*b) == t_COL): 0;
2430 2376810 : *aco = lg(a) - 1;
2431 2376810 : if (!*aco) /* a empty */
2432 : {
2433 70 : if (*b && lg(*b) != 1) pari_err_DIM("gauss");
2434 70 : *li = 0; return 0;
2435 : }
2436 2376740 : *li = nbrows(a);
2437 2376739 : if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
2438 2376739 : if (*b)
2439 : {
2440 2112427 : switch(typ(*b))
2441 : {
2442 121627 : case t_MAT:
2443 121627 : if (lg(*b) == 1) return 0;
2444 121627 : *b = RgM_shallowcopy(*b);
2445 121629 : break;
2446 1990801 : case t_COL:
2447 1990801 : *b = mkmat( leafcopy(*b) );
2448 1990799 : break;
2449 0 : default: pari_err_TYPE("gauss",*b);
2450 : }
2451 2112428 : if (nbrows(*b) != *li) pari_err_DIM("gauss");
2452 : }
2453 : else
2454 264312 : *b = matid(*li);
2455 2376740 : return 1;
2456 : }
2457 :
2458 : static GEN
2459 2051 : RgM_inv_FpM(GEN a, GEN p)
2460 : {
2461 : ulong pp;
2462 2051 : a = RgM_Fp_init(a, p, &pp);
2463 2051 : switch(pp)
2464 : {
2465 35 : case 0:
2466 35 : a = FpM_inv(a,p);
2467 35 : if (a) a = FpM_to_mod(a, p);
2468 35 : break;
2469 189 : case 2:
2470 189 : a = F2m_inv(a);
2471 189 : if (a) a = F2m_to_mod(a);
2472 189 : break;
2473 1827 : default:
2474 1827 : a = Flm_inv_sp(a, NULL, pp);
2475 1827 : if (a) a = Flm_to_mod(a, pp);
2476 : }
2477 2051 : return a;
2478 : }
2479 :
2480 : static GEN
2481 42 : RgM_inv_FqM(GEN x, GEN pol, GEN p)
2482 : {
2483 42 : pari_sp av = avma;
2484 42 : GEN b, T = RgX_to_FpX(pol, p);
2485 42 : if (signe(T) == 0) pari_err_OP("^",x,gen_m1);
2486 42 : b = FqM_inv(RgM_to_FqM(x, T, p), T, p);
2487 42 : if (!b) return gc_NULL(av);
2488 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
2489 : }
2490 :
2491 : /* Returns gen_0 instead of NULL for 'no fast algorithm'. NULL is already
2492 : * reserved for 'not invertible' */
2493 : static GEN
2494 526325 : RgM_inv_fast(GEN x, pivot_fun *fun, GEN *data)
2495 : {
2496 : GEN p, pol;
2497 526325 : long pa, t = RgM_type(x, &p,&pol,&pa);
2498 526323 : set_pivot_fun(fun, data, t, x, p);
2499 526323 : switch(t)
2500 : {
2501 48496 : case t_INT: /* Fall back */
2502 48496 : case t_FRAC: return QM_inv(x);
2503 147 : case t_FFELT: return FFM_inv(x, pol);
2504 2051 : case t_INTMOD: return RgM_inv_FpM(x, p);
2505 42 : case RgX_type_code(t_POLMOD, t_INTMOD):
2506 42 : return RgM_inv_FqM(x, pol, p);
2507 475587 : default: return gen_0;
2508 : }
2509 : }
2510 :
2511 : static GEN
2512 63 : RgM_RgC_solve_FpC(GEN a, GEN b, GEN p)
2513 : {
2514 63 : pari_sp av = avma;
2515 : ulong pp;
2516 63 : a = RgM_Fp_init(a, p, &pp);
2517 63 : switch(pp)
2518 : {
2519 14 : case 0:
2520 14 : b = RgC_to_FpC(b, p);
2521 14 : a = FpM_FpC_gauss(a,b,p);
2522 14 : return a ? gerepileupto(av, FpC_to_mod(a, p)): NULL;
2523 28 : case 2:
2524 28 : b = RgV_to_F2v(b);
2525 28 : a = F2m_F2c_gauss(a,b);
2526 28 : return a ? gerepileupto(av, F2c_to_mod(a)): NULL;
2527 21 : default:
2528 21 : b = RgV_to_Flv(b, pp);
2529 21 : a = Flm_Flc_gauss(a, b, pp);
2530 21 : return a ? gerepileupto(av, Flc_to_mod(a, pp)): NULL;
2531 : }
2532 : }
2533 :
2534 : static GEN
2535 105 : RgM_solve_FpM(GEN a, GEN b, GEN p)
2536 : {
2537 105 : pari_sp av = avma;
2538 : ulong pp;
2539 105 : a = RgM_Fp_init(a, p, &pp);
2540 105 : switch(pp)
2541 : {
2542 35 : case 0:
2543 35 : b = RgM_to_FpM(b, p);
2544 35 : a = FpM_gauss(a,b,p);
2545 35 : return a ? gerepileupto(av, FpM_to_mod(a, p)): NULL;
2546 28 : case 2:
2547 28 : b = RgM_to_F2m(b);
2548 28 : a = F2m_gauss(a,b);
2549 28 : return a ? gerepileupto(av, F2m_to_mod(a)): NULL;
2550 42 : default:
2551 42 : b = RgM_to_Flm(b, pp);
2552 42 : a = Flm_gauss(a,b,pp);
2553 42 : return a ? gerepileupto(av, Flm_to_mod(a, pp)): NULL;
2554 : }
2555 : }
2556 :
2557 : /* Gaussan Elimination. If a is square, return a^(-1)*b;
2558 : * if a has more rows than columns and b is NULL, return c such that c a = Id.
2559 : * a is a (not necessarily square) matrix
2560 : * b is a matrix or column vector, NULL meaning: take the identity matrix,
2561 : * effectively returning the inverse of a
2562 : * If a and b are empty, the result is the empty matrix.
2563 : *
2564 : * li: number of rows of a and b
2565 : * aco: number of columns of a
2566 : * bco: number of columns of b (if matrix)
2567 : */
2568 : static GEN
2569 1692208 : RgM_solve_basecase(GEN a, GEN b, pivot_fun pivot, GEN data)
2570 : {
2571 1692208 : pari_sp av = avma;
2572 : long i, j, k, li, bco, aco;
2573 : int iscol;
2574 : GEN p, u;
2575 :
2576 1692208 : if (lg(a)-1 == 2 && nbrows(a) == 2)
2577 : { /* 2x2 matrix, start by inverting a */
2578 1026229 : GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
2579 1026229 : GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
2580 1026229 : GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
2581 1026231 : if (gequal0(D)) return NULL;
2582 1026233 : ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
2583 1026232 : ainv = RgM_Rg_mul(ainv, ginv(D));
2584 1026219 : if (b) ainv = gmul(ainv, b);
2585 1026224 : return gerepileupto(av, ainv);
2586 : }
2587 665978 : if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
2588 665985 : a = RgM_shallowcopy(a);
2589 665985 : bco = lg(b)-1;
2590 665985 : if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
2591 :
2592 665985 : p = NULL; /* gcc -Wall */
2593 2291861 : for (i=1; i<=aco; i++)
2594 : {
2595 : /* k is the line where we find the pivot */
2596 2291852 : k = pivot(a, data, i, NULL);
2597 2291874 : if (k > li) return NULL;
2598 2291859 : if (k != i)
2599 : { /* exchange the lines s.t. k = i */
2600 1795137 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
2601 1739214 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
2602 : }
2603 2291859 : p = gcoeff(a,i,i);
2604 2291859 : if (i == aco) break;
2605 :
2606 5118323 : for (k=i+1; k<=li; k++)
2607 : {
2608 3492458 : GEN m = gcoeff(a,k,i);
2609 3492458 : if (!gequal0(m))
2610 : {
2611 2837488 : m = gdiv(m,p);
2612 12115778 : for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
2613 11861817 : for (j=1; j<=bco; j++) _submul(gel(b,j),k,i,m);
2614 : }
2615 : }
2616 1625865 : if (gc_needed(av,1))
2617 : {
2618 12 : if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
2619 12 : gerepileall(av,2, &a,&b);
2620 : }
2621 : }
2622 :
2623 665977 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
2624 665977 : u = cgetg(bco+1,t_MAT);
2625 1931695 : for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
2626 665959 : return gc_GEN(av, iscol? gel(u,1): u);
2627 : }
2628 :
2629 : /* Returns gen_0 instead of NULL for 'no fast algorithm'. NULL is already
2630 : * reserved for 'not invertible' */
2631 : static GEN
2632 1177449 : RgM_RgC_solve_fast(GEN x, GEN y, pivot_fun *fun, GEN *data)
2633 : {
2634 : GEN p, pol;
2635 1177449 : long pa, t = RgM_RgC_type(x, y, &p,&pol,&pa);
2636 1177445 : set_pivot_fun(fun, data, t, x, p);
2637 1177446 : switch(t)
2638 : {
2639 9238 : case t_INT: return ZM_gauss(x, y);
2640 7 : case t_FRAC: return QM_gauss(x, y);
2641 63 : case t_INTMOD: return RgM_RgC_solve_FpC(x, y, p);
2642 42 : case t_FFELT: return FFM_FFC_gauss(x, y, pol);
2643 1168096 : default: return gen_0;
2644 : }
2645 : }
2646 : static GEN
2647 48782 : RgM_solve_fast(GEN x, GEN y, pivot_fun *fun, GEN *data)
2648 : {
2649 : GEN p, pol;
2650 48782 : long pa, t = RgM_type2(x, y, &p,&pol,&pa);
2651 48780 : set_pivot_fun(fun, data, t, x, p);
2652 48780 : switch(t)
2653 : {
2654 77 : case t_INT: return ZM_gauss(x, y);
2655 14 : case t_FRAC: return QM_gauss(x, y);
2656 105 : case t_INTMOD: return RgM_solve_FpM(x, y, p);
2657 56 : case t_FFELT: return FFM_gauss(x, y, pol);
2658 48528 : default: return gen_0;
2659 : }
2660 : }
2661 :
2662 : GEN
2663 1226231 : RgM_solve(GEN a, GEN b)
2664 : {
2665 1226231 : pari_sp av = avma;
2666 : pivot_fun fun;
2667 : GEN u, data;
2668 1226231 : if (!b) return RgM_inv(a);
2669 48782 : u = typ(b)==t_MAT ? RgM_solve_fast(a, b, &fun, &data)
2670 1226231 : : RgM_RgC_solve_fast(a, b, &fun, &data);
2671 1226226 : if (!u) return gc_NULL(av);
2672 1226121 : if (u != gen_0) return u;
2673 1216624 : return RgM_solve_basecase(a, b, fun, data);
2674 : }
2675 :
2676 : GEN
2677 28 : RgM_div(GEN a, GEN b)
2678 : {
2679 28 : pari_sp av = avma;
2680 28 : GEN u = RgM_solve(shallowtrans(b), shallowtrans(a));
2681 28 : if (!u) return gc_NULL(av);
2682 21 : return gc_GEN(av, shallowtrans(u));
2683 : }
2684 :
2685 : GEN
2686 526325 : RgM_inv(GEN a)
2687 : {
2688 : pivot_fun fun;
2689 526325 : GEN data, b = RgM_inv_fast(a, &fun, &data);
2690 526308 : return b==gen_0? RgM_solve_basecase(a, NULL, fun, data): b;
2691 : }
2692 :
2693 : /* assume dim A >= 1, A invertible + upper triangular */
2694 : static GEN
2695 3230732 : RgM_inv_upper_ind(GEN A, long index)
2696 : {
2697 3230732 : long n = lg(A)-1, i = index, j;
2698 3230732 : GEN u = zerocol(n);
2699 3230734 : gel(u,i) = ginv(gcoeff(A,i,i));
2700 6535339 : for (i--; i>0; i--)
2701 : {
2702 3304614 : pari_sp av = avma;
2703 3304614 : GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
2704 14642789 : for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
2705 3304595 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
2706 : }
2707 3230725 : return u;
2708 : }
2709 : GEN
2710 1617388 : RgM_inv_upper(GEN A)
2711 : {
2712 : long i, l;
2713 1617388 : GEN B = cgetg_copy(A, &l);
2714 4848108 : for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
2715 1617380 : return B;
2716 : }
2717 :
2718 : static GEN
2719 4516846 : split_realimag_col(GEN z, long r1, long r2)
2720 : {
2721 4516846 : long i, ru = r1+r2;
2722 4516846 : GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
2723 12537845 : for (i=1; i<=r1; i++) {
2724 8021001 : GEN a = gel(z,i);
2725 8021001 : if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
2726 8021001 : gel(x,i) = a;
2727 : }
2728 7223827 : for ( ; i<=ru; i++) {
2729 2706983 : GEN b, a = gel(z,i);
2730 2706983 : if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
2731 2706983 : gel(x,i) = a;
2732 2706983 : gel(y,i) = b;
2733 : }
2734 4516844 : return x;
2735 : }
2736 : GEN
2737 2570062 : split_realimag(GEN x, long r1, long r2)
2738 : {
2739 2570062 : if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
2740 4502686 : pari_APPLY_same(split_realimag_col(gel(x,i), r1, r2));
2741 : }
2742 :
2743 : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
2744 : * r1 first lines of M,y are real. Solve the system obtained by splitting
2745 : * real and imaginary parts. */
2746 : GEN
2747 1215494 : RgM_solve_realimag(GEN M, GEN y)
2748 : {
2749 1215494 : long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
2750 1215494 : return RgM_solve(split_realimag(M, r1,r2),
2751 : split_realimag(y, r1,r2));
2752 : }
2753 :
2754 : GEN
2755 434 : gauss(GEN a, GEN b)
2756 : {
2757 : GEN z;
2758 434 : long t = typ(b);
2759 434 : if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
2760 434 : if (t!=t_COL && t!=t_MAT) pari_err_TYPE("gauss",b);
2761 434 : z = RgM_solve(a,b);
2762 434 : if (!z) pari_err_INV("gauss",a);
2763 329 : return z;
2764 : }
2765 :
2766 : /* #C = n, C[z[i]] = K[i], complete by 0s */
2767 : static GEN
2768 14 : RgC_inflate(GEN K, GEN z, long n)
2769 : {
2770 14 : GEN c = zerocol(n);
2771 14 : long j, l = lg(K);
2772 28 : for (j = 1; j < l; j++) gel(c, z[j]) = gel(K, j);
2773 14 : return c;
2774 : }
2775 : /* in place: C[i] *= cB / v[i] */
2776 : static void
2777 6314 : QC_normalize(GEN C, GEN v, GEN cB)
2778 : {
2779 6314 : long l = lg(C), i;
2780 47838 : for (i = 1; i < l; i++)
2781 : {
2782 41524 : GEN c = cB, k = gel(C,i), d = gel(v,i);
2783 41524 : if (d)
2784 : {
2785 24619 : if (isintzero(d)) { gel(C,i) = gen_0; continue; }
2786 24619 : c = div_content(c, d);
2787 : }
2788 41524 : gel(C,i) = c? gmul(k,c): k;
2789 : }
2790 6314 : }
2791 :
2792 : /* same as above, M rational; if flag = 1, call indexrank and return 1 sol */
2793 : GEN
2794 6307 : QM_gauss_i(GEN M, GEN B, long flag)
2795 : {
2796 6307 : pari_sp av = avma;
2797 : long i, l, n;
2798 6307 : int col = typ(B) == t_COL;
2799 6307 : GEN K, cB, N = cgetg_copy(M, &l), v = cgetg(l, t_VEC), z2 = NULL;
2800 :
2801 47859 : for (i = 1; i < l; i++)
2802 41552 : gel(N,i) = Q_primitive_part(gel(M,i), &gel(v,i));
2803 6307 : if (flag)
2804 : {
2805 329 : GEN z = ZM_indexrank(N), z1 = gel(z,1);
2806 329 : z2 = gel(z,2);
2807 329 : N = shallowmatextract(N, z1, z2);
2808 329 : B = col? vecpermute(B,z1): rowpermute(B,z1);
2809 329 : if (lg(z2) == l) z2 = NULL; else v = vecpermute(v, z2);
2810 : }
2811 6307 : B = Q_primitive_part(B, &cB);
2812 6307 : K = ZM_gauss(N, B); if (!K) return gc_NULL(av);
2813 6307 : n = l - 1;
2814 6307 : if (col)
2815 : {
2816 6279 : QC_normalize(K, v, cB);
2817 6279 : if (z2) K = RgC_inflate(K, z2, n);
2818 : }
2819 : else
2820 : {
2821 28 : long lK = lg(K);
2822 63 : for (i = 1; i < lK; i++)
2823 : {
2824 35 : QC_normalize(gel(K,i), v, cB);
2825 35 : if (z2) gel(K,i) = RgC_inflate(gel(K,i), z2, n);
2826 : }
2827 : }
2828 6307 : return gc_GEN(av, K);
2829 : }
2830 : GEN
2831 5978 : QM_gauss(GEN M, GEN B) { return QM_gauss_i(M, B, 0); }
2832 :
2833 : static GEN
2834 794267 : ZM_inv_slice(GEN A, GEN P, GEN *mod)
2835 : {
2836 794267 : pari_sp av = avma;
2837 794267 : long i, n = lg(P)-1;
2838 : GEN H, T;
2839 794267 : if (n == 1)
2840 : {
2841 761622 : ulong p = uel(P,1);
2842 761622 : GEN Hp, a = ZM_to_Flm(A, p);
2843 761623 : Hp = Flm_adjoint(a, p);
2844 761621 : Hp = gerepileupto(av, Flm_to_ZM(Hp));
2845 761623 : *mod = utoipos(p); return Hp;
2846 : }
2847 32645 : T = ZV_producttree(P);
2848 32645 : A = ZM_nv_mod_tree(A, P, T);
2849 32645 : H = cgetg(n+1, t_VEC);
2850 182006 : for(i=1; i <= n; i++)
2851 149361 : gel(H,i) = Flm_adjoint(gel(A, i), uel(P,i));
2852 32645 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
2853 32645 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2854 : }
2855 :
2856 : static GEN
2857 720846 : RgM_true_Hadamard(GEN a)
2858 : {
2859 720846 : pari_sp av = avma;
2860 720846 : long n = lg(a)-1, i;
2861 : GEN B;
2862 720846 : if (n == 0) return gen_1;
2863 720846 : a = RgM_gtofp(a, LOWDEFAULTPREC);
2864 720848 : B = gnorml2(gel(a,1));
2865 2988630 : for (i = 2; i <= n; i++) B = gmul(B, gnorml2(gel(a,i)));
2866 720843 : return gc_INT(av, ceil_safe(sqrtr(B)));
2867 : }
2868 :
2869 : GEN
2870 794267 : ZM_inv_worker(GEN P, GEN A)
2871 : {
2872 794267 : GEN V = cgetg(3, t_VEC);
2873 794267 : gel(V,1) = ZM_inv_slice(A, P, &gel(V,2));
2874 794268 : return V;
2875 : }
2876 :
2877 : static GEN
2878 43533 : ZM_inv0(GEN A, GEN *pden)
2879 : {
2880 43533 : if (pden) *pden = gen_1;
2881 43533 : (void)A; return cgetg(1, t_MAT);
2882 : }
2883 : static GEN
2884 644319 : ZM_inv1(GEN A, GEN *pden)
2885 : {
2886 644319 : GEN a = gcoeff(A,1,1);
2887 644319 : long s = signe(a);
2888 644319 : if (!s) return NULL;
2889 644319 : if (pden) *pden = absi(a);
2890 644319 : retmkmat(mkcol(s == 1? gen_1: gen_m1));
2891 : }
2892 : static GEN
2893 725225 : ZM_inv2(GEN A, GEN *pden)
2894 : {
2895 : GEN a, b, c, d, D, cA;
2896 : long s;
2897 725225 : A = Q_primitive_part(A, &cA);
2898 725220 : a = gcoeff(A,1,1); b = gcoeff(A,1,2);
2899 725220 : c = gcoeff(A,2,1); d = gcoeff(A,2,2);
2900 725220 : D = subii(mulii(a,d), mulii(b,c)); /* left on stack */
2901 725216 : s = signe(D);
2902 725216 : if (!s) return NULL;
2903 725202 : if (s < 0) D = negi(D);
2904 725206 : if (pden) *pden = mul_denom(D, cA);
2905 725203 : if (s > 0)
2906 683577 : retmkmat2(mkcol2(icopy(d), negi(c)), mkcol2(negi(b), icopy(a)));
2907 : else
2908 41626 : retmkmat2(mkcol2(negi(d), icopy(c)), mkcol2(icopy(b), negi(a)));
2909 : }
2910 :
2911 : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
2912 : * not available. Return H primitive such that M*H = den*Id */
2913 : GEN
2914 0 : ZM_inv_ratlift(GEN M, GEN *pden)
2915 : {
2916 0 : pari_sp av2, av = avma;
2917 : GEN Hp, q, H;
2918 : ulong p;
2919 0 : long m = lg(M)-1;
2920 : forprime_t S;
2921 : pari_timer ti;
2922 :
2923 0 : if (m == 0) return ZM_inv0(M,pden);
2924 0 : if (m == 1 && nbrows(M)==1) return ZM_inv1(M,pden);
2925 0 : if (m == 2 && nbrows(M)==2) return ZM_inv2(M,pden);
2926 :
2927 0 : if (DEBUGLEVEL>5) timer_start(&ti);
2928 0 : init_modular_big(&S);
2929 0 : av2 = avma;
2930 0 : H = NULL;
2931 0 : while ((p = u_forprime_next(&S)))
2932 : {
2933 : GEN Mp, B, Hr;
2934 0 : Mp = ZM_to_Flm(M,p);
2935 0 : Hp = Flm_inv_sp(Mp, NULL, p);
2936 0 : if (!Hp) continue;
2937 0 : if (!H)
2938 : {
2939 0 : H = ZM_init_CRT(Hp, p);
2940 0 : q = utoipos(p);
2941 : }
2942 : else
2943 0 : ZM_incremental_CRT(&H, Hp, &q, p);
2944 0 : B = sqrti(shifti(q,-1));
2945 0 : Hr = FpM_ratlift(H,q,B,B,NULL);
2946 0 : if (DEBUGLEVEL>5)
2947 0 : timer_printf(&ti,"ZM_inv mod %lu (ratlift=%ld)", p,!!Hr);
2948 0 : if (Hr) {/* DONE ? */
2949 0 : GEN Hl = Q_remove_denom(Hr, pden);
2950 0 : if (ZM_isscalar(ZM_mul(Hl, M), *pden)) { H = Hl; break; }
2951 : }
2952 :
2953 0 : if (gc_needed(av,2))
2954 : {
2955 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
2956 0 : gerepileall(av2, 2, &H, &q);
2957 : }
2958 : }
2959 0 : if (!*pden) *pden = gen_1;
2960 0 : return gc_all(av, 2, &H, pden);
2961 : }
2962 :
2963 : GEN
2964 76806 : FpM_ratlift_worker(GEN A, GEN mod, GEN B)
2965 : {
2966 : long l, i;
2967 76806 : GEN H = cgetg_copy(A, &l);
2968 161616 : for (i = 1; i < l; i++)
2969 : {
2970 84814 : GEN c = FpC_ratlift(gel(A,i), mod, B, B, NULL);
2971 84811 : gel(H,i) = c? c: gen_0;
2972 : }
2973 76802 : return H;
2974 : }
2975 : static int
2976 765984 : can_ratlift(GEN x, GEN mod, GEN B)
2977 : {
2978 765984 : pari_sp av = avma;
2979 : GEN a, b;
2980 765984 : return gc_bool(av, Fp_ratlift(x, mod, B, B, &a,&b));
2981 : }
2982 : static GEN
2983 2738333 : FpM_ratlift_parallel(GEN A, GEN mod, GEN B)
2984 : {
2985 2738333 : pari_sp av = avma;
2986 : GEN worker;
2987 2738333 : long i, l = lg(A), m = mt_nbthreads();
2988 2738333 : int test = !!B;
2989 :
2990 2738333 : if (l == 1 || lgcols(A) == 1) return gcopy(A);
2991 2738334 : if (!B) B = sqrti(shifti(mod,-1));
2992 2738281 : if (m == 1 || l == 2 || lgcols(A) < 10)
2993 : {
2994 2730498 : A = FpM_ratlift(A, mod, B, B, NULL);
2995 2730549 : return A? A: gc_NULL(av);
2996 : }
2997 : /* test one coefficient first */
2998 7783 : if (test && !can_ratlift(gcoeff(A,1,1), mod, B)) return gc_NULL(av);
2999 7665 : worker = snm_closure(is_entry("_FpM_ratlift_worker"), mkvec2(mod,B));
3000 7665 : A = gen_parapply_slice(worker, A, m);
3001 84030 : for (i = 1; i < l; i++) if (typ(gel(A,i)) != t_COL) return gc_NULL(av);
3002 6657 : return A;
3003 : }
3004 :
3005 : static GEN
3006 758722 : ZM_adj_ratlift(GEN A, GEN H, GEN mod, GEN T)
3007 : {
3008 758722 : pari_sp av = avma;
3009 : GEN B, D, g;
3010 758722 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3011 758722 : if (T) D = mulii(T, D);
3012 758722 : g = gcdii(D, mod);
3013 758722 : if (!equali1(g))
3014 : {
3015 14 : mod = diviiexact(mod, g);
3016 14 : H = FpM_red(H, mod);
3017 : }
3018 758722 : D = Fp_inv(Fp_red(D, mod), mod);
3019 : /* test 1 coeff first */
3020 758722 : B = sqrti(shifti(mod,-1));
3021 758718 : if (!can_ratlift(Fp_mul(D, gcoeff(A,1,1), mod), mod, B)) return gc_NULL(av);
3022 737553 : H = FpM_Fp_mul(H, D, mod);
3023 737549 : H = FpM_ratlift_parallel(H, mod, B);
3024 737554 : return H? H: gc_NULL(av);
3025 : }
3026 :
3027 : /* if (T) return T A^(-1) in Mn(Q), else B in Mn(Z) such that A B = den*Id */
3028 : static GEN
3029 2133929 : ZM_inv_i(GEN A, GEN *pden, GEN T)
3030 : {
3031 2133929 : pari_sp av = avma;
3032 2133929 : long m = lg(A)-1, n, k1 = 1, k2;
3033 2133929 : GEN H = NULL, D, H1 = NULL, mod1 = NULL, worker;
3034 : ulong bnd, mask;
3035 : forprime_t S;
3036 : pari_timer ti;
3037 :
3038 2133929 : if (m == 0) return ZM_inv0(A,pden);
3039 2090396 : if (pden) *pden = gen_1;
3040 2090396 : if (nbrows(A) < m) return NULL;
3041 2090388 : if (m == 1 && nbrows(A)==1 && !T) return ZM_inv1(A,pden);
3042 1446070 : if (m == 2 && nbrows(A)==2 && !T) return ZM_inv2(A,pden);
3043 :
3044 720845 : if (DEBUGLEVEL>=5) timer_start(&ti);
3045 720845 : init_modular_big(&S);
3046 720846 : bnd = expi(RgM_true_Hadamard(A));
3047 720846 : worker = snm_closure(is_entry("_ZM_inv_worker"), mkvec(A));
3048 720849 : gen_inccrt("ZM_inv_r", worker, NULL, k1, 0, &S, &H1, &mod1, nmV_chinese_center, FpM_center);
3049 720848 : n = (bnd+1)/expu(S.p)+1;
3050 720848 : if (DEBUGLEVEL>=5) timer_printf(&ti,"inv (%ld/%ld primes)", k1, n);
3051 720848 : mask = quadratic_prec_mask(n);
3052 720848 : for (k2 = 0;;)
3053 66171 : {
3054 : GEN Hr;
3055 787019 : if (k2 > 0)
3056 : {
3057 58845 : gen_inccrt("ZM_inv_r", worker, NULL, k2, 0, &S, &H1, &mod1,nmV_chinese_center,FpM_center);
3058 58845 : k1 += k2;
3059 58845 : if (DEBUGLEVEL>=5) timer_printf(&ti,"CRT (%ld/%ld primes)", k1, n);
3060 : }
3061 787019 : if (mask == 1) break;
3062 758722 : k2 = (mask&1UL) ? k1-1: k1;
3063 758722 : mask >>= 1;
3064 :
3065 758722 : Hr = ZM_adj_ratlift(A, H1, mod1, T);
3066 758722 : if (DEBUGLEVEL>=5) timer_printf(&ti,"ratlift (%ld/%ld primes)", k1, n);
3067 758722 : if (Hr) {/* DONE ? */
3068 696426 : GEN Hl = Q_primpart(Hr), R = ZM_mul(Hl, A), d = gcoeff(R,1,1);
3069 696425 : if (gsigne(d) < 0) { d = gneg(d); Hl = ZM_neg(Hl); }
3070 696425 : if (DEBUGLEVEL>=5) timer_printf(&ti,"mult (%ld/%ld primes)", k1, n);
3071 696425 : if (equali1(d))
3072 : {
3073 595840 : if (ZM_isidentity(R)) { H = Hl; break; }
3074 : }
3075 100585 : else if (ZM_isscalar(R, d))
3076 : {
3077 96712 : if (T) T = gdiv(T,d);
3078 89635 : else if (pden) *pden = d;
3079 96712 : H = Hl; break;
3080 : }
3081 : }
3082 : }
3083 720848 : if (!H)
3084 : {
3085 : GEN d;
3086 28297 : H = H1;
3087 28297 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3088 28297 : if (signe(D)==0) pari_err_INV("ZM_inv", A);
3089 28297 : if (T) T = gdiv(T, D);
3090 : else
3091 : {
3092 27148 : d = gcdii(Q_content_safe(H), D);
3093 27148 : if (signe(D) < 0) d = negi(d);
3094 27148 : if (!equali1(d))
3095 : {
3096 15426 : H = ZM_Z_divexact(H, d);
3097 15426 : D = diviiexact(D, d);
3098 : }
3099 27149 : if (pden) *pden = D;
3100 : }
3101 : }
3102 720849 : if (T && !isint1(T)) H = ZM_Q_mul(H, T);
3103 720849 : return gc_all(av, pden? 2: 1, &H, pden);
3104 : }
3105 : GEN
3106 2068604 : ZM_inv(GEN A, GEN *pden) { return ZM_inv_i(A, pden, NULL); }
3107 :
3108 : /* same as above, M rational */
3109 : GEN
3110 65325 : QM_inv(GEN M)
3111 : {
3112 65325 : pari_sp av = avma;
3113 : GEN den, dM, K;
3114 65325 : M = Q_remove_denom(M, &dM);
3115 65325 : K = ZM_inv_i(M, &den, dM);
3116 65325 : if (!K) return gc_NULL(av);
3117 65304 : if (den && !equali1(den)) K = ZM_Q_mul(K, ginv(den));
3118 65290 : return gerepileupto(av, K);
3119 : }
3120 :
3121 : static GEN
3122 105428 : ZM_ker_filter(GEN A, GEN P)
3123 : {
3124 105428 : long i, j, l = lg(A), n = 1, d = lg(gmael(A,1,1));
3125 105428 : GEN B, Q, D = gmael(A,1,2);
3126 215597 : for (i=2; i<l; i++)
3127 : {
3128 110169 : GEN Di = gmael(A,i,2);
3129 110169 : long di = lg(gmael(A,i,1));
3130 110169 : int c = vecsmall_lexcmp(D, Di);
3131 110169 : if (di==d && c==0) n++;
3132 45588 : else if (d > di || (di==d && c>0))
3133 37680 : { n = 1; d = di; D = Di; }
3134 : }
3135 105428 : B = cgetg(n+1, t_VEC);
3136 105428 : Q = cgetg(n+1, typ(P));
3137 321025 : for (i=1, j=1; i<l; i++)
3138 : {
3139 215597 : if (lg(gmael(A,i,1))==d && vecsmall_lexcmp(D, gmael(A,i,2))==0)
3140 : {
3141 170009 : gel(B,j) = gmael(A,i,1);
3142 170009 : Q[j] = P[i];
3143 170009 : j++;
3144 : }
3145 : }
3146 105428 : return mkvec3(B,Q,D);
3147 : }
3148 :
3149 : static GEN
3150 69755 : ZM_ker_chinese(GEN A, GEN P, GEN *mod)
3151 : {
3152 69755 : GEN BQD = ZM_ker_filter(A, P);
3153 69755 : return mkvec2(nmV_chinese_center(gel(BQD,1), gel(BQD,2), mod), gel(BQD,3));
3154 : }
3155 :
3156 : static GEN
3157 133560 : ZM_ker_slice(GEN A, GEN P, GEN *mod)
3158 : {
3159 133560 : pari_sp av = avma;
3160 133560 : long i, n = lg(P)-1;
3161 : GEN BQD, B, Q, D, H, HD, T;
3162 133560 : if (n == 1)
3163 : {
3164 97887 : ulong p = uel(P,1);
3165 97887 : GEN K = Flm_ker_sp(ZM_to_Flm(A, p), p, 2);
3166 97887 : *mod = utoipos(p); return mkvec2(Flm_to_ZM(gel(K,1)), gel(K,2));
3167 : }
3168 35673 : T = ZV_producttree(P);
3169 35673 : A = ZM_nv_mod_tree(A, P, T);
3170 35673 : H = cgetg(n+1, t_VEC);
3171 111525 : for(i=1 ; i <= n; i++)
3172 75852 : gel(H,i) = Flm_ker_sp(gel(A, i), P[i], 2);
3173 35673 : BQD = ZM_ker_filter(H, P);
3174 35673 : B = gel(BQD,1); Q = gel(BQD,2); D = gel(BQD, 3);
3175 35673 : if (lg(Q) != lg(P)) T = ZV_producttree(Q);
3176 35673 : H = nmV_chinese_center_tree_seq(B, Q, T, ZV_chinesetree(Q,T));
3177 35673 : *mod = gmael(T, lg(T)-1, 1);
3178 35673 : HD = mkvec2(H, D);
3179 35673 : return gc_all(av, 2, &HD, mod);
3180 : }
3181 :
3182 : GEN
3183 133560 : ZM_ker_worker(GEN P, GEN A)
3184 : {
3185 133560 : GEN V = cgetg(3, t_VEC);
3186 133560 : gel(V,1) = ZM_ker_slice(A, P, &gel(V,2));
3187 133560 : return V;
3188 : }
3189 :
3190 : /* assume lg(A) > 1 */
3191 : static GEN
3192 66629 : ZM_ker_i(GEN A)
3193 : {
3194 : pari_sp av;
3195 66629 : long k, m = lg(A)-1;
3196 66629 : GEN HD = NULL, mod = gen_1, worker;
3197 : forprime_t S;
3198 :
3199 66629 : if (m >= 2*nbrows(A))
3200 : {
3201 3059 : GEN v = ZM_indexrank(A), y = gel(v,2), z = indexcompl(y, m);
3202 : GEN B, A1, A1i, d;
3203 3059 : A = rowpermute(A, gel(v,1)); /* same kernel */
3204 3059 : A1 = vecpermute(A, y); /* maximal rank submatrix */
3205 3059 : B = vecpermute(A, z);
3206 3059 : A1i = ZM_inv(A1, &d);
3207 3059 : if (!d) d = gen_1;
3208 3059 : B = vconcat(ZM_mul(ZM_neg(A1i), B), scalarmat_shallow(d, lg(B)-1));
3209 3059 : if (!gequal(y, identity_perm(lg(y)-1)))
3210 685 : B = rowpermute(B, perm_inv(shallowconcat(y,z)));
3211 3059 : return vec_Q_primpart(B);
3212 : }
3213 63570 : init_modular_big(&S);
3214 63570 : worker = snm_closure(is_entry("_ZM_ker_worker"), mkvec(A));
3215 63570 : av = avma;
3216 63570 : for (k = 1;; k <<= 1)
3217 65542 : {
3218 : pari_timer ti;
3219 : GEN H, Hr;
3220 129112 : gen_inccrt_i("ZM_ker", worker, NULL, (k+1)>>1, 0,
3221 : &S, &HD, &mod, ZM_ker_chinese, NULL);
3222 129112 : gerepileall(av, 2, &HD, &mod);
3223 146258 : H = gel(HD, 1); if (lg(H) == 1) return H;
3224 82688 : if (DEBUGLEVEL >= 4) timer_start(&ti);
3225 82688 : Hr = FpM_ratlift_parallel(H, mod, NULL);
3226 82688 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: ratlift (%ld)",!!Hr);
3227 82688 : if (Hr)
3228 : {
3229 : GEN MH;
3230 71743 : Hr = vec_Q_primpart(Hr);
3231 71743 : MH = ZM_mul(A, Hr);
3232 71743 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: QM_mul");
3233 71743 : if (ZM_equal0(MH)) return Hr;
3234 : }
3235 : }
3236 : }
3237 :
3238 : GEN
3239 49269 : ZM_ker(GEN M)
3240 : {
3241 49269 : pari_sp av = avma;
3242 49269 : long l = lg(M)-1;
3243 49269 : if (l==0) return cgetg(1, t_MAT);
3244 49269 : if (lgcols(M)==1) return matid(l);
3245 49269 : return gc_GEN(av, ZM_ker_i(M));
3246 : }
3247 :
3248 : static GEN
3249 2018601 : ZM_gauss_slice(GEN A, GEN B, GEN P, GEN *mod)
3250 : {
3251 2018601 : pari_sp av = avma;
3252 2018601 : long i, n = lg(P)-1;
3253 : GEN H, T;
3254 2018601 : if (n == 1)
3255 : {
3256 1946445 : ulong p = uel(P,1);
3257 1946445 : GEN Hp = Flm_gauss(ZM_to_Flm(A, p) , ZM_to_Flm(B, p) ,p);
3258 1946457 : if (!Hp) { *mod=gen_1; return zeromat(lg(A)-1,lg(B)-1); }
3259 1946457 : Hp = gerepileupto(av, Flm_to_ZM(Hp));
3260 1946452 : *mod = utoipos(p); return Hp;
3261 : }
3262 72156 : T = ZV_producttree(P);
3263 72156 : A = ZM_nv_mod_tree(A, P, T);
3264 72156 : B = ZM_nv_mod_tree(B, P, T);
3265 72156 : H = cgetg(n+1, t_VEC);
3266 452760 : for(i=1; i <= n; i++)
3267 : {
3268 380604 : GEN Hi = Flm_gauss(gel(A, i), gel(B,i), uel(P,i));
3269 380604 : gel(H,i) = Hi ? Hi: zero_Flm(lg(A)-1,lg(B)-1);
3270 380604 : if (!Hi) uel(P,i)=1;
3271 : }
3272 72156 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
3273 72156 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3274 : }
3275 :
3276 : GEN
3277 2018603 : ZM_gauss_worker(GEN P, GEN A, GEN B)
3278 : {
3279 2018603 : GEN V = cgetg(3, t_VEC);
3280 2018601 : gel(V,1) = ZM_gauss_slice(A, B, P, &gel(V,2));
3281 2018607 : return V;
3282 : }
3283 :
3284 : /* assume lg(A) > 1 */
3285 : static GEN
3286 1710829 : ZM_gauss_i(GEN A, GEN B)
3287 : {
3288 : pari_sp av;
3289 : long k, m, ncol;
3290 : int iscol;
3291 1710829 : GEN y, y1, y2, Hr, H = NULL, mod = gen_1, worker;
3292 : forprime_t S;
3293 1710829 : if (!init_gauss(A, &B, &m, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
3294 1710757 : init_modular_big(&S);
3295 1710759 : y = ZM_indexrank(A); y1 = gel(y,1); y2 = gel(y,2);
3296 1710772 : if (lg(y2)-1 != m) return NULL;
3297 1710744 : A = rowpermute(A, y1);
3298 1710744 : B = rowpermute(B, y1);
3299 : /* a is square and invertible */
3300 1710739 : ncol = lg(B);
3301 1710739 : worker = snm_closure(is_entry("_ZM_gauss_worker"), mkvec2(A,B));
3302 1710743 : av = avma;
3303 1710743 : for (k = 1;; k <<= 1)
3304 207358 : {
3305 : pari_timer ti;
3306 1918101 : gen_inccrt_i("ZM_gauss", worker, NULL, (k+1)>>1 , m,
3307 : &S, &H, &mod, nmV_chinese_center, FpM_center);
3308 1918092 : gerepileall(av, 2, &H, &mod);
3309 1918105 : if (DEBUGLEVEL >= 4) timer_start(&ti);
3310 1918105 : Hr = FpM_ratlift_parallel(H, mod, NULL);
3311 1918095 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_gauss: ratlift (%ld)",!!Hr);
3312 1918100 : if (Hr)
3313 : {
3314 : GEN MH, c;
3315 1765049 : MH = ZM_mul(A, Q_remove_denom(Hr, &c));
3316 1765035 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_gauss: QM_mul");
3317 1765043 : if (ZM_equal(MH, c ? ZM_Z_mul(B, c): B)) break;
3318 : }
3319 : }
3320 1710737 : return iscol ? gel(Hr, 1): Hr;
3321 : }
3322 :
3323 : GEN
3324 1710829 : ZM_gauss(GEN A, GEN B)
3325 : {
3326 1710829 : pari_sp av = avma;
3327 1710829 : GEN C = ZM_gauss_i(A,B);
3328 1710831 : return C ? gc_GEN(av, C): NULL;
3329 : }
3330 :
3331 : GEN
3332 18235 : QM_ker(GEN M)
3333 : {
3334 18235 : pari_sp av = avma;
3335 18235 : long l = lg(M)-1;
3336 18235 : if (l==0) return cgetg(1, t_MAT);
3337 18193 : if (lgcols(M)==1) return matid(l);
3338 17276 : return gc_GEN(av, ZM_ker_i(row_Q_primpart(M)));
3339 : }
3340 :
3341 : /* x a ZM. Return a multiple of the determinant of the lattice generated by
3342 : * the columns of x. From Algorithm 2.2.6 in GTM138 */
3343 : GEN
3344 49964 : detint(GEN A)
3345 : {
3346 49964 : if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
3347 49964 : RgM_check_ZM(A, "detint");
3348 49964 : return ZM_detmult(A);
3349 : }
3350 : GEN
3351 166009 : ZM_detmult(GEN A)
3352 : {
3353 166009 : pari_sp av1, av = avma;
3354 : GEN B, c, v, piv;
3355 166009 : long rg, i, j, k, m, n = lg(A) - 1;
3356 :
3357 166009 : if (!n) return gen_1;
3358 166009 : m = nbrows(A);
3359 166009 : if (n < m) return gen_0;
3360 165932 : c = zero_zv(m);
3361 165932 : av1 = avma;
3362 165932 : B = zeromatcopy(m,m);
3363 165932 : v = cgetg(m+1, t_COL);
3364 165931 : piv = gen_1; rg = 0;
3365 718481 : for (k=1; k<=n; k++)
3366 : {
3367 718467 : GEN pivprec = piv;
3368 718467 : long t = 0;
3369 5335418 : for (i=1; i<=m; i++)
3370 : {
3371 4616953 : pari_sp av2 = avma;
3372 : GEN vi;
3373 4616953 : if (c[i]) continue;
3374 :
3375 2667960 : vi = mulii(piv, gcoeff(A,i,k));
3376 28320857 : for (j=1; j<=m; j++)
3377 25652897 : if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
3378 2667960 : if (!t && signe(vi)) t = i;
3379 2667960 : gel(v,i) = gc_INT(av2, vi);
3380 : }
3381 718465 : if (!t) continue;
3382 : /* at this point c[t] = 0 */
3383 :
3384 718374 : if (++rg >= m) { /* full rank; mostly done */
3385 165917 : GEN det = gel(v,t); /* last on stack */
3386 165917 : if (++k > n)
3387 165786 : det = absi(det);
3388 : else
3389 : {
3390 : /* improve further; at this point c[i] is set for all i != t */
3391 131 : gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
3392 418 : for ( ; k<=n; k++)
3393 286 : det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
3394 : }
3395 165918 : return gc_INT(av, det);
3396 : }
3397 :
3398 552457 : piv = gel(v,t);
3399 4450535 : for (i=1; i<=m; i++)
3400 : {
3401 : GEN mvi;
3402 3898078 : if (c[i] || i == t) continue;
3403 :
3404 1949040 : gcoeff(B,t,i) = mvi = negi(gel(v,i));
3405 22981595 : for (j=1; j<=m; j++)
3406 21032555 : if (c[j]) /* implies j != t */
3407 : {
3408 5711493 : pari_sp av2 = avma;
3409 5711493 : GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
3410 5711492 : if (rg > 1) z = diviiexact(z, pivprec);
3411 5711492 : gcoeff(B,j,i) = gc_INT(av2, z);
3412 : }
3413 : }
3414 552457 : c[t] = k;
3415 552457 : if (gc_needed(av,1))
3416 : {
3417 0 : if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
3418 0 : gerepileall(av1, 2, &piv,&B); v = zerovec(m);
3419 : }
3420 : }
3421 14 : return gc_const(av, gen_0);
3422 : }
3423 :
3424 : /* Reduce x modulo (invertible) y */
3425 : GEN
3426 9119 : closemodinvertible(GEN x, GEN y)
3427 : {
3428 9119 : return gmul(y, ground(RgM_solve(y,x)));
3429 : }
3430 : GEN
3431 7 : reducemodinvertible(GEN x, GEN y)
3432 : {
3433 7 : return gsub(x, closemodinvertible(x,y));
3434 : }
3435 : GEN
3436 0 : reducemodlll(GEN x,GEN y)
3437 : {
3438 0 : return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
3439 : }
3440 :
3441 : /*******************************************************************/
3442 : /* */
3443 : /* KERNEL of an m x n matrix */
3444 : /* return n - rk(x) linearly independent vectors */
3445 : /* */
3446 : /*******************************************************************/
3447 : static GEN
3448 28 : RgM_deplin_i(GEN x0)
3449 : {
3450 28 : pari_sp av = avma, av2;
3451 28 : long i, j, k, nl, nc = lg(x0)-1;
3452 : GEN D, x, y, c, l, d, ck;
3453 :
3454 28 : if (!nc) return NULL;
3455 28 : nl = nbrows(x0);
3456 28 : c = zero_zv(nl);
3457 28 : l = cgetg(nc+1, t_VECSMALL); /* not initialized */
3458 28 : av2 = avma;
3459 28 : x = RgM_shallowcopy(x0);
3460 28 : d = const_vec(nl, gen_1); /* pivot list */
3461 28 : ck = NULL; /* gcc -Wall */
3462 98 : for (k=1; k<=nc; k++)
3463 : {
3464 91 : ck = gel(x,k);
3465 196 : for (j=1; j<k; j++)
3466 : {
3467 105 : GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
3468 420 : for (i=1; i<=nl; i++)
3469 315 : if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
3470 : }
3471 :
3472 91 : i = gauss_get_pivot_NZ(x, NULL, k, c);
3473 91 : if (i > nl) break;
3474 70 : if (gc_needed(av,1))
3475 : {
3476 0 : if (DEBUGMEM>1) pari_warn(warnmem,"deplin k = %ld/%ld",k,nc);
3477 0 : gerepileall(av2, 2, &x, &d);
3478 0 : ck = gel(x,k);
3479 : }
3480 70 : gel(d,k) = gel(ck,i);
3481 70 : c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
3482 : }
3483 28 : if (k > nc) return gc_NULL(av);
3484 21 : if (k == 1) { set_avma(av); return scalarcol_shallow(gen_1,nc); }
3485 21 : y = cgetg(nc+1,t_COL);
3486 21 : gel(y,1) = gcopy(gel(ck, l[1]));
3487 49 : for (D=gel(d,1),j=2; j<k; j++)
3488 : {
3489 28 : gel(y,j) = gmul(gel(ck, l[j]), D);
3490 28 : D = gmul(D, gel(d,j));
3491 : }
3492 21 : gel(y,j) = gneg(D);
3493 21 : for (j++; j<=nc; j++) gel(y,j) = gen_0;
3494 21 : y = primitive_part(y, &c);
3495 21 : return c? gerepileupto(av, y): gc_GEN(av, y);
3496 : }
3497 : static GEN
3498 0 : RgV_deplin(GEN v)
3499 : {
3500 0 : pari_sp av = avma;
3501 0 : long n = lg(v)-1;
3502 0 : GEN y, p = NULL;
3503 0 : if (n <= 1)
3504 : {
3505 0 : if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
3506 0 : return cgetg(1, t_COL);
3507 : }
3508 0 : if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
3509 0 : v = primpart(mkvec2(gel(v,1),gel(v,2)));
3510 0 : if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
3511 0 : y = zerocol(n);
3512 0 : gel(y,1) = gneg(gel(v,2));
3513 0 : gel(y,2) = gcopy(gel(v,1));
3514 0 : return gerepileupto(av, y);
3515 :
3516 : }
3517 :
3518 : static GEN
3519 105 : RgM_deplin_FpM(GEN x, GEN p)
3520 : {
3521 105 : pari_sp av = avma;
3522 : ulong pp;
3523 105 : x = RgM_Fp_init3(x, p, &pp);
3524 105 : switch(pp)
3525 : {
3526 35 : case 0:
3527 35 : x = FpM_ker_gen(x,p,1);
3528 35 : if (!x) return gc_NULL(av);
3529 21 : x = FpC_center(x,p,shifti(p,-1));
3530 21 : break;
3531 14 : case 2:
3532 14 : x = F2m_ker_sp(x,1);
3533 14 : if (!x) return gc_NULL(av);
3534 7 : x = F2c_to_ZC(x); break;
3535 0 : case 3:
3536 0 : x = F3m_ker_sp(x,1);
3537 0 : if (!x) return gc_NULL(av);
3538 0 : x = F3c_to_ZC(x); break;
3539 56 : default:
3540 56 : x = Flm_ker_sp(x,pp,1);
3541 56 : if (!x) return gc_NULL(av);
3542 35 : x = Flv_center(x, pp, pp>>1);
3543 35 : x = zc_to_ZC(x);
3544 35 : break;
3545 : }
3546 63 : return gerepileupto(av, x);
3547 : }
3548 :
3549 : /* FIXME: implement direct modular ZM_deplin ? */
3550 : static GEN
3551 119 : QM_deplin(GEN M)
3552 : {
3553 119 : pari_sp av = avma;
3554 119 : long l = lg(M)-1;
3555 : GEN k;
3556 119 : if (l==0) return NULL;
3557 84 : if (lgcols(M)==1) return col_ei(l, 1);
3558 84 : k = ZM_ker_i(row_Q_primpart(M));
3559 84 : if (lg(k)== 1) return gc_NULL(av);
3560 70 : return gc_GEN(av, gel(k,1));
3561 : }
3562 :
3563 : static GEN
3564 49 : RgM_deplin_FqM(GEN x, GEN pol, GEN p)
3565 : {
3566 49 : pari_sp av = avma;
3567 49 : GEN b, T = RgX_to_FpX(pol, p);
3568 49 : if (signe(T) == 0) pari_err_OP("deplin",x,pol);
3569 49 : b = FqM_deplin(RgM_to_FqM(x, T, p), T, p);
3570 49 : if (!b) return gc_NULL(av);
3571 35 : return gerepileupto(av, b);
3572 : }
3573 :
3574 : /* Returns gen_0 instead of NULL for 'no fast algorithm'. NULL is already
3575 : * reserved for 'not invertible' */
3576 : static GEN
3577 385 : RgM_deplin_fast(GEN x)
3578 : {
3579 : GEN p, pol;
3580 385 : long pa, t = RgM_type(x, &p,&pol,&pa);
3581 385 : switch(t)
3582 : {
3583 119 : case t_INT: /* fall through */
3584 119 : case t_FRAC: return QM_deplin(x);
3585 84 : case t_FFELT: return FFM_deplin(x, pol);
3586 105 : case t_INTMOD: return RgM_deplin_FpM(x, p);
3587 49 : case RgX_type_code(t_POLMOD, t_INTMOD):
3588 49 : return RgM_deplin_FqM(x, pol, p);
3589 28 : default: return gen_0;
3590 : }
3591 : }
3592 :
3593 : static GEN
3594 385 : RgM_deplin(GEN x)
3595 : {
3596 385 : GEN z = RgM_deplin_fast(x);
3597 385 : if (z!= gen_0) return z;
3598 28 : return RgM_deplin_i(x);
3599 : }
3600 :
3601 : GEN
3602 385 : deplin(GEN x)
3603 : {
3604 385 : switch(typ(x))
3605 : {
3606 385 : case t_MAT:
3607 : {
3608 385 : GEN z = RgM_deplin(x);
3609 385 : if (z) return z;
3610 147 : return cgetg(1, t_COL);
3611 : }
3612 0 : case t_VEC: return RgV_deplin(x);
3613 0 : default: pari_err_TYPE("deplin",x);
3614 : }
3615 : return NULL;/*LCOV_EXCL_LINE*/
3616 : }
3617 :
3618 : /*******************************************************************/
3619 : /* */
3620 : /* GAUSS REDUCTION OF MATRICES (m lines x n cols) */
3621 : /* (kernel, image, complementary image, rank) */
3622 : /* */
3623 : /*******************************************************************/
3624 : /* return the transform of x under a standard Gauss pivot.
3625 : * x0 is a reference point when guessing whether x[i,j] ~ 0
3626 : * (iff x[i,j] << x0[i,j])
3627 : * Set r = dim ker(x). d[k] contains the index of the first nonzero pivot
3628 : * in column k */
3629 : static GEN
3630 1271 : gauss_pivot_ker(GEN x, GEN *dd, long *rr, pivot_fun pivot, GEN data)
3631 : {
3632 : GEN c, d, p;
3633 : pari_sp av;
3634 : long i, j, k, r, t, n, m;
3635 :
3636 1271 : n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
3637 1271 : m=nbrows(x); r=0;
3638 1271 : x = RgM_shallowcopy(x);
3639 1271 : c = zero_zv(m);
3640 1271 : d = cgetg(n+1,t_VECSMALL);
3641 1271 : av=avma;
3642 7475 : for (k=1; k<=n; k++)
3643 : {
3644 6204 : j = pivot(x, data, k, c);
3645 6204 : if (j > m)
3646 : {
3647 1463 : r++; d[k]=0;
3648 6496 : for(j=1; j<k; j++)
3649 5033 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
3650 : }
3651 : else
3652 : { /* pivot for column k on row j */
3653 4741 : c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
3654 4741 : gcoeff(x,j,k) = gen_m1;
3655 : /* x[j,] /= - x[j,k] */
3656 24169 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3657 42136 : for (t=1; t<=m; t++)
3658 37395 : if (t!=j)
3659 : { /* x[t,] -= 1 / x[j,k] x[j,] */
3660 32654 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3661 32654 : if (gequal0(p)) continue;
3662 86934 : for (i=k+1; i<=n; i++)
3663 69470 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
3664 17464 : if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
3665 : }
3666 : }
3667 : }
3668 1271 : *dd=d; *rr=r; return x;
3669 : }
3670 :
3671 : /* r = dim ker(x).
3672 : * Returns d:
3673 : * d[k] != 0 contains the index of a nonzero pivot in column k
3674 : * d[k] == 0 if column k is a linear combination of the (k-1) first ones */
3675 : GEN
3676 167668 : RgM_pivots(GEN x0, long *rr, pivot_fun pivot, GEN data)
3677 : {
3678 : GEN x, c, d, p;
3679 167668 : long i, j, k, r, t, m, n = lg(x0)-1;
3680 : pari_sp av;
3681 :
3682 167668 : if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
3683 152425 : if (!n) { *rr = 0; return NULL; }
3684 :
3685 152425 : d = cgetg(n+1, t_VECSMALL);
3686 152425 : x = RgM_shallowcopy(x0);
3687 152424 : m = nbrows(x); r = 0;
3688 152424 : c = zero_zv(m);
3689 152442 : av = avma;
3690 929065 : for (k=1; k<=n; k++)
3691 : {
3692 776641 : j = pivot(x, data, k, c);
3693 776644 : if (j > m) { r++; d[k] = 0; }
3694 : else
3695 : {
3696 291762 : c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
3697 1883375 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3698 :
3699 1054764 : for (t=1; t<=m; t++)
3700 763023 : if (!c[t]) /* no pivot on that line yet */
3701 : {
3702 257220 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3703 4119762 : for (i=k+1; i<=n; i++)
3704 3862537 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
3705 257225 : if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
3706 : }
3707 2175244 : for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
3708 : }
3709 : }
3710 152424 : *rr = r; return gc_const((pari_sp)d, d);
3711 : }
3712 :
3713 : static long
3714 4213910 : ZM_count_0_cols(GEN M)
3715 : {
3716 4213910 : long i, l = lg(M), n = 0;
3717 18110589 : for (i = 1; i < l; i++)
3718 13896678 : if (ZV_equal0(gel(M,i))) n++;
3719 4213911 : return n;
3720 : }
3721 :
3722 : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
3723 : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
3724 : GEN
3725 4227593 : ZM_pivots(GEN M0, long *rr)
3726 : {
3727 4227593 : GEN d, dbest = NULL;
3728 : long m, mm, n, nn, i, imax, rmin, rbest, zc;
3729 4227593 : int beenthere = 0;
3730 4227593 : pari_sp av, av0 = avma;
3731 : forprime_t S;
3732 :
3733 4227593 : rbest = n = lg(M0)-1;
3734 4227593 : if (n == 0) { *rr = 0; return NULL; }
3735 4213909 : zc = ZM_count_0_cols(M0);
3736 4213904 : if (n == zc) { *rr = zc; return zero_zv(n); }
3737 :
3738 4213772 : m = nbrows(M0);
3739 4213769 : rmin = maxss(zc, n-m);
3740 4213768 : init_modular_small(&S);
3741 4213796 : if (n <= m) { nn = n; mm = m; } else { nn = m; mm = n; }
3742 4213796 : imax = (nn < 16)? 1: (nn < 64)? 2: 3; /* heuristic */
3743 :
3744 : for(;;)
3745 0 : {
3746 : GEN row, col, M, KM, IM, RHS, X, cX;
3747 : long rk;
3748 4237029 : for (av = avma, i = 0;; set_avma(av), i++)
3749 23233 : {
3750 4237029 : ulong p = u_forprime_next(&S);
3751 : long rp;
3752 4237023 : if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
3753 4237023 : d = Flm_pivots(ZM_to_Flm(M0, p), p, &rp, 1);
3754 4237016 : if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
3755 45079 : if (rp < rbest) { /* save best r so far */
3756 21871 : rbest = rp;
3757 21871 : guncloneNULL(dbest);
3758 21871 : dbest = gclone(d);
3759 21871 : if (beenthere) break;
3760 : }
3761 45079 : if (!beenthere && i >= imax) break;
3762 : }
3763 21846 : beenthere = 1;
3764 : /* Dubious case: there is (probably) a non trivial kernel */
3765 21846 : indexrank_all(m,n, rbest, dbest, &row, &col);
3766 21846 : M = rowpermute(vecpermute(M0, col), row);
3767 21846 : rk = n - rbest; /* (probable) dimension of image */
3768 21846 : if (n > m) M = shallowtrans(M);
3769 21846 : IM = vecslice(M,1,rk);
3770 21846 : KM = vecslice(M,rk+1, nn);
3771 21846 : M = rowslice(IM, 1,rk); /* square maximal rank */
3772 21846 : X = ZM_gauss(M, rowslice(KM, 1,rk));
3773 21846 : RHS = rowslice(KM,rk+1,mm);
3774 21846 : M = rowslice(IM,rk+1,mm);
3775 21846 : X = Q_remove_denom(X, &cX);
3776 21846 : if (cX) RHS = ZM_Z_mul(RHS, cX);
3777 21846 : if (ZM_equal(ZM_mul(M, X), RHS)) { d = vecsmall_copy(dbest); goto END; }
3778 0 : set_avma(av);
3779 : }
3780 4213783 : END:
3781 4213783 : *rr = rbest; guncloneNULL(dbest);
3782 4213782 : return gerepileuptoleaf(av0, d);
3783 : }
3784 :
3785 : /* compute ker(x) */
3786 : static GEN
3787 1271 : ker_aux(GEN x, pivot_fun fun, GEN data)
3788 : {
3789 1271 : pari_sp av = avma;
3790 : GEN d,y;
3791 : long i,j,k,r,n;
3792 :
3793 1271 : x = gauss_pivot_ker(x,&d,&r, fun, data);
3794 1271 : if (!r) retgc_const(av, cgetg(1, t_MAT));
3795 1211 : n = lg(x)-1; y=cgetg(r+1,t_MAT);
3796 2674 : for (j=k=1; j<=r; j++,k++)
3797 : {
3798 1463 : GEN p = cgetg(n+1,t_COL);
3799 :
3800 5586 : gel(y,j) = p; while (d[k]) k++;
3801 6496 : for (i=1; i<k; i++)
3802 5033 : if (d[i])
3803 : {
3804 4641 : GEN p1=gcoeff(x,d[i],k);
3805 4641 : gel(p,i) = gcopy(p1); gunclone(p1);
3806 : }
3807 : else
3808 392 : gel(p,i) = gen_0;
3809 2541 : gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
3810 : }
3811 1211 : return gerepileupto(av,y);
3812 : }
3813 :
3814 : static GEN
3815 553 : RgM_ker_FpM(GEN x, GEN p)
3816 : {
3817 553 : pari_sp av = avma;
3818 : ulong pp;
3819 553 : x = RgM_Fp_init3(x, p, &pp);
3820 553 : switch(pp)
3821 : {
3822 35 : case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
3823 21 : case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
3824 77 : case 3: x = F3m_to_mod(F3m_ker_sp(x,0)); break;
3825 420 : default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
3826 : }
3827 553 : return gerepileupto(av, x);
3828 : }
3829 :
3830 : static GEN
3831 91 : RgM_ker_FqM(GEN x, GEN pol, GEN p)
3832 : {
3833 91 : pari_sp av = avma;
3834 91 : GEN b, T = RgX_to_FpX(pol, p);
3835 91 : if (signe(T) == 0) pari_err_OP("ker",x,pol);
3836 84 : b = FqM_ker(RgM_to_FqM(x, T, p), T, p);
3837 84 : return gerepileupto(av, FqM_to_mod(b, T, p));
3838 : }
3839 :
3840 : static GEN
3841 10668 : RgM_ker_fast(GEN x, pivot_fun *fun, GEN *data)
3842 : {
3843 : GEN p, pol;
3844 10668 : long pa, t = RgM_type(x, &p,&pol,&pa);
3845 10668 : set_pivot_fun(fun, data, t, x, p);
3846 10668 : switch(t)
3847 : {
3848 9079 : case t_INT: /* fall through */
3849 9079 : case t_FRAC: return QM_ker(x);
3850 63 : case t_FFELT: return FFM_ker(x, pol);
3851 553 : case t_INTMOD: return RgM_ker_FpM(x, p);
3852 91 : case RgX_type_code(t_POLMOD, t_INTMOD):
3853 91 : return RgM_ker_FqM(x, pol, p);
3854 882 : default: return NULL;
3855 : }
3856 : }
3857 :
3858 : GEN
3859 10668 : ker(GEN x)
3860 : {
3861 : pivot_fun fun;
3862 10668 : GEN data, b = RgM_ker_fast(x, &fun, &data);
3863 10661 : if (b) return b;
3864 882 : return ker_aux(x, fun, data);
3865 : }
3866 :
3867 : GEN
3868 46221 : matker0(GEN x,long flag)
3869 : {
3870 46221 : if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
3871 46221 : if (!flag) return ker(x);
3872 45934 : RgM_check_ZM(x, "matker");
3873 45934 : return ZM_ker(x);
3874 : }
3875 :
3876 : static GEN
3877 525 : RgM_image_FpM(GEN x, GEN p)
3878 : {
3879 525 : pari_sp av = avma;
3880 : ulong pp;
3881 525 : x = RgM_Fp_init(x, p, &pp);
3882 525 : switch(pp)
3883 : {
3884 28 : case 0: x = FpM_to_mod(FpM_image(x,p),p); break;
3885 7 : case 2: x = F2m_to_mod(F2m_image(x)); break;
3886 490 : default:x = Flm_to_mod(Flm_image(x,pp), pp); break;
3887 : }
3888 525 : return gerepileupto(av, x);
3889 : }
3890 :
3891 : static GEN
3892 35 : RgM_image_FqM(GEN x, GEN pol, GEN p)
3893 : {
3894 35 : pari_sp av = avma;
3895 35 : GEN b, T = RgX_to_FpX(pol, p);
3896 35 : if (signe(T) == 0) pari_err_OP("image",x,pol);
3897 28 : b = FqM_image(RgM_to_FqM(x, T, p), T, p);
3898 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
3899 : }
3900 :
3901 : GEN
3902 6181 : QM_image_shallow(GEN A)
3903 : {
3904 6181 : A = vec_Q_primpart(A);
3905 6181 : return vecpermute(A, ZM_indeximage(A));
3906 : }
3907 : GEN
3908 5411 : QM_image(GEN A)
3909 : {
3910 5411 : pari_sp av = avma;
3911 5411 : return gc_GEN(av, QM_image_shallow(A));
3912 : }
3913 :
3914 : static GEN
3915 6034 : RgM_image_fast(GEN x, pivot_fun *fun, GEN *data)
3916 : {
3917 : GEN p, pol;
3918 6034 : long pa, t = RgM_type(x, &p,&pol,&pa);
3919 6034 : set_pivot_fun(fun, data, t, x, p);
3920 6034 : switch(t)
3921 : {
3922 5411 : case t_INT: /* fall through */
3923 5411 : case t_FRAC: return QM_image(x);
3924 49 : case t_FFELT: return FFM_image(x, pol);
3925 525 : case t_INTMOD: return RgM_image_FpM(x, p);
3926 35 : case RgX_type_code(t_POLMOD, t_INTMOD):
3927 35 : return RgM_image_FqM(x, pol, p);
3928 14 : default: return NULL;
3929 : }
3930 : }
3931 :
3932 : GEN
3933 6034 : image(GEN x)
3934 : {
3935 : pivot_fun fun;
3936 : GEN d, M, data;
3937 : long r;
3938 :
3939 6034 : if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
3940 6034 : M = RgM_image_fast(x, &fun, &data);
3941 6027 : if (M) return M;
3942 14 : d = RgM_pivots(x, &r, fun, data); /* d left on stack for efficiency */
3943 14 : return image_from_pivot(x,d,r);
3944 : }
3945 :
3946 : static GEN
3947 84 : imagecompl_aux(GEN d, long r)
3948 : {
3949 84 : GEN y = cgetg(r+1,t_VECSMALL);
3950 : long j, i;
3951 126 : for (i = j = 1; j<=r; i++)
3952 42 : if (!d[i]) y[j++] = i;
3953 84 : return y;
3954 : }
3955 : GEN
3956 84 : imagecompl(GEN x)
3957 : {
3958 84 : pari_sp av = avma;
3959 : GEN data, d;
3960 : long r;
3961 : pivot_fun fun;
3962 :
3963 84 : if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
3964 84 : init_pivot_list(x); set_pivot_fun_all(&fun, &data, x);
3965 84 : d = RgM_pivots(x, &r, fun, data); /* if (!d) then r = 0 */
3966 84 : set_avma(av); return imagecompl_aux(d, r);
3967 : }
3968 : GEN
3969 0 : ZM_imagecompl(GEN x)
3970 : {
3971 0 : pari_sp av = avma;
3972 : GEN d;
3973 : long r;
3974 :
3975 0 : init_pivot_list(x);
3976 0 : d = ZM_pivots(x, &r); /* if (!d) then r = 0 */
3977 0 : set_avma(av); return imagecompl_aux(d, r);
3978 : }
3979 :
3980 : static GEN
3981 28 : RgM_RgC_invimage_FpC(GEN A, GEN y, GEN p)
3982 : {
3983 28 : pari_sp av = avma;
3984 : ulong pp;
3985 : GEN x;
3986 28 : A = RgM_Fp_init(A,p,&pp);
3987 28 : switch(pp)
3988 : {
3989 7 : case 0:
3990 7 : y = RgC_to_FpC(y,p);
3991 7 : x = FpM_FpC_invimage(A, y, p);
3992 7 : return x ? gerepileupto(av, FpC_to_mod(x,p)): NULL;
3993 7 : case 2:
3994 7 : y = RgV_to_F2v(y);
3995 7 : x = F2m_F2c_invimage(A, y);
3996 7 : return x ? gerepileupto(av, F2c_to_mod(x)): NULL;
3997 14 : default:
3998 14 : y = RgV_to_Flv(y,pp);
3999 14 : x = Flm_Flc_invimage(A, y, pp);
4000 14 : return x ? gerepileupto(av, Flc_to_mod(x,pp)): NULL;
4001 : }
4002 : }
4003 :
4004 : /* Returns gen_0 instead of NULL for 'no fast algorithm'. NULL is already
4005 : * reserved for 'not invertible' */
4006 : static GEN
4007 3654 : RgM_RgC_invimage_fast(GEN x, GEN y)
4008 : {
4009 : GEN p, pol;
4010 3654 : long pa, t = RgM_RgC_type(x, y, &p,&pol,&pa);
4011 3654 : switch(t)
4012 : {
4013 28 : case t_INTMOD: return RgM_RgC_invimage_FpC(x, y, p);
4014 63 : case t_FFELT: return FFM_FFC_invimage(x, y, pol);
4015 3563 : default: return gen_0;
4016 : }
4017 : }
4018 :
4019 : GEN
4020 3759 : RgM_RgC_invimage(GEN A, GEN y)
4021 : {
4022 3759 : pari_sp av = avma;
4023 3759 : long i, l = lg(A);
4024 : GEN M, x, t;
4025 3759 : if (l==1) return NULL;
4026 3654 : if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
4027 3654 : M = RgM_RgC_invimage_fast(A, y);
4028 3654 : if (!M) return gc_NULL(av);
4029 3633 : if (M != gen_0) return M;
4030 3563 : M = ker(shallowconcat(A, y));
4031 3563 : i = lg(M)-1;
4032 3563 : if (!i) return gc_NULL(av);
4033 :
4034 3304 : x = gel(M,i); t = gel(x,l);
4035 3304 : if (gequal0(t)) return gc_NULL(av);
4036 :
4037 1862 : t = gneg_i(t); setlg(x,l);
4038 1862 : return gerepileupto(av, RgC_Rg_div(x, t));
4039 : }
4040 :
4041 : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
4042 : * if no solution exist */
4043 : GEN
4044 3920 : inverseimage(GEN m, GEN v)
4045 : {
4046 : GEN y;
4047 3920 : if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
4048 3920 : switch(typ(v))
4049 : {
4050 3682 : case t_COL:
4051 3682 : y = RgM_RgC_invimage(m,v);
4052 3682 : return y? y: cgetg(1,t_COL);
4053 238 : case t_MAT:
4054 238 : y = RgM_invimage(m, v);
4055 238 : return y? y: cgetg(1,t_MAT);
4056 : }
4057 0 : pari_err_TYPE("inverseimage",v);
4058 : return NULL;/*LCOV_EXCL_LINE*/
4059 : }
4060 :
4061 : static GEN
4062 84 : RgM_invimage_FpM(GEN A, GEN B, GEN p)
4063 : {
4064 84 : pari_sp av = avma;
4065 : ulong pp;
4066 : GEN x;
4067 84 : A = RgM_Fp_init(A,p,&pp);
4068 84 : switch(pp)
4069 : {
4070 35 : case 0:
4071 35 : B = RgM_to_FpM(B,p);
4072 35 : x = FpM_invimage_gen(A, B, p);
4073 35 : return x ? gerepileupto(av, FpM_to_mod(x, p)): x;
4074 7 : case 2:
4075 7 : B = RgM_to_F2m(B);
4076 7 : x = F2m_invimage_i(A, B);
4077 7 : return x ? gerepileupto(av, F2m_to_mod(x)): x;
4078 42 : default:
4079 42 : B = RgM_to_Flm(B,pp);
4080 42 : x = Flm_invimage_i(A, B, pp);
4081 42 : return x ? gerepileupto(av, Flm_to_mod(x, pp)): x;
4082 : }
4083 : }
4084 :
4085 : /* Returns gen_0 instead of NULL for 'no fast algorithm'. NULL is already
4086 : * reserved for 'not invertible' */
4087 : static GEN
4088 364 : RgM_invimage_fast(GEN x, GEN y)
4089 : {
4090 : GEN p, pol;
4091 364 : long pa, t = RgM_type2(x, y, &p,&pol,&pa);
4092 364 : switch(t)
4093 : {
4094 84 : case t_INTMOD: return RgM_invimage_FpM(x, y, p);
4095 105 : case t_FFELT: return FFM_invimage(x, y, pol);
4096 175 : default: return gen_0;
4097 : }
4098 : }
4099 :
4100 : /* find Z such that A Z = B. Return NULL if no solution */
4101 : GEN
4102 364 : RgM_invimage(GEN A, GEN B)
4103 : {
4104 364 : pari_sp av = avma;
4105 : GEN d, x, X, Y;
4106 364 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
4107 364 : X = RgM_invimage_fast(A, B);
4108 364 : if (!X) return gc_NULL(av);
4109 252 : if (X != gen_0) return X;
4110 175 : x = ker(shallowconcat(RgM_neg(A), B));
4111 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
4112 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
4113 : * Y has at least nB columns and full rank */
4114 175 : nY = lg(x)-1;
4115 175 : if (nY < nB) return gc_NULL(av);
4116 161 : Y = rowslice(x, nA+1, nA+nB); /* nB rows */
4117 161 : d = cgetg(nB+1, t_VECSMALL);
4118 721 : for (i = nB, j = nY; i >= 1; i--, j--)
4119 : {
4120 805 : for (; j>=1; j--)
4121 756 : if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
4122 609 : if (!j) return gc_NULL(av);
4123 : }
4124 : /* reduce to the case Y square, upper triangular with 1s on diagonal */
4125 112 : Y = vecpermute(Y, d);
4126 112 : x = vecpermute(x, d);
4127 112 : X = rowslice(x, 1, nA);
4128 112 : return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
4129 : }
4130 :
4131 : static GEN
4132 70 : RgM_suppl_FpM(GEN x, GEN p)
4133 : {
4134 70 : pari_sp av = avma;
4135 : ulong pp;
4136 70 : x = RgM_Fp_init(x, p, &pp);
4137 70 : switch(pp)
4138 : {
4139 21 : case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
4140 14 : case 2: x = F2m_to_mod(F2m_suppl(x)); break;
4141 35 : default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
4142 : }
4143 70 : return gerepileupto(av, x);
4144 : }
4145 :
4146 : static GEN
4147 175 : RgM_suppl_fast(GEN x, pivot_fun *fun, GEN *data)
4148 : {
4149 : GEN p, pol;
4150 175 : long pa, t = RgM_type(x,&p,&pol,&pa);
4151 175 : set_pivot_fun(fun, data, t, x, p);
4152 175 : switch(t)
4153 : {
4154 70 : case t_INTMOD: return RgM_suppl_FpM(x, p);
4155 35 : case t_FFELT: return FFM_suppl(x, pol);
4156 70 : default: return NULL;
4157 : }
4158 : }
4159 :
4160 : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
4161 : * whose first k columns are given by x. If rank(x) < k, undefined result. */
4162 : GEN
4163 175 : suppl(GEN x)
4164 : {
4165 175 : pari_sp av = avma;
4166 : pivot_fun fun;
4167 : GEN d, M, data;
4168 : long r;
4169 175 : if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
4170 175 : M = RgM_suppl_fast(x, &fun, &data);
4171 175 : if (M) return M;
4172 70 : init_suppl(x);
4173 70 : d = RgM_pivots(x, &r, fun, data); set_avma(av);
4174 70 : return get_suppl(x,d,nbrows(x),r,&col_ei);
4175 : }
4176 :
4177 : GEN
4178 7 : image2(GEN x)
4179 : {
4180 7 : pari_sp av = avma;
4181 : long k, n, i;
4182 : GEN A, B;
4183 :
4184 7 : if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
4185 7 : if (lg(x) == 1) return cgetg(1,t_MAT);
4186 7 : A = ker(x); k = lg(A)-1;
4187 7 : if (!k) { set_avma(av); return gcopy(x); }
4188 7 : A = suppl(A); n = lg(A)-1;
4189 7 : B = cgetg(n-k+1, t_MAT);
4190 21 : for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
4191 7 : return gerepileupto(av, B);
4192 : }
4193 :
4194 : GEN
4195 217 : matimage0(GEN x,long flag)
4196 : {
4197 217 : switch(flag)
4198 : {
4199 210 : case 0: return image(x);
4200 7 : case 1: return image2(x);
4201 0 : default: pari_err_FLAG("matimage");
4202 : }
4203 : return NULL; /* LCOV_EXCL_LINE */
4204 : }
4205 :
4206 : static long
4207 126 : RgM_rank_FpM(GEN x, GEN p)
4208 : {
4209 126 : pari_sp av = avma;
4210 : ulong pp;
4211 : long r;
4212 126 : x = RgM_Fp_init(x,p,&pp);
4213 126 : switch(pp)
4214 : {
4215 28 : case 0: r = FpM_rank(x,p); break;
4216 63 : case 2: r = F2m_rank(x); break;
4217 35 : default:r = Flm_rank(x,pp); break;
4218 : }
4219 126 : return gc_long(av, r);
4220 : }
4221 :
4222 : static long
4223 49 : RgM_rank_FqM(GEN x, GEN pol, GEN p)
4224 : {
4225 49 : pari_sp av = avma;
4226 : long r;
4227 49 : GEN T = RgX_to_FpX(pol, p);
4228 49 : if (signe(T) == 0) pari_err_OP("rank",x,pol);
4229 42 : r = FqM_rank(RgM_to_FqM(x, T, p), T, p);
4230 42 : return gc_long(av,r);
4231 : }
4232 :
4233 : static long
4234 371 : RgM_rank_fast(GEN x, pivot_fun *fun, GEN *data)
4235 : {
4236 : GEN p, pol;
4237 371 : long pa, t = RgM_type(x,&p,&pol,&pa);
4238 371 : set_pivot_fun(fun, data, t, x, p);
4239 371 : switch(t)
4240 : {
4241 98 : case t_INT: return ZM_rank(x);
4242 21 : case t_FRAC: return QM_rank(x);
4243 126 : case t_INTMOD: return RgM_rank_FpM(x, p);
4244 70 : case t_FFELT: return FFM_rank(x, pol);
4245 49 : case RgX_type_code(t_POLMOD, t_INTMOD):
4246 49 : return RgM_rank_FqM(x, pol, p);
4247 7 : default: return -1;
4248 : }
4249 : }
4250 :
4251 : long
4252 371 : rank(GEN x)
4253 : {
4254 371 : pari_sp av = avma;
4255 : pivot_fun fun;
4256 : GEN data;
4257 : long r;
4258 :
4259 371 : if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
4260 371 : r = RgM_rank_fast(x, &fun, &data);
4261 364 : if (r >= 0) return r;
4262 7 : (void)RgM_pivots(x, &r, fun, data);
4263 7 : return gc_long(av, lg(x)-1 - r);
4264 : }
4265 :
4266 : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
4267 : * followed by the missing indices */
4268 : static GEN
4269 43692 : perm_complete(GEN d, long n)
4270 : {
4271 43692 : GEN y = cgetg(n+1, t_VECSMALL);
4272 43692 : long i, j = 1, k = n, l = lg(d);
4273 43692 : pari_sp av = avma;
4274 43692 : char *T = stack_calloc(n+1);
4275 214790 : for (i = 1; i < l; i++) T[d[i]] = 1;
4276 418608 : for (i = 1; i <= n; i++)
4277 374916 : if (T[i]) y[j++] = i; else y[k--] = i;
4278 43692 : return gc_const(av, y);
4279 : }
4280 :
4281 : /* n = dim x, r = dim Ker(x), d from RgM_pivots */
4282 : static GEN
4283 6181 : indeximage0(long n, long r, GEN d)
4284 : {
4285 : long i, j;
4286 : GEN v;
4287 :
4288 6181 : r = n - r; /* now r = dim Im(x) */
4289 6181 : v = cgetg(r+1,t_VECSMALL);
4290 34419 : if (d) for (i=j=1; j<=n; j++)
4291 28238 : if (d[j]) v[i++] = j;
4292 6181 : return v;
4293 : }
4294 : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from RgM_pivots */
4295 : static void
4296 21846 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
4297 : {
4298 21846 : GEN IR = indexrank0(n, r, d);
4299 21846 : *prow = perm_complete(gel(IR,1), m);
4300 21846 : *pcol = perm_complete(gel(IR,2), n);
4301 21846 : }
4302 :
4303 : static GEN
4304 28 : RgM_indexrank_FpM(GEN x, GEN p)
4305 : {
4306 28 : pari_sp av = avma;
4307 : ulong pp;
4308 : GEN r;
4309 28 : x = RgM_Fp_init(x,p,&pp);
4310 28 : switch(pp)
4311 : {
4312 7 : case 0: r = FpM_indexrank(x,p); break;
4313 7 : case 2: r = F2m_indexrank(x); break;
4314 14 : default: r = Flm_indexrank(x,pp); break;
4315 : }
4316 28 : return gerepileupto(av, r);
4317 : }
4318 :
4319 : static GEN
4320 0 : RgM_indexrank_FqM(GEN x, GEN pol, GEN p)
4321 : {
4322 0 : pari_sp av = avma;
4323 0 : GEN r, T = RgX_to_FpX(pol, p);
4324 0 : if (signe(T) == 0) pari_err_OP("indexrank",x,pol);
4325 0 : r = FqM_indexrank(RgM_to_FqM(x, T, p), T, p);
4326 0 : return gerepileupto(av, r);
4327 : }
4328 :
4329 : static GEN
4330 77599 : RgM_indexrank_fast(GEN x, pivot_fun *fun, GEN *data)
4331 : {
4332 : GEN p, pol;
4333 77599 : long pa, t = RgM_type(x,&p,&pol,&pa);
4334 77600 : set_pivot_fun(fun, data, t, x, p);
4335 77600 : switch(t)
4336 : {
4337 406 : case t_INT: return ZM_indexrank(x);
4338 1344 : case t_FRAC: return QM_indexrank(x);
4339 28 : case t_INTMOD: return RgM_indexrank_FpM(x, p);
4340 21 : case t_FFELT: return FFM_indexrank(x, pol);
4341 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
4342 0 : return RgM_indexrank_FqM(x, pol, p);
4343 75801 : default: return NULL;
4344 : }
4345 : }
4346 :
4347 : GEN
4348 77599 : indexrank(GEN x)
4349 : {
4350 : pari_sp av;
4351 : pivot_fun fun;
4352 : long r;
4353 : GEN d, data;
4354 77599 : if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
4355 77599 : d = RgM_indexrank_fast(x, &fun, &data);
4356 77600 : if (d) return d;
4357 75801 : av = avma;
4358 75801 : init_pivot_list(x); d = RgM_pivots(x, &r, fun, data);
4359 75800 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4360 : }
4361 :
4362 : GEN
4363 6181 : ZM_indeximage(GEN x) {
4364 6181 : pari_sp av = avma;
4365 : long r;
4366 : GEN d;
4367 6181 : init_pivot_list(x); d = ZM_pivots(x,&r);
4368 6181 : set_avma(av); return indeximage0(lg(x)-1, r, d);
4369 : }
4370 : long
4371 2225351 : ZM_rank(GEN x) {
4372 2225351 : pari_sp av = avma;
4373 : long r;
4374 2225351 : (void)ZM_pivots(x,&r);
4375 2225353 : return gc_long(av, lg(x)-1-r);
4376 : }
4377 : GEN
4378 1742190 : ZM_indexrank(GEN x) {
4379 1742190 : pari_sp av = avma;
4380 : long r;
4381 : GEN d;
4382 1742190 : init_pivot_list(x); d = ZM_pivots(x,&r);
4383 1742196 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4384 : }
4385 :
4386 : long
4387 21 : QM_rank(GEN x)
4388 : {
4389 21 : pari_sp av = avma;
4390 21 : return gc_long(av, ZM_rank(Q_primpart(x)));
4391 : }
4392 :
4393 : GEN
4394 1344 : QM_indexrank(GEN x)
4395 : {
4396 1344 : pari_sp av = avma;
4397 1344 : return gerepileupto(av, ZM_indexrank(Q_primpart(x)));
4398 : }
4399 :
4400 : /*******************************************************************/
4401 : /* */
4402 : /* ZabM */
4403 : /* */
4404 : /*******************************************************************/
4405 :
4406 : static GEN
4407 1276 : FpXM_ratlift(GEN a, GEN q)
4408 : {
4409 : GEN B, y;
4410 1276 : long i, j, l = lg(a), n;
4411 1276 : B = sqrti(shifti(q,-1));
4412 1276 : y = cgetg(l, t_MAT);
4413 1276 : if (l==1) return y;
4414 1276 : n = lgcols(a);
4415 3059 : for (i=1; i<l; i++)
4416 : {
4417 2404 : GEN yi = cgetg(n, t_COL);
4418 32311 : for (j=1; j<n; j++)
4419 : {
4420 30528 : GEN v = FpX_ratlift(gmael(a,i,j), q, B, B, NULL);
4421 30528 : if (!v) return NULL;
4422 29907 : gel(yi, j) = RgX_renormalize(v);
4423 : }
4424 1783 : gel(y,i) = yi;
4425 : }
4426 655 : return y;
4427 : }
4428 :
4429 : static GEN
4430 4485 : FlmV_recover_pre(GEN a, GEN M, ulong p, ulong pi, long sv)
4431 : {
4432 4485 : GEN a1 = gel(a,1);
4433 4485 : long i, j, k, l = lg(a1), n, lM = lg(M);
4434 4485 : GEN v = cgetg(lM, t_VECSMALL);
4435 4485 : GEN y = cgetg(l, t_MAT);
4436 4485 : if (l==1) return y;
4437 4485 : n = lgcols(a1);
4438 22521 : for (i=1; i<l; i++)
4439 : {
4440 18036 : GEN yi = cgetg(n, t_COL);
4441 347613 : for (j=1; j<n; j++)
4442 : {
4443 4676595 : for (k=1; k<lM; k++) uel(v,k) = umael(gel(a,k),i,j);
4444 329577 : gel(yi, j) = Flm_Flc_mul_pre_Flx(M, v, p, pi, sv);
4445 : }
4446 18036 : gel(y,i) = yi;
4447 : }
4448 4485 : return y;
4449 : }
4450 :
4451 : static GEN
4452 0 : FlkM_inv(GEN M, GEN P, ulong p)
4453 : {
4454 0 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4455 0 : GEN R = Flx_roots_pre(P, p, pi);
4456 0 : long l = lg(R), i;
4457 0 : GEN W = Flv_invVandermonde(R, 1UL, p);
4458 0 : GEN V = cgetg(l, t_VEC);
4459 0 : for(i=1; i<l; i++)
4460 : {
4461 0 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, PI);
4462 0 : GEN H = Flm_inv_sp(FlxM_eval_powers_pre(M, pows, p, pi), NULL, p);
4463 0 : if (!H) return NULL;
4464 0 : gel(V, i) = H;
4465 : }
4466 0 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4467 : }
4468 :
4469 : static GEN
4470 3209 : FlkM_adjoint(GEN M, GEN P, ulong p)
4471 : {
4472 3209 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4473 3209 : GEN R = Flx_roots_pre(P, p, pi);
4474 3209 : long l = lg(R), i;
4475 3209 : GEN W = Flv_invVandermonde(R, 1UL, p);
4476 3209 : GEN V = cgetg(l, t_VEC);
4477 15577 : for(i=1; i<l; i++)
4478 : {
4479 12368 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, PI);
4480 12368 : gel(V, i) = Flm_adjoint(FlxM_eval_powers_pre(M, pows, p, pi), p);
4481 : }
4482 3209 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4483 : }
4484 :
4485 : static GEN
4486 1985 : ZabM_inv_slice(GEN A, GEN Q, GEN P, GEN *mod)
4487 : {
4488 1985 : pari_sp av = avma;
4489 1985 : long i, n = lg(P)-1, w = varn(Q);
4490 : GEN H, T;
4491 1985 : if (n == 1)
4492 : {
4493 1554 : ulong p = uel(P,1);
4494 1554 : GEN Qp = ZX_to_Flx(Q, p);
4495 1554 : GEN Ap = ZXM_to_FlxM(A, p, get_Flx_var(Qp));
4496 1554 : GEN Hp = FlkM_adjoint(Ap, Qp, p);
4497 1554 : Hp = gerepileupto(av, FlxM_to_ZXM(Hp));
4498 1554 : *mod = utoipos(p); return Hp;
4499 : }
4500 431 : T = ZV_producttree(P);
4501 431 : A = ZXM_nv_mod_tree(A, P, T, w);
4502 431 : Q = ZX_nv_mod_tree(Q, P, T);
4503 431 : H = cgetg(n+1, t_VEC);
4504 2086 : for(i=1; i <= n; i++)
4505 : {
4506 1655 : ulong p = P[i];
4507 1655 : GEN a = gel(A,i), q = gel(Q, i);
4508 1655 : gel(H,i) = FlkM_adjoint(a, q, p);
4509 : }
4510 431 : H = nxMV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
4511 431 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
4512 : }
4513 :
4514 : GEN
4515 1985 : ZabM_inv_worker(GEN P, GEN A, GEN Q)
4516 : {
4517 1985 : GEN V = cgetg(3, t_VEC);
4518 1985 : gel(V,1) = ZabM_inv_slice(A, Q, P, &gel(V,2));
4519 1985 : return V;
4520 : }
4521 :
4522 : static GEN
4523 5509 : vecnorml1(GEN x)
4524 60508 : { pari_APPLY_same(gnorml1_fake(gel(x,i))); }
4525 :
4526 : static GEN
4527 1827 : ZabM_true_Hadamard(GEN a)
4528 : {
4529 1827 : pari_sp av = avma;
4530 1827 : long n = lg(a)-1, i;
4531 : GEN B;
4532 1827 : if (n == 0) return gen_1;
4533 1827 : if (n == 1) return gnorml1_fake(gcoeff(a,1,1));
4534 1183 : B = gen_1;
4535 6692 : for (i = 1; i <= n; i++)
4536 5509 : B = gmul(B, gnorml2(RgC_gtofp(vecnorml1(gel(a,i)),DEFAULTPREC)));
4537 1183 : return gc_INT(av, ceil_safe(sqrtr_abs(B)));
4538 : }
4539 :
4540 : GEN
4541 1827 : ZabM_inv(GEN A, GEN Q, long n, GEN *pt_den)
4542 : {
4543 1827 : pari_sp av = avma;
4544 : forprime_t S;
4545 : GEN bnd, H, D, d, mod, worker;
4546 1827 : if (lg(A) == 1)
4547 : {
4548 0 : if (pt_den) *pt_den = gen_1;
4549 0 : return cgetg(1, t_MAT);
4550 : }
4551 1827 : bnd = ZabM_true_Hadamard(A);
4552 1827 : worker = snm_closure(is_entry("_ZabM_inv_worker"), mkvec2(A, Q));
4553 1827 : u_forprime_arith_init(&S, HIGHBIT+1, ULONG_MAX, 1, n);
4554 1827 : H = gen_crt("ZabM_inv", worker, &S, NULL, expi(bnd), 0, &mod,
4555 : nxMV_chinese_center, FpXM_center);
4556 1827 : D = RgMrow_RgC_mul(H, gel(A,1), 1);
4557 1827 : D = ZX_rem(D, Q);
4558 1827 : d = Z_content(mkvec2(H, D));
4559 1827 : if (d)
4560 : {
4561 518 : D = ZX_Z_divexact(D, d);
4562 518 : H = Q_div_to_int(H, d);
4563 : }
4564 1827 : if (!pt_den) return gerepileupto(av, H);
4565 1827 : *pt_den = D; return gc_all(av, 2, &H, pt_den);
4566 : }
4567 :
4568 : GEN
4569 0 : ZabM_inv_ratlift(GEN M, GEN P, long n, GEN *pden)
4570 : {
4571 0 : pari_sp av2, av = avma;
4572 : GEN q, H;
4573 0 : ulong m = LONG_MAX>>1;
4574 0 : ulong p= 1 + m - (m % n);
4575 0 : long lM = lg(M);
4576 0 : if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
4577 :
4578 0 : av2 = avma;
4579 0 : H = NULL;
4580 : for(;;)
4581 0 : {
4582 : GEN Hp, Pp, Mp, Hr;
4583 0 : do p += n; while(!uisprime(p));
4584 0 : Pp = ZX_to_Flx(P, p);
4585 0 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4586 0 : Hp = FlkM_inv(Mp, Pp, p);
4587 0 : if (!Hp) continue;
4588 0 : if (!H)
4589 : {
4590 0 : H = ZXM_init_CRT(Hp, degpol(P)-1, p);
4591 0 : q = utoipos(p);
4592 : }
4593 : else
4594 0 : ZXM_incremental_CRT(&H, Hp, &q, p);
4595 0 : Hr = FpXM_ratlift(H, q);
4596 0 : if (DEBUGLEVEL>5) err_printf("ZabM_inv mod %ld (ratlift=%ld)\n", p,!!Hr);
4597 0 : if (Hr) {/* DONE ? */
4598 0 : GEN Hl = Q_remove_denom(Hr, pden);
4599 0 : GEN MH = ZXQM_mul(Hl, M, P);
4600 0 : if (*pden)
4601 0 : { if (RgM_isscalar(MH, *pden)) { H = Hl; break; }}
4602 : else
4603 0 : { if (RgM_isidentity(MH)) { H = Hl; *pden = gen_1; break; } }
4604 : }
4605 :
4606 0 : if (gc_needed(av,2))
4607 : {
4608 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_inv");
4609 0 : gerepileall(av2, 2, &H, &q);
4610 : }
4611 : }
4612 0 : return gc_all(av, 2, &H, pden);
4613 : }
4614 :
4615 : static GEN
4616 1276 : FlkM_ker(GEN M, GEN P, ulong p)
4617 : {
4618 1276 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4619 1276 : GEN R = Flx_roots_pre(P, p, pi);
4620 1276 : long l = lg(R), i, dP = degpol(P), r;
4621 : GEN M1, K, D;
4622 1276 : GEN W = Flv_invVandermonde(R, 1UL, p);
4623 1276 : GEN V = cgetg(l, t_VEC);
4624 1276 : M1 = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,1), dP, p, PI), p, pi);
4625 1276 : K = Flm_ker_sp(M1, p, 2);
4626 1276 : r = lg(gel(K,1)); D = gel(K,2);
4627 1276 : gel(V, 1) = gel(K,1);
4628 2652 : for(i=2; i<l; i++)
4629 : {
4630 1376 : GEN Mi = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,i), dP, p, PI), p, pi);
4631 1376 : GEN K = Flm_ker_sp(Mi, p, 2);
4632 1376 : if (lg(gel(K,1)) != r || !zv_equal(D, gel(K,2))) return NULL;
4633 1376 : gel(V, i) = gel(K,1);
4634 : }
4635 1276 : return mkvec2(FlmV_recover_pre(V, W, p, pi, P[1]), D);
4636 : }
4637 :
4638 : static int
4639 655 : ZabM_ker_check(GEN M, GEN H, ulong p, GEN P, long n)
4640 : {
4641 : GEN pow;
4642 655 : long j, l = lg(H);
4643 : ulong pi, r;
4644 3899 : do p += n; while(!uisprime(p));
4645 655 : pi = get_Fl_red(p);
4646 655 : P = ZX_to_Flx(P, p);
4647 655 : r = Flx_oneroot_pre(P, p, pi);
4648 655 : pow = Fl_powers_pre(r, degpol(P),p, (p & HIGHMASK)? pi: 0);
4649 655 : M = ZXM_to_FlxM(M, p, P[1]); M = FlxM_eval_powers_pre(M, pow, p, pi);
4650 655 : H = ZXM_to_FlxM(H, p, P[1]); H = FlxM_eval_powers_pre(H, pow, p, pi);
4651 2178 : for (j = 1; j < l; j++)
4652 1555 : if (!zv_equal0(Flm_Flc_mul_pre(M, gel(H,j), p, pi))) return 0;
4653 623 : return 1;
4654 : }
4655 :
4656 : GEN
4657 623 : ZabM_ker(GEN M, GEN P, long n)
4658 : {
4659 623 : pari_sp av = avma;
4660 : pari_timer ti;
4661 623 : GEN q, H = NULL, D = NULL;
4662 623 : ulong m = LONG_MAX>>1;
4663 623 : ulong p = 1 + m - (m % n);
4664 :
4665 623 : if (DEBUGLEVEL>5) timer_start(&ti);
4666 : for(;;)
4667 653 : {
4668 : GEN Kp, Hp, Dp, Pp, Mp, Hr;
4669 22341 : do p += n; while(!uisprime(p));
4670 1276 : Pp = ZX_to_Flx(P, p);
4671 1276 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4672 1276 : Kp = FlkM_ker(Mp, Pp, p);
4673 1276 : if (!Kp) continue;
4674 1276 : Hp = gel(Kp,1); Dp = gel(Kp,2);
4675 1276 : if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
4676 1276 : if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
4677 : {
4678 623 : H = ZXM_init_CRT(Hp, degpol(P)-1, p); D = Dp;
4679 623 : q = utoipos(p);
4680 : }
4681 : else
4682 653 : ZXM_incremental_CRT(&H, Hp, &q, p);
4683 1276 : Hr = FpXM_ratlift(H, q);
4684 1276 : if (DEBUGLEVEL>5) timer_printf(&ti,"ZabM_ker mod %ld (ratlift=%ld)", p,!!Hr);
4685 1276 : if (Hr) {/* DONE ? */
4686 655 : GEN Hl = vec_Q_primpart(Hr);
4687 655 : if (ZabM_ker_check(M, Hl, p, P, n)) { H = Hl; break; }
4688 : }
4689 :
4690 653 : if (gc_needed(av,2))
4691 : {
4692 4 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_ker");
4693 4 : gerepileall(av, 3, &H, &D, &q);
4694 : }
4695 : }
4696 623 : return gc_GEN(av, H);
4697 : }
4698 :
4699 : GEN
4700 2387 : ZabM_indexrank(GEN M, GEN P, long n)
4701 : {
4702 2387 : pari_sp av = avma;
4703 2387 : ulong m = LONG_MAX>>1;
4704 2387 : ulong p = 1+m-(m%n), D = degpol(P);
4705 2387 : long lM = lg(M), lmax = 0, c = 0;
4706 : GEN v;
4707 : for(;;)
4708 735 : {
4709 : GEN R, Pp, Mp, K;
4710 : ulong pi;
4711 : long l;
4712 61415 : do p += n; while (!uisprime(p));
4713 3122 : pi = (p & HIGHMASK)? get_Fl_red(p): 0;
4714 3122 : Pp = ZX_to_Flx(P, p);
4715 3122 : R = Flx_roots_pre(Pp, p, pi);
4716 3122 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4717 3122 : K = FlxM_eval_powers_pre(Mp, Fl_powers_pre(uel(R,1), D,p,pi), p,pi);
4718 3122 : v = Flm_indexrank(K, p);
4719 3122 : l = lg(gel(v,2));
4720 3122 : if (l == lM) break;
4721 980 : if (lmax >= 0 && l > lmax) { lmax = l; c = 0; } else c++;
4722 980 : if (c > 2)
4723 : { /* probably not maximal rank, expensive check */
4724 245 : lM -= lg(ZabM_ker(M, P, n))-1; /* actual rank (+1) */
4725 245 : if (lmax == lM) break;
4726 0 : lmax = -1; /* disable check */
4727 : }
4728 : }
4729 2387 : return gerepileupto(av, v);
4730 : }
4731 :
4732 : #if 0
4733 : GEN
4734 : ZabM_gauss(GEN M, GEN P, long n, GEN *den)
4735 : {
4736 : pari_sp av = avma;
4737 : GEN v, S, W;
4738 : v = ZabM_indexrank(M, P, n);
4739 : S = shallowmatextract(M,gel(v,1),gel(v,2));
4740 : W = ZabM_inv(S, P, n, den);
4741 : return gc_all(av,2,&W,den);
4742 : }
4743 : #endif
4744 :
4745 : GEN
4746 140 : ZabM_pseudoinv(GEN M, GEN P, long n, GEN *pv, GEN *den)
4747 : {
4748 140 : GEN v = ZabM_indexrank(M, P, n);
4749 140 : if (pv) *pv = v;
4750 140 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4751 140 : return ZabM_inv(M, P, n, den);
4752 : }
4753 : GEN
4754 5040 : ZM_pseudoinv(GEN M, GEN *pv, GEN *den)
4755 : {
4756 5040 : GEN v = ZM_indexrank(M);
4757 5040 : if (pv) *pv = v;
4758 5040 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4759 5040 : return ZM_inv(M, den);
4760 : }
4761 :
4762 : /*******************************************************************/
4763 : /* */
4764 : /* Structured Elimination */
4765 : /* */
4766 : /*******************************************************************/
4767 :
4768 : static void
4769 95967 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
4770 : {
4771 95967 : long lc = lg(c), k;
4772 95967 : iscol[i] = 0; (*rcol)--;
4773 891425 : for (k = 1; k < lc; ++k)
4774 : {
4775 795458 : Wrow[c[k]]--;
4776 795458 : if (Wrow[c[k]]==0) (*rrow)--;
4777 : }
4778 95967 : }
4779 :
4780 : static void
4781 7641 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long idx, long *rcol, long *rrow)
4782 : {
4783 : long i, j;
4784 7641 : long nbcol = lg(iscol)-1, last;
4785 : do
4786 : {
4787 9568 : last = 0;
4788 16914392 : for (i = 1; i <= nbcol; ++i)
4789 16904824 : if (iscol[i])
4790 : {
4791 9073526 : GEN c = idx ? gmael(M, i, idx): gel(M,i);
4792 9073526 : long lc = lg(c);
4793 83822176 : for (j = 1; j < lc; ++j)
4794 74766687 : if (Wrow[c[j]] == 1)
4795 : {
4796 18037 : rem_col(c, i, iscol, Wrow, rcol, rrow);
4797 18037 : last=1; break;
4798 : }
4799 : }
4800 9568 : } while (last);
4801 7641 : }
4802 :
4803 : static GEN
4804 7448 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
4805 : {
4806 7448 : long nbcol = lg(iscol)-1;
4807 : long i, j, m, last;
4808 : GEN per;
4809 20552 : for (m = 2, last=0; !last ; m++)
4810 : {
4811 25077349 : for (i = 1; i <= nbcol; ++i)
4812 : {
4813 25064245 : wcol[i] = 0;
4814 25064245 : if (iscol[i])
4815 : {
4816 13863572 : GEN c = gmael(M, i, 1);
4817 13863572 : long lc = lg(c);
4818 123886864 : for (j = 1; j < lc; ++j)
4819 110023292 : if (Wrow[c[j]] == m) { wcol[i]++; last = 1; }
4820 : }
4821 : }
4822 : }
4823 7448 : per = vecsmall_indexsort(wcol);
4824 7448 : *w = wcol[per[nbcol]];
4825 7448 : return per;
4826 : }
4827 :
4828 : /* M is a RgMs with nbrow rows, A a list of row indices.
4829 : Eliminate rows of M with a single entry that do not belong to A,
4830 : and the corresponding columns. Also eliminate columns until #colums=#rows.
4831 : Return pcol and prow:
4832 : pcol is a map from the new columns indices to the old one.
4833 : prow is a map from the old rows indices to the new one (0 if removed).
4834 : */
4835 :
4836 : void
4837 147 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4838 : {
4839 147 : long i, j, k, lA = lg(A);
4840 147 : GEN prow = cgetg(nbrow+1, t_VECSMALL);
4841 147 : GEN pcol = zero_zv(nbcol);
4842 147 : pari_sp av = avma;
4843 147 : long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
4844 147 : GEN iscol = const_vecsmall(nbcol, 1);
4845 147 : GEN Wrow = zero_zv(nbrow);
4846 147 : GEN wcol = cgetg(nbcol+1, t_VECSMALL);
4847 147 : pari_sp av2 = avma;
4848 110397 : for (i = 1; i <= nbcol; ++i)
4849 : {
4850 110250 : GEN F = gmael(M, i, 1);
4851 110250 : long l = lg(F)-1;
4852 924676 : for (j = 1; j <= l; ++j) Wrow[F[j]]++;
4853 : }
4854 147 : for (j = 1; j < lA; ++j)
4855 : {
4856 0 : if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
4857 0 : Wrow[A[j]] = -1;
4858 : }
4859 228354 : for (i = 1; i <= nbrow; ++i)
4860 228207 : if (Wrow[i]) rrow++;
4861 147 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow);
4862 147 : if (rcol < rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
4863 7595 : while (rcol > rrow)
4864 : {
4865 : long w;
4866 7448 : GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
4867 85378 : for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
4868 77930 : rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
4869 7448 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow); set_avma(av2);
4870 : }
4871 110397 : for (j = 1, i = 1; i <= nbcol; ++i)
4872 110250 : if (iscol[i]) pcol[j++] = i;
4873 147 : setlg(pcol,j);
4874 228354 : for (k = 1, i = 1; i <= nbrow; ++i) prow[i] = Wrow[i]? k++: 0;
4875 147 : *p_col = pcol; *p_row = prow; set_avma(av);
4876 : }
4877 :
4878 : void
4879 0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4880 0 : { RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row); }
4881 :
4882 : GEN
4883 46 : F2Ms_colelim(GEN M, long nbrow)
4884 : {
4885 46 : long i,j, nbcol = lg(M)-1, rcol = nbcol, rrow = 0;
4886 46 : GEN pcol = zero_zv(nbcol);
4887 46 : pari_sp av = avma;
4888 46 : GEN iscol = const_vecsmall(nbcol, 1), Wrow = zero_zv(nbrow);
4889 77470 : for (i = 1; i <= nbcol; ++i)
4890 : {
4891 77424 : GEN F = gel(M, i);
4892 77424 : long l = lg(F)-1;
4893 1431938 : for (j = 1; j <= l; ++j) Wrow[F[j]]++;
4894 : }
4895 46 : rem_singleton(M, iscol, Wrow, 0, &rcol, &rrow);
4896 77470 : for (j = 1, i = 1; i <= nbcol; ++i)
4897 77424 : if (iscol[i]) pcol[j++] = i;
4898 46 : fixlg(pcol,j); return gc_const(av, pcol);
4899 : }
4900 :
4901 : /*******************************************************************/
4902 : /* */
4903 : /* EIGENVECTORS */
4904 : /* (independent eigenvectors, sorted by increasing eigenvalue) */
4905 : /* */
4906 : /*******************************************************************/
4907 : /* assume x is square of dimension > 0 */
4908 : static int
4909 53 : RgM_is_symmetric_cx(GEN x, long bit)
4910 : {
4911 53 : pari_sp av = avma;
4912 53 : long i, j, l = lg(x);
4913 239 : for (i = 1; i < l; i++)
4914 708 : for (j = 1; j < i; j++)
4915 : {
4916 522 : GEN a = gcoeff(x,i,j), b = gcoeff(x,j,i), c = gsub(a,b);
4917 522 : if (!gequal0(c) && gexpo(c) - gexpo(a) > -bit) return gc_long(av,0);
4918 : }
4919 21 : return gc_long(av,1);
4920 : }
4921 : static GEN
4922 53 : eigen_err(int exact, GEN x, long flag, long prec)
4923 : {
4924 53 : pari_sp av = avma;
4925 : GEN y;
4926 53 : if (RgM_is_symmetric_cx(x, prec - 10))
4927 : { /* approximately symmetric: recover */
4928 21 : x = jacobi(x, prec); if (flag) return x;
4929 14 : return gc_GEN(av, gel(x,2));
4930 : }
4931 32 : if (!exact) x = bestappr(x, NULL);
4932 32 : y = mateigen(x, flag, precdbl(prec));
4933 32 : if (exact)
4934 18 : y = gprec_wtrunc(y, prec);
4935 14 : else if (flag)
4936 7 : y = mkvec2(RgV_gtofp(gel(y,1), prec), RgM_gtofp(gel(y,2), prec));
4937 : else
4938 7 : y = RgM_gtofp(y, prec);
4939 32 : return gc_GEN(av, y);
4940 : }
4941 : GEN
4942 144 : mateigen(GEN x, long flag, long prec)
4943 : {
4944 : GEN y, R, T;
4945 144 : long k, l, ex, n = lg(x);
4946 : int exact;
4947 144 : pari_sp av = avma;
4948 :
4949 144 : if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
4950 144 : if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
4951 144 : if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
4952 144 : if (n == 1)
4953 : {
4954 14 : if (flag) retmkvec2(cgetg(1,t_COL), cgetg(1,t_MAT));
4955 7 : return cgetg(1,t_MAT);
4956 : }
4957 130 : if (n == 2)
4958 : {
4959 14 : if (flag) retmkvec2(mkcolcopy(gcoeff(x,1,1)), matid(1));
4960 7 : return matid(1);
4961 : }
4962 :
4963 116 : ex = 16 - prec;
4964 116 : T = charpoly(x,0);
4965 116 : exact = RgX_is_QX(T);
4966 116 : if (exact)
4967 : {
4968 74 : T = ZX_radical( Q_primpart(T) );
4969 74 : R = nfrootsQ(T); settyp(R, t_COL);
4970 74 : if (lg(R)-1 < degpol(T))
4971 : { /* add missing complex roots */
4972 60 : GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
4973 60 : R = shallowconcat(R, r);
4974 : }
4975 : }
4976 : else
4977 : {
4978 42 : GEN r1, v = vectrunc_init(lg(T));
4979 : long e;
4980 42 : R = cleanroots(T,prec);
4981 42 : r1 = NULL;
4982 266 : for (k = 1; k < lg(R); k++)
4983 : {
4984 224 : GEN r2 = gel(R,k), r = grndtoi(r2, &e);
4985 224 : if (e < ex) r2 = r;
4986 224 : if (r1)
4987 : {
4988 182 : r = gsub(r1,r2);
4989 182 : if (gequal0(r) || gexpo(r) < ex) continue;
4990 : }
4991 182 : vectrunc_append(v, r2);
4992 182 : r1 = r2;
4993 : }
4994 42 : R = v;
4995 : }
4996 : /* R = distinct complex roots of charpoly(x) */
4997 116 : l = lg(R); y = cgetg(l, t_VEC);
4998 452 : for (k = 1; k < l; k++)
4999 : {
5000 389 : GEN M = RgM_Rg_sub_shallow(x, gel(R,k));
5001 389 : GEN F = ker_aux(M, gauss_get_pivot_max, x);
5002 389 : long d = lg(F)-1;
5003 389 : if (!d) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5004 336 : gel(y,k) = F;
5005 336 : if (flag) gel(R,k) = const_col(d, gel(R,k));
5006 : }
5007 63 : y = shallowconcat1(y);
5008 63 : if (lg(y) > n) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5009 : /* lg(y) < n if x is not diagonalizable */
5010 63 : if (flag) y = mkvec2(shallowconcat1(R), y);
5011 63 : return gc_GEN(av,y);
5012 : }
5013 : GEN
5014 0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
5015 :
5016 : /*******************************************************************/
5017 : /* */
5018 : /* DETERMINANT */
5019 : /* */
5020 : /*******************************************************************/
5021 :
5022 : GEN
5023 26593 : det0(GEN a,long flag)
5024 : {
5025 26593 : switch(flag)
5026 : {
5027 26579 : case 0: return det(a);
5028 14 : case 1: return det2(a);
5029 0 : default: pari_err_FLAG("matdet");
5030 : }
5031 : return NULL; /* LCOV_EXCL_LINE */
5032 : }
5033 :
5034 : /* M a 2x2 matrix, returns det(M) */
5035 : static GEN
5036 94283 : RgM_det2(GEN M)
5037 : {
5038 94283 : pari_sp av = avma;
5039 94283 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5040 94283 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5041 94283 : return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
5042 : }
5043 : /* M a 2x2 ZM, returns det(M) */
5044 : static GEN
5045 8722 : ZM_det2(GEN M)
5046 : {
5047 8722 : pari_sp av = avma;
5048 8722 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5049 8722 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5050 8722 : return gc_INT(av, subii(mulii(a,d), mulii(b, c)));
5051 : }
5052 : /* M a 3x3 ZM, return det(M) */
5053 : static GEN
5054 100472 : ZM_det3(GEN M)
5055 : {
5056 100472 : pari_sp av = avma;
5057 100472 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
5058 100472 : GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
5059 100472 : GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
5060 100472 : GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
5061 100472 : if (signe(g))
5062 : {
5063 66202 : t = mulii(subii(mulii(b,f), mulii(c,e)), g);
5064 66202 : D = addii(D, t);
5065 : }
5066 100472 : if (signe(h))
5067 : {
5068 77604 : t = mulii(subii(mulii(c,d), mulii(a,f)), h);
5069 77604 : D = addii(D, t);
5070 : }
5071 100472 : return gc_INT(av, D);
5072 : }
5073 :
5074 : static GEN
5075 58263 : det_simple_gauss(GEN a, pivot_fun pivot, GEN data)
5076 : {
5077 58263 : pari_sp av = avma;
5078 58263 : long i,j,k, s = 1, nbco = lg(a)-1;
5079 58263 : GEN p, x = gen_1;
5080 :
5081 58263 : a = RgM_shallowcopy(a);
5082 342099 : for (i=1; i<nbco; i++)
5083 : {
5084 283842 : k = pivot(a, data, i, NULL);
5085 283844 : if (k > nbco) return gc_GEN(av, gcoeff(a,i,i));
5086 283837 : if (k != i)
5087 : { /* exchange the lines s.t. k = i */
5088 1159169 : for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
5089 119043 : s = -s;
5090 : }
5091 283837 : p = gcoeff(a,i,i);
5092 :
5093 283837 : x = gmul(x,p);
5094 1786807 : for (k=i+1; k<=nbco; k++)
5095 : {
5096 1502976 : GEN m = gcoeff(a,i,k);
5097 1502976 : if (gequal0(m)) continue;
5098 :
5099 1067616 : m = gdiv(m,p);
5100 9951206 : for (j=i+1; j<=nbco; j++)
5101 8883592 : gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
5102 : }
5103 283831 : if (gc_needed(av,2))
5104 : {
5105 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5106 0 : gerepileall(av,2, &a,&x);
5107 : }
5108 : }
5109 58257 : if (s < 0) x = gneg_i(x);
5110 58257 : return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
5111 : }
5112 :
5113 : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
5114 : * Gauss-Bareiss. */
5115 : static GEN
5116 462 : det_bareiss(GEN a)
5117 : {
5118 462 : pari_sp av = avma;
5119 462 : long nbco = lg(a)-1,i,j,k,s = 1;
5120 : GEN p, pprec;
5121 :
5122 462 : a = RgM_shallowcopy(a);
5123 1337 : for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
5124 : {
5125 882 : int diveuc = (gequal1(pprec)==0);
5126 : GEN ci;
5127 :
5128 882 : p = gcoeff(a,i,i);
5129 882 : if (gequal0(p))
5130 : {
5131 14 : k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
5132 7 : if (k>nbco) return gc_GEN(av, p);
5133 0 : swap(gel(a,k), gel(a,i)); s = -s;
5134 0 : p = gcoeff(a,i,i);
5135 : }
5136 875 : ci = gel(a,i);
5137 2373 : for (k=i+1; k<=nbco; k++)
5138 : {
5139 1498 : GEN ck = gel(a,k), m = gel(ck,i);
5140 1498 : if (gequal0(m))
5141 : {
5142 7 : if (gequal1(p))
5143 : {
5144 0 : if (diveuc)
5145 0 : gel(a,k) = gdiv(gel(a,k), pprec);
5146 : }
5147 : else
5148 42 : for (j=i+1; j<=nbco; j++)
5149 : {
5150 35 : GEN p1 = gmul(p, gel(ck,j));
5151 35 : if (diveuc) p1 = gdiv(p1,pprec);
5152 35 : gel(ck,j) = p1;
5153 : }
5154 : }
5155 : else
5156 4662 : for (j=i+1; j<=nbco; j++)
5157 : {
5158 3171 : pari_sp av2 = avma;
5159 3171 : GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
5160 3171 : if (diveuc) p1 = gdiv(p1,pprec);
5161 3171 : gel(ck,j) = gerepileupto(av2, p1);
5162 : }
5163 1498 : if (gc_needed(av,2))
5164 : {
5165 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5166 0 : gerepileall(av,2, &a,&pprec);
5167 0 : ci = gel(a,i);
5168 0 : p = gcoeff(a,i,i);
5169 : }
5170 : }
5171 : }
5172 455 : p = gcoeff(a,nbco,nbco);
5173 455 : p = (s < 0)? gneg(p): gcopy(p);
5174 455 : return gerepileupto(av, p);
5175 : }
5176 :
5177 : /* count nonzero entries in col j, at most 'max' of them.
5178 : * Return their indices */
5179 : static GEN
5180 1470 : col_count_non_zero(GEN a, long j, long max)
5181 : {
5182 1470 : GEN v = cgetg(max+1, t_VECSMALL);
5183 1470 : GEN c = gel(a,j);
5184 1470 : long i, l = lg(a), k = 1;
5185 5614 : for (i = 1; i < l; i++)
5186 5376 : if (!gequal0(gel(c,i)))
5187 : {
5188 5110 : if (k > max) return NULL; /* fail */
5189 3878 : v[k++] = i;
5190 : }
5191 238 : setlg(v, k); return v;
5192 : }
5193 : /* count nonzero entries in row i, at most 'max' of them.
5194 : * Return their indices */
5195 : static GEN
5196 1456 : row_count_non_zero(GEN a, long i, long max)
5197 : {
5198 1456 : GEN v = cgetg(max+1, t_VECSMALL);
5199 1456 : long j, l = lg(a), k = 1;
5200 5558 : for (j = 1; j < l; j++)
5201 5334 : if (!gequal0(gcoeff(a,i,j)))
5202 : {
5203 5096 : if (k > max) return NULL; /* fail */
5204 3864 : v[k++] = j;
5205 : }
5206 224 : setlg(v, k); return v;
5207 : }
5208 :
5209 : static GEN det_develop(GEN a, long max, double bound);
5210 : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
5211 : static GEN
5212 406 : coeff_det(GEN a, long i, long j, long max, double bound)
5213 : {
5214 406 : GEN c = gcoeff(a, i, j);
5215 406 : c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
5216 406 : if (odd(i+j)) c = gneg(c);
5217 406 : return c;
5218 : }
5219 : /* a square t_MAT, 'bound' a rough upper bound for the number of
5220 : * multiplications we are willing to pay while developing rows/columns before
5221 : * switching to Gaussian elimination */
5222 : static GEN
5223 658 : det_develop(GEN M, long max, double bound)
5224 : {
5225 658 : pari_sp av = avma;
5226 658 : long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
5227 658 : GEN best = NULL;
5228 :
5229 658 : if (bound < 1.) return det_bareiss(M); /* too costly now */
5230 :
5231 434 : switch(n)
5232 : {
5233 0 : case 0: return gen_1;
5234 0 : case 1: return gcopy(gcoeff(M,1,1));
5235 14 : case 2: return RgM_det2(M);
5236 : }
5237 420 : if (max > ((n+2)>>1)) max = (n+2)>>1;
5238 1876 : for (j = 1; j <= n; j++)
5239 : {
5240 1470 : pari_sp av2 = avma;
5241 1470 : GEN v = col_count_non_zero(M, j, max);
5242 : long lv;
5243 1470 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5244 182 : if (lv == 1) { set_avma(av); return gen_0; }
5245 182 : if (lv == 2) {
5246 14 : set_avma(av);
5247 14 : return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
5248 : }
5249 168 : best = v; lbest = lv; best_col = j;
5250 : }
5251 1862 : for (i = 1; i <= n; i++)
5252 : {
5253 1456 : pari_sp av2 = avma;
5254 1456 : GEN v = row_count_non_zero(M, i, max);
5255 : long lv;
5256 1456 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5257 0 : if (lv == 1) { set_avma(av); return gen_0; }
5258 0 : if (lv == 2) {
5259 0 : set_avma(av);
5260 0 : return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
5261 : }
5262 0 : best = v; lbest = lv; best_row = i;
5263 : }
5264 406 : if (best_row)
5265 : {
5266 0 : double d = lbest-1;
5267 0 : GEN s = NULL;
5268 : long k;
5269 0 : bound /= d*d*d;
5270 0 : for (k = 1; k < lbest; k++)
5271 : {
5272 0 : GEN c = coeff_det(M, best_row, best[k], max, bound);
5273 0 : s = s? gadd(s, c): c;
5274 : }
5275 0 : return gerepileupto(av, s);
5276 : }
5277 406 : if (best_col)
5278 : {
5279 168 : double d = lbest-1;
5280 168 : GEN s = NULL;
5281 : long k;
5282 168 : bound /= d*d*d;
5283 560 : for (k = 1; k < lbest; k++)
5284 : {
5285 392 : GEN c = coeff_det(M, best[k], best_col, max, bound);
5286 392 : s = s? gadd(s, c): c;
5287 : }
5288 168 : return gerepileupto(av, s);
5289 : }
5290 238 : return det_bareiss(M);
5291 : }
5292 :
5293 : /* area of parallelogram bounded by (v1,v2) */
5294 : static GEN
5295 64380 : parallelogramarea(GEN v1, GEN v2)
5296 64380 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
5297 :
5298 : /* Square of Hadamard bound for det(a), a square matrix.
5299 : * Slight improvement: instead of using the column norms, use the area of
5300 : * the parallelogram formed by pairs of consecutive vectors */
5301 : GEN
5302 20017 : RgM_Hadamard(GEN a)
5303 : {
5304 20017 : pari_sp av = avma;
5305 20017 : long n = lg(a)-1, i;
5306 : GEN B;
5307 20017 : if (n == 0) return gen_1;
5308 20017 : if (n == 1) return gsqr(gcoeff(a,1,1));
5309 20017 : a = RgM_gtofp(a, LOWDEFAULTPREC);
5310 20017 : B = gen_1;
5311 84397 : for (i = 1; i <= n/2; i++)
5312 64380 : B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
5313 20017 : if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
5314 20017 : return gc_INT(av, ceil_safe(B));
5315 : }
5316 :
5317 : /* If B=NULL, assume B=A' */
5318 : static GEN
5319 20872 : ZM_det_slice(GEN A, GEN P, GEN *mod)
5320 : {
5321 20872 : pari_sp av = avma;
5322 20872 : long i, n = lg(P)-1;
5323 : GEN H, T;
5324 20872 : if (n == 1)
5325 : {
5326 0 : ulong Hp, p = uel(P,1);
5327 0 : GEN a = ZM_to_Flm(A, p);
5328 0 : Hp = Flm_det_sp(a, p);
5329 0 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
5330 : }
5331 20872 : T = ZV_producttree(P);
5332 20872 : A = ZM_nv_mod_tree(A, P, T);
5333 20872 : H = cgetg(n+1, t_VECSMALL);
5334 87546 : for(i=1; i <= n; i++)
5335 : {
5336 66674 : ulong p = P[i];
5337 66674 : GEN a = gel(A,i);
5338 66674 : H[i] = Flm_det_sp(a, p);
5339 : }
5340 20872 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
5341 20872 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
5342 : }
5343 :
5344 : GEN
5345 20872 : ZM_det_worker(GEN P, GEN A)
5346 : {
5347 20872 : GEN V = cgetg(3, t_VEC);
5348 20872 : gel(V,1) = ZM_det_slice(A, P, &gel(V,2));
5349 20872 : return V;
5350 : }
5351 :
5352 : GEN
5353 130793 : ZM_det(GEN M)
5354 : {
5355 : pari_sp av, av2;
5356 130793 : long n = lg(M)-1;
5357 : ulong p, Dp;
5358 : forprime_t S;
5359 : pari_timer ti;
5360 : GEN H, mod, h, q, worker;
5361 : #ifdef LONG_IS_64BIT
5362 112116 : const ulong PMAX = 18446744073709551557UL;
5363 : #else
5364 18677 : const ulong PMAX = 4294967291UL;
5365 : #endif
5366 :
5367 130793 : switch(n)
5368 : {
5369 7 : case 0: return gen_1;
5370 1575 : case 1: return icopy(gcoeff(M,1,1));
5371 8722 : case 2: return ZM_det2(M);
5372 100472 : case 3: return ZM_det3(M);
5373 : }
5374 20017 : if (DEBUGLEVEL>=4) timer_start(&ti);
5375 20017 : av = avma; h = RgM_Hadamard(M); /* |D| <= sqrt(h) */
5376 20017 : if (!signe(h)) { set_avma(av); return gen_0; }
5377 20017 : h = sqrti(h);
5378 20017 : if (lgefint(h) == 3 && (ulong)h[2] <= (PMAX >> 1))
5379 : { /* h < p/2 => direct result */
5380 7233 : p = PMAX;
5381 7233 : Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5382 7233 : set_avma(av);
5383 7233 : if (!Dp) return gen_0;
5384 7233 : return (Dp <= (p>>1))? utoipos(Dp): utoineg(p - Dp);
5385 : }
5386 12784 : q = gen_1; Dp = 1;
5387 12784 : init_modular_big(&S);
5388 12784 : p = 0; /* -Wall */
5389 12784 : while (cmpii(q, h) <= 0 && (p = u_forprime_next(&S)))
5390 : {
5391 12784 : av2 = avma; Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5392 12784 : set_avma(av2);
5393 12784 : if (Dp) break;
5394 0 : q = muliu(q, p);
5395 : }
5396 12784 : if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
5397 12784 : if (!Dp) { set_avma(av); return gen_0; }
5398 12784 : worker = snm_closure(is_entry("_ZM_det_worker"), mkvec(M));
5399 12784 : H = gen_crt("ZM_det", worker, &S, NULL, expi(h)+1, 0, &mod,
5400 : ZV_chinese, NULL);
5401 : /* H = det(M) modulo mod, (mod,D) = 1; |det(M) / D| <= h */
5402 12784 : H = Fp_center(H, mod, shifti(mod,-1));
5403 12784 : return gc_INT(av, H);
5404 : }
5405 :
5406 : static GEN
5407 1519 : RgM_det_FpM(GEN a, GEN p)
5408 : {
5409 1519 : pari_sp av = avma;
5410 : ulong pp, d;
5411 1519 : a = RgM_Fp_init(a,p,&pp);
5412 1519 : switch(pp)
5413 : {
5414 70 : case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
5415 14 : case 2: d = F2m_det_sp(a); break;
5416 1435 : default:d = Flm_det_sp(a, pp); break;
5417 : }
5418 1449 : set_avma(av); return mkintmodu(d, pp);
5419 : }
5420 :
5421 : static GEN
5422 42 : RgM_det_FqM(GEN x, GEN pol, GEN p)
5423 : {
5424 42 : pari_sp av = avma;
5425 42 : GEN b, T = RgX_to_FpX(pol, p);
5426 42 : if (signe(T) == 0) pari_err_OP("%",x,pol);
5427 42 : b = FqM_det(RgM_to_FqM(x, T, p), T, p);
5428 42 : if (!b) return gc_NULL(av);
5429 42 : return gc_GEN(av, mkpolmod(FpX_to_mod(b, p), FpX_to_mod(T, p)));
5430 : }
5431 :
5432 : static GEN
5433 33907 : RgM_det_fast(GEN x, pivot_fun *fun, GEN *data)
5434 : {
5435 : GEN p, pol;
5436 33907 : long pa, t = RgM_type(x, &p,&pol,&pa);
5437 33907 : set_pivot_fun(fun, data, t, x, p);
5438 33907 : switch(t)
5439 : {
5440 301 : case t_INT: return ZM_det(x);
5441 203 : case t_FRAC: return QM_det(x);
5442 63 : case t_FFELT: return FFM_det(x, pol);
5443 1519 : case t_INTMOD: return RgM_det_FpM(x, p);
5444 42 : case RgX_type_code(t_POLMOD, t_INTMOD): return RgM_det_FqM(x, pol, p);
5445 31779 : default: return NULL;
5446 : }
5447 : }
5448 :
5449 : static long
5450 252 : det_init_max(long n)
5451 : {
5452 252 : if (n > 100) return 0;
5453 252 : if (n > 50) return 1;
5454 252 : if (n > 30) return 4;
5455 252 : return 7;
5456 : }
5457 :
5458 : GEN
5459 246213 : det(GEN a)
5460 : {
5461 246213 : long n = lg(a)-1;
5462 : double B;
5463 : GEN data, b;
5464 : pivot_fun fun;
5465 :
5466 246213 : if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
5467 246213 : if (!n) return gen_1;
5468 246171 : if (n != nbrows(a)) pari_err_DIM("det");
5469 246164 : if (n == 1) return gcopy(gcoeff(a,1,1));
5470 69060 : if (n == 2) return RgM_det2(a);
5471 33907 : b = RgM_det_fast(a, &fun, &data);
5472 33907 : if (b) return b;
5473 31779 : if (data) return det_simple_gauss(a, fun, data);
5474 252 : B = (double)n; return det_develop(a, det_init_max(n), B*B*B);
5475 : }
5476 :
5477 : GEN
5478 134553 : det2(GEN a)
5479 : {
5480 134553 : long n = lg(a)-1;
5481 : GEN data;
5482 : pivot_fun fun;
5483 :
5484 134553 : if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
5485 134553 : if (!n) return gen_1;
5486 134553 : if (n != nbrows(a)) pari_err_DIM("det2");
5487 134553 : if (n == 1) return gcopy(gcoeff(a,1,1));
5488 85852 : if (n == 2) return RgM_det2(a);
5489 26736 : set_pivot_fun_all(&fun, &data, a);
5490 26736 : return det_simple_gauss(a, fun, data);
5491 : }
5492 :
5493 : GEN
5494 203 : QM_det(GEN M)
5495 : {
5496 203 : pari_sp av = avma;
5497 203 : GEN cM, pM = Q_primitive_part(M, &cM);
5498 203 : GEN b = ZM_det(pM);
5499 203 : if (cM) b = gmul(b, gpowgs(cM, lg(M)-1));
5500 203 : return gerepileupto(av, b);
5501 : }
|
