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,"gauss_pivot. 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 1892 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
124 : {
125 1892 : pari_sp av0 = avma, av, tetpil;
126 : GEN y, c, d;
127 : long i, j, k, r, t, n, m;
128 :
129 1892 : n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
130 1892 : m=nbrows(x); r=0;
131 1892 : x = RgM_shallowcopy(x);
132 1892 : c = zero_zv(m);
133 1892 : d=new_chunk(n+1);
134 1892 : av=avma;
135 6855 : for (k=1; k<=n; k++)
136 : {
137 15162 : for (j=1; j<=m; j++)
138 13010 : if (!c[j])
139 : {
140 9238 : gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
141 9238 : if (!ff->equal0(gcoeff(x,j,k))) break;
142 : }
143 4998 : if (j>m)
144 : {
145 2152 : 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 2117 : r++; d[k]=0;
153 5087 : for(j=1; j<k; j++)
154 2970 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
155 : }
156 : else
157 : {
158 2846 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
159 2846 : c[j] = k; d[k] = j;
160 2846 : gcoeff(x,j,k) = ff->s(E,-1);
161 6853 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
162 16212 : for (t=1; t<=m; t++)
163 : {
164 13366 : if (t==j) continue;
165 :
166 10520 : piv = ff->red(E,gcoeff(x,t,k));
167 10520 : if (ff->equal0(piv)) continue;
168 :
169 3091 : gcoeff(x,t,k) = ff->s(E,0);
170 7549 : for (i=k+1; i<=n; i++)
171 4458 : gcoeff(x,t,i) = ff->red(E, ff->add(E, gcoeff(x,t,i),
172 4458 : ff->mul(E,piv,gcoeff(x,j,i))));
173 3091 : if (gc_needed(av,1))
174 0 : gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
175 : }
176 : }
177 : }
178 1857 : if (deplin) return gc_NULL(av0);
179 :
180 1829 : tetpil=avma; y=cgetg(r+1,t_MAT);
181 3946 : for (j=k=1; j<=r; j++,k++)
182 : {
183 2117 : GEN C = cgetg(n+1,t_COL);
184 2117 : GEN g0 = ff->s(E,0), g1 = ff->s(E,1);
185 4078 : gel(y,j) = C; while (d[k]) k++;
186 5087 : for (i=1; i<k; i++)
187 2970 : if (d[i])
188 : {
189 2420 : GEN p1=gcoeff(x,d[i],k);
190 2420 : gel(C,i) = ff->red(E,p1); gunclone(p1);
191 : }
192 : else
193 550 : gel(C,i) = g0;
194 2962 : gel(C,k) = g1; for (i=k+1; i<=n; i++) gel(C,i) = g0;
195 : }
196 1829 : return gerepile(av0,tetpil,y);
197 : }
198 :
199 : GEN
200 1891 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
201 : {
202 : pari_sp av;
203 : GEN c, d;
204 1891 : long i, j, k, r, t, m, n = lg(x)-1;
205 :
206 1891 : if (!n) { *rr = 0; return NULL; }
207 :
208 1891 : m=nbrows(x); r=0;
209 1891 : d = cgetg(n+1, t_VECSMALL);
210 1891 : x = RgM_shallowcopy(x);
211 1891 : c = zero_zv(m);
212 1891 : av=avma;
213 6830 : for (k=1; k<=n; k++)
214 : {
215 13696 : for (j=1; j<=m; j++)
216 13088 : if (!c[j])
217 : {
218 9200 : gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
219 9200 : if (!ff->equal0(gcoeff(x,j,k))) break;
220 : }
221 4939 : if (j>m) { r++; d[k]=0; }
222 : else
223 : {
224 4331 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
225 4331 : GEN g0 = ff->s(E,0);
226 4331 : c[j] = k; d[k] = j;
227 9032 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
228 26263 : for (t=1; t<=m; t++)
229 : {
230 21932 : if (c[t]) continue; /* already a pivot on that line */
231 :
232 13493 : piv = ff->red(E,gcoeff(x,t,k));
233 13493 : if (ff->equal0(piv)) continue;
234 5346 : gcoeff(x,t,k) = g0;
235 9841 : for (i=k+1; i<=n; i++)
236 4495 : gcoeff(x,t,i) = ff->red(E, ff->add(E,gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i))));
237 5346 : if (gc_needed(av,1))
238 0 : gerepile_gauss(x,k,t,av,j,c);
239 : }
240 13363 : for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
241 : }
242 : }
243 1891 : *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 145618 : _gen_addmul(GEN b, long k, long i, GEN m, void *E, const struct bb_field *ff)
292 : {
293 145618 : gel(b,i) = ff->red(E,gel(b,i));
294 145618 : gel(b,k) = ff->add(E,gel(b,k), ff->mul(E,m, gel(b,i)));
295 145618 : }
296 :
297 : static GEN
298 61717 : _gen_get_col(GEN a, GEN b, long li, void *E, const struct bb_field *ff)
299 : {
300 61717 : GEN u = cgetg(li+1,t_COL);
301 61717 : pari_sp av = avma;
302 : long i, j;
303 :
304 61717 : gel(u,li) = gerepileupto(av, ff->red(E,ff->mul(E,gel(b,li), gcoeff(a,li,li))));
305 318715 : for (i=li-1; i>0; i--)
306 : {
307 256998 : pari_sp av = avma;
308 256998 : GEN m = gel(b,i);
309 1018244 : 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 256998 : m = ff->red(E, m);
311 256998 : gel(u,i) = gerepileupto(av, ff->red(E,ff->mul(E,m, gcoeff(a,i,i))));
312 : }
313 61717 : return u;
314 : }
315 :
316 : GEN
317 13591 : gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff)
318 : {
319 : long i, j, k, li, bco, aco;
320 13591 : GEN u, g0 = ff->s(E,0);
321 13591 : pari_sp av = avma;
322 13591 : a = RgM_shallowcopy(a);
323 13591 : b = RgM_shallowcopy(b);
324 13591 : aco = lg(a)-1; bco = lg(b)-1; li = nbrows(a);
325 60436 : for (i=1; i<=aco; i++)
326 : {
327 : GEN invpiv;
328 73062 : for (k = i; k <= li; k++)
329 : {
330 73020 : GEN piv = ff->red(E,gcoeff(a,k,i));
331 73020 : if (!ff->equal0(piv)) { gcoeff(a,k,i) = ff->inv(E,piv); break; }
332 12626 : gcoeff(a,k,i) = g0;
333 : }
334 : /* found a pivot on line k */
335 60436 : if (k > li) return NULL;
336 60394 : if (k != i)
337 : { /* swap lines so that k = i */
338 52266 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
339 69813 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
340 : }
341 60394 : if (i == aco) break;
342 :
343 46845 : invpiv = gcoeff(a,i,i); /* 1/piv mod p */
344 176142 : for (k=i+1; k<=li; k++)
345 : {
346 129297 : GEN m = ff->red(E,gcoeff(a,k,i)); gcoeff(a,k,i) = g0;
347 129297 : if (ff->equal0(m)) continue;
348 :
349 17297 : m = ff->red(E,ff->neg(E,ff->mul(E,m, invpiv)));
350 71259 : for (j=i+1; j<=aco; j++) _gen_addmul(gel(a,j),k,i,m,E,ff);
351 108953 : for (j=1 ; j<=bco; j++) _gen_addmul(gel(b,j),k,i,m,E,ff);
352 : }
353 46845 : 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 13549 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
361 13549 : u = cgetg(bco+1,t_MAT);
362 75266 : for (j=1; j<=bco; j++) gel(u,j) = _gen_get_col(a, gel(b,j), aco, E, ff);
363 13549 : return u;
364 : }
365 :
366 : /* compatible t_MAT * t_COL, lgA = lg(A) = lg(B) > 1, l = lgcols(A) */
367 : static GEN
368 644407 : gen_matcolmul_i(GEN A, GEN B, ulong lgA, ulong l,
369 : void *E, const struct bb_field *ff)
370 : {
371 644407 : GEN C = cgetg(l, t_COL);
372 : ulong i;
373 4217273 : for (i = 1; i < l; i++) {
374 3572866 : pari_sp av = avma;
375 3572866 : GEN e = ff->mul(E, gcoeff(A, i, 1), gel(B, 1));
376 : ulong k;
377 15412742 : for(k = 2; k < lgA; k++)
378 11839876 : e = ff->add(E, e, ff->mul(E, gcoeff(A, i, k), gel(B, k)));
379 3572866 : gel(C, i) = gerepileupto(av, ff->red(E, e));
380 : }
381 644407 : return C;
382 : }
383 :
384 : GEN
385 206341 : gen_matcolmul(GEN A, GEN B, void *E, const struct bb_field *ff)
386 : {
387 206341 : ulong lgA = lg(A);
388 206341 : if (lgA != (ulong)lg(B))
389 0 : pari_err_OP("operation 'gen_matcolmul'", A, B);
390 206341 : if (lgA == 1)
391 0 : return cgetg(1, t_COL);
392 206341 : return gen_matcolmul_i(A, B, lgA, lgcols(A), E, ff);
393 : }
394 :
395 : static GEN
396 85114 : 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 85114 : GEN C = cgetg(lb, t_MAT);
401 523180 : for(j = 1; j < lb; j++)
402 438066 : gel(C, j) = gen_matcolmul_i(A, gel(B, j), la, l, E, ff);
403 85114 : 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 85114 : gen_matmul_i(GEN A, GEN B, long l, long la, long lb,
583 : void *E, const struct bb_field *ff)
584 : {
585 85114 : if (l <= gen_matmul_sw_bound
586 7 : || la <= gen_matmul_sw_bound
587 0 : || lb <= gen_matmul_sw_bound)
588 85114 : 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 85114 : gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
595 : {
596 85114 : ulong lgA, lgB = lg(B);
597 85114 : if (lgB == 1)
598 0 : return cgetg(1, t_MAT);
599 85114 : lgA = lg(A);
600 85114 : if (lgA != (ulong)lgcols(B))
601 0 : pari_err_OP("operation 'gen_matmul'", A, B);
602 85114 : if (lgA == 1)
603 0 : return zeromat(0, lgB - 1);
604 85114 : return gen_matmul_i(A, B, lgcols(A), lgA, lgB, E, ff);
605 : }
606 :
607 : static GEN
608 19326 : gen_colneg(GEN x, void *E, const struct bb_field *ff)
609 74127 : { pari_APPLY_same(ff->neg(E, gel(x,i))); }
610 :
611 : static GEN
612 4015 : gen_matneg(GEN x, void *E, const struct bb_field *ff)
613 23271 : { pari_APPLY_same(gen_colneg(gel(x,i), E, ff)); }
614 :
615 : static GEN
616 354675 : gen_colscalmul(GEN x, GEN b, void *E, const struct bb_field *ff)
617 826358 : { pari_APPLY_same(ff->red(E, ff->mul(E, gel(x,i), b))); }
618 :
619 : static GEN
620 54434 : gen_matscalmul(GEN x, GEN b, void *E, const struct bb_field *ff)
621 409109 : { pari_APPLY_same(gen_colscalmul(gel(x,i), b, E, ff)); }
622 :
623 : static GEN
624 681846 : gen_colsub(GEN x, GEN y, void *E, const struct bb_field *ff)
625 2432956 : { pari_APPLY_same(ff->add(E, gel(x,i), ff->neg(E, gel(y,i)))); }
626 :
627 : static GEN
628 78576 : gen_matsub(GEN x, GEN y, void *E, const struct bb_field *ff)
629 760422 : { pari_APPLY_same(gen_colsub(gel(x,i), gel(y,i), E, ff)); }
630 :
631 : static GEN
632 14885 : gen_zerocol(long n, void* data, const struct bb_field *R)
633 14885 : { return const_col(n, R->s(data, 0)); }
634 :
635 : static GEN
636 14885 : gen_zeromat(long m, long n, void* data, const struct bb_field *R)
637 : {
638 14885 : GEN M = const_vec(n, gen_zerocol(m, data, R));
639 14885 : 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 gerepilecopy(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 365101 : image_from_pivot(GEN x, GEN d, long r)
738 : {
739 : GEN y;
740 : long j, k;
741 :
742 365101 : if (!d) return gcopy(x);
743 : /* d left on stack for efficiency */
744 360379 : r = lg(x)-1 - r; /* = dim Im(x) */
745 360379 : y = cgetg(r+1,t_MAT);
746 2070499 : for (j=k=1; j<=r; k++)
747 1710118 : if (d[k]) gel(y,j++) = gcopy(gel(x,k));
748 360381 : return y;
749 : }
750 :
751 : /* r = dim Ker x, n = nbrows(x) */
752 : static GEN
753 269719 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
754 : {
755 : pari_sp av;
756 : GEN y, c;
757 269719 : long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
758 :
759 269719 : if (rx == n && r == 0) return gcopy(x);
760 199271 : y = cgetg(n+1, t_MAT);
761 199275 : av = avma; c = zero_zv(n);
762 : /* c = lines containing pivots (could get it from gauss_pivot, but cheap)
763 : * In theory r = 0 and d[j] > 0 for all j, but why take chances? */
764 840742 : for (k = j = 1; j<=rx; j++)
765 641466 : if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
766 1202773 : for (j=1; j<=n; j++)
767 1003497 : if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
768 199276 : set_avma(av);
769 :
770 199276 : rx -= r;
771 840672 : for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
772 561378 : for ( ; j<=n; j++) gel(y,j) = ei(n, y[j]);
773 199276 : return y;
774 : }
775 :
776 : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
777 : static GEN
778 1966868 : indexrank0(long n, long r, GEN d)
779 : {
780 1966868 : GEN p1, p2, res = cgetg(3,t_VEC);
781 : long i, j;
782 :
783 1966865 : r = n - r; /* now r = dim Im(x) */
784 1966865 : p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
785 1966862 : p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
786 1966861 : if (d)
787 : {
788 7896051 : for (i=0,j=1; j<=n; j++)
789 5932698 : if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
790 1963353 : vecsmall_sort(p1);
791 : }
792 1966869 : 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 17592 : indexcompl(GEN v, long n)
820 : {
821 17592 : long i, j, k, m = lg(v) - 1;
822 17592 : GEN w = cgetg(n - m + 1, t_VECSMALL);
823 165740 : for (i = j = k = 1; i <= n; i++)
824 148148 : if (j <= m && v[j] == i) j++; else w[k++] = i;
825 17592 : return w;
826 : }
827 :
828 : static GEN
829 4085 : gen_solve_upper_1(GEN U, GEN B, void *E, const struct bb_field *ff)
830 4085 : { 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 = gerepilecopy(av, X);
872 2658 : return X;
873 : }
874 :
875 : static GEN
876 6180 : gen_lsolve_upper_2(GEN U, GEN B, void *E, const struct bb_field *ff)
877 : {
878 6180 : GEN a = gcoeff(U, 1, 1), b = gcoeff(U, 1, 2), d = gcoeff(U, 2, 2);
879 6180 : GEN D = ff->red(E, ff->mul(E, a, d)), Dinv = ff->inv(E, D);
880 6180 : GEN ainv = ff->red(E, ff->mul(E, d, Dinv)), dinv = ff->red(E, ff->mul(E, a, Dinv));
881 6180 : GEN B1 = vecslice(B, 1, 1);
882 6180 : GEN B2 = vecslice(B, 2, 2);
883 6180 : GEN X1 = gen_matscalmul(B1, ainv, E, ff);
884 6180 : GEN X2 = gen_matscalmul(gen_matsub(B2, gen_matscalmul(X1, b, E, ff), E, ff), dinv, E, ff);
885 6180 : return shallowconcat(X1, X2);
886 : }
887 :
888 : /* solve X*U = B, U upper triangular and invertible */
889 : static GEN
890 14256 : gen_lsolve_upper(GEN U, GEN B, void *E, const struct bb_field *ff,
891 : GEN (*mul)(void *E, GEN a, GEN))
892 : {
893 14256 : long n = lg(U) - 1, n1;
894 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
895 14256 : pari_sp av = avma;
896 :
897 14256 : if (n == 0) return B;
898 14256 : if (n == 1) return gen_solve_upper_1(U, B, E, ff);
899 11097 : if (n == 2) return gen_lsolve_upper_2(U, B, E, ff);
900 4917 : n1 = (n + 1)/2;
901 4917 : U2 = vecslice(U, n1 + 1, n);
902 4917 : U11 = matslice(U, 1,n1, 1,n1);
903 4917 : U12 = rowslice(U2, 1, n1);
904 4917 : U22 = rowslice(U2, n1 + 1, n);
905 4917 : B1 = vecslice(B, 1, n1);
906 4917 : B2 = vecslice(B, n1 + 1, n);
907 4917 : X1 = gen_lsolve_upper(U11, B1, E, ff, mul);
908 4917 : B2 = gen_matsub(B2, mul(E, X1, U12), E, ff);
909 4917 : if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
910 4917 : X2 = gen_lsolve_upper(U22, B2, E, ff, mul);
911 4917 : X = shallowconcat(X1, X2);
912 4917 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
913 4917 : return X;
914 : }
915 :
916 : static GEN
917 17183 : gen_rsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
918 : {
919 17183 : GEN X1 = rowslice(A, 1, 1);
920 17183 : GEN X2 = gen_matsub(rowslice(A, 2, 2), gen_matscalmul(X1, gcoeff(L, 2, 1), E, ff), E, ff);
921 17183 : 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 40198 : 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 40198 : long m = lg(L) - 1, m1, n;
931 : GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
932 40198 : pari_sp av = avma;
933 :
934 40198 : if (m == 0) return zeromat(0, lg(A) - 1);
935 40198 : if (m == 1) return rowslice(A, 1, 1);
936 31768 : if (m == 2) return gen_rsolve_lower_unit_2(L, A, E, ff);
937 14585 : m1 = (m + 1)/2;
938 14585 : n = nbrows(L);
939 14585 : L1 = vecslice(L, 1, m1);
940 14585 : L11 = rowslice(L1, 1, m1);
941 14585 : L21 = rowslice(L1, m1 + 1, n);
942 14585 : A1 = rowslice(A, 1, m1);
943 14585 : X1 = gen_rsolve_lower_unit(L11, A1, E, ff, mul);
944 14585 : A2 = rowslice(A, m1 + 1, n);
945 14585 : A2 = gen_matsub(A2, mul(E, L21, X1), E, ff);
946 14585 : if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
947 14585 : L22 = matslice(L, m1+1,n, m1+1,m);
948 14585 : X2 = gen_rsolve_lower_unit(L22, A2, E, ff, mul);
949 14585 : X = vconcat(X1, X2);
950 14585 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
951 14585 : return X;
952 : }
953 :
954 : static GEN
955 7858 : gen_lsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
956 : {
957 7858 : GEN X2 = vecslice(A, 2, 2);
958 7858 : GEN X1 = gen_matsub(vecslice(A, 1, 1),
959 7858 : gen_matscalmul(X2, gcoeff(L, 2, 1), E, ff), E, ff);
960 7858 : 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 20118 : 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 20118 : long m = lg(L) - 1, m1;
970 : GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
971 20118 : pari_sp av = avma;
972 :
973 20118 : if (m <= 1) return A;
974 15913 : if (m == 2) return gen_lsolve_lower_unit_2(L, A, E, ff);
975 8055 : m1 = (m + 1)/2;
976 8055 : L2 = vecslice(L, m1 + 1, m);
977 8055 : L22 = rowslice(L2, m1 + 1, m);
978 8055 : A2 = vecslice(A, m1 + 1, m);
979 8055 : X2 = gen_lsolve_lower_unit(L22, A2, E, ff, mul);
980 8055 : if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
981 8055 : L1 = vecslice(L, 1, m1);
982 8055 : L21 = rowslice(L1, m1 + 1, m);
983 8055 : A1 = vecslice(A, 1, m1);
984 8055 : A1 = gen_matsub(A1, mul(E, X2, L21), E, ff);
985 8055 : L11 = rowslice(L1, 1, m1);
986 8055 : if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
987 8055 : X1 = gen_lsolve_lower_unit(L11, A1, E, ff, mul);
988 8055 : X = shallowconcat(X1, X2);
989 8055 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
990 8055 : return X;
991 : }
992 :
993 : /* destroy A */
994 : static long
995 22429 : gen_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, void *E, const struct bb_field *ff)
996 : {
997 22429 : long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
998 : pari_sp av;
999 : GEN u, v;
1000 :
1001 22429 : if (P) *P = identity_perm(n);
1002 22429 : *R = cgetg(m + 1, t_VECSMALL);
1003 22429 : av = avma;
1004 58432 : for (j = 1, pr = 0; j <= n; j++)
1005 : {
1006 134204 : for (pr++, pc = 0; pr <= m; pr++)
1007 : {
1008 657259 : for (k = j; k <= n; k++)
1009 : {
1010 540534 : v = ff->red(E, gcoeff(A, pr, k));
1011 540534 : gcoeff(A, pr, k) = v;
1012 540534 : if (!pc && !ff->equal0(v)) pc = k;
1013 : }
1014 116725 : if (pc) break;
1015 : }
1016 53482 : if (!pc) break;
1017 36003 : (*R)[j] = pr;
1018 36003 : if (pc != j)
1019 : {
1020 5073 : swap(gel(A, j), gel(A, pc));
1021 5073 : if (P) lswap((*P)[j], (*P)[pc]);
1022 : }
1023 36003 : u = ff->inv(E, gcoeff(A, pr, j));
1024 179286 : for (i = pr + 1; i <= m; i++)
1025 : {
1026 143283 : v = ff->red(E, ff->mul(E, gcoeff(A, i, j), u));
1027 143283 : gcoeff(A, i, j) = v;
1028 143283 : v = ff->neg(E, v);
1029 462790 : for (k = j + 1; k <= n; k++)
1030 319507 : gcoeff(A, i, k) = ff->add(E, gcoeff(A, i, k),
1031 319507 : ff->red(E, ff->mul(E, gcoeff(A, pr, k), v)));
1032 : }
1033 36003 : if (gc_needed(av, 2)) A = gerepilecopy(av, A);
1034 : }
1035 22429 : setlg(*R, j);
1036 22429 : *C = vecslice(A, 1, j - 1);
1037 22429 : if (U) *U = rowpermute(A, *R);
1038 22429 : return j - 1;
1039 : }
1040 :
1041 : static const long gen_CUP_LIMIT = 5;
1042 :
1043 : static long
1044 11480 : 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 11480 : 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 11480 : pari_sp av = avma;
1051 :
1052 11480 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1053 : /* destroy A; not called at the outermost recursion level */
1054 6573 : return gen_CUP_basecase(A, R, C, U, P, E, ff);
1055 4907 : m1 = (minss(m, n) + 1)/2;
1056 4907 : A1 = rowslice(A, 1, m1);
1057 4907 : A2 = rowslice(A, m1 + 1, m);
1058 4907 : r1 = gen_CUP(A1, &R1, &C1, &U1, &P1, E, ff, mul);
1059 4907 : if (r1 == 0)
1060 : {
1061 485 : r2 = gen_CUP(A2, &R2, &C2, &U2, &P2, E, ff, mul);
1062 485 : *R = cgetg(r2 + 1, t_VECSMALL);
1063 790 : for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
1064 485 : *C = vconcat(gen_zeromat(m1, r2, E, ff), C2);
1065 485 : *U = U2;
1066 485 : *P = P2;
1067 485 : r = r2;
1068 : }
1069 : else
1070 : {
1071 4422 : U11 = vecslice(U1, 1, r1);
1072 4422 : U12 = vecslice(U1, r1 + 1, n);
1073 4422 : T21 = vecslicepermute(A2, P1, 1, r1);
1074 4422 : T22 = vecslicepermute(A2, P1, r1 + 1, n);
1075 4422 : C21 = gen_lsolve_upper(U11, T21, E, ff, mul);
1076 4422 : if (gc_needed(av, 1))
1077 0 : gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
1078 4422 : B2 = gen_matsub(T22, mul(E, C21, U12), E, ff);
1079 4422 : r2 = gen_CUP(B2, &R2, &C2, &U2, &P2, E, ff, mul);
1080 4422 : r = r1 + r2;
1081 4422 : *R = cgetg(r + 1, t_VECSMALL);
1082 19941 : for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
1083 21029 : for ( ; i <= r; i++) (*R)[i] = R2[i - r1] + m1;
1084 4422 : *C = shallowconcat(vconcat(C1, C21),
1085 : vconcat(gen_zeromat(m1, r2, E, ff), C2));
1086 4422 : *U = shallowconcat(vconcat(U11, gen_zeromat(r2, r1, E, ff)),
1087 : vconcat(vecpermute(U12, P2), U2));
1088 :
1089 4422 : *P = cgetg(n + 1, t_VECSMALL);
1090 19941 : for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
1091 50825 : for ( ; i <= n; i++) (*P)[i] = P1[P2[i - r1] + r1];
1092 : }
1093 4907 : if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
1094 4907 : return r;
1095 : }
1096 :
1097 : /* column echelon form */
1098 : static long
1099 27879 : gen_echelon(GEN A, GEN *R, GEN *C, void *E, const struct bb_field *ff,
1100 : GEN (*mul)(void*, GEN, GEN))
1101 : {
1102 27879 : 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 27879 : pari_sp av = avma;
1106 :
1107 27879 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1108 15856 : return gen_CUP_basecase(shallowcopy(A), R, C, NULL, NULL, E, ff);
1109 :
1110 12023 : n1 = (n + 1)/2;
1111 12023 : A1 = vecslice(A, 1, n1);
1112 12023 : A2 = vecslice(A, n1 + 1, n);
1113 12023 : r1 = gen_echelon(A1, &R1, &C1, E, ff, mul);
1114 12023 : if (!r1) return gen_echelon(A2, R, C, E, ff, mul);
1115 10700 : if (r1 == m) { *R = R1; *C = C1; return r1; }
1116 10462 : R1c = indexcompl(R1, m);
1117 10462 : C11 = rowpermute(C1, R1);
1118 10462 : C21 = rowpermute(C1, R1c);
1119 10462 : A12 = rowpermute(A2, R1);
1120 10462 : A22 = rowpermute(A2, R1c);
1121 10462 : M12 = gen_rsolve_lower_unit(C11, A12, E, ff, mul);
1122 10462 : B2 = gen_matsub(A22, mul(E, C21, M12), E, ff);
1123 10462 : r2 = gen_echelon(B2, &R2, &C2, E, ff, mul);
1124 10462 : if (!r2) { *R = R1; *C = C1; r = r1; }
1125 : else
1126 : {
1127 5535 : R2 = perm_mul(R1c, R2);
1128 5535 : C2 = rowpermute(vconcat(gen_zeromat(r1, r2, E, ff), C2),
1129 : perm_inv(vecsmall_concat(R1, R1c)));
1130 5535 : r = r1 + r2;
1131 5535 : *R = cgetg(r + 1, t_VECSMALL);
1132 5535 : *C = cgetg(r + 1, t_MAT);
1133 39631 : for (j = j1 = j2 = 1; j <= r; j++)
1134 34096 : if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
1135 : {
1136 20045 : gel(*C, j) = gel(C1, j1);
1137 20045 : (*R)[j] = R1[j1++];
1138 : }
1139 : else
1140 : {
1141 14051 : gel(*C, j) = gel(C2, j2);
1142 14051 : (*R)[j] = R2[j2++];
1143 : }
1144 : }
1145 10462 : if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
1146 10462 : return r;
1147 : }
1148 :
1149 : static GEN
1150 904 : 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 904 : long i, n = lg(x) - 1, r;
1155 904 : GEN R, C, U, P, d = zero_zv(n);
1156 904 : av = avma;
1157 904 : r = gen_CUP(x, &R, &C, &U, &P, E, ff, mul);
1158 6921 : for(i = 1; i <= r; i++)
1159 6017 : d[P[i]] = R[i];
1160 904 : set_avma(av);
1161 904 : *rr = n - r;
1162 904 : 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 gerepilecopy(av, X);
1217 : }
1218 :
1219 : static GEN
1220 3938 : gen_ker_echelon(GEN x, void *E, const struct bb_field *ff,
1221 : GEN (*mul)(void*, GEN, GEN))
1222 : {
1223 3938 : pari_sp av = avma;
1224 : GEN R, Rc, C, C1, C2, S, K;
1225 3938 : long n = lg(x) - 1, r;
1226 3938 : r = gen_echelon(shallowtrans(x), &R, &C, E, ff, mul);
1227 3938 : Rc = indexcompl(R, n);
1228 3938 : C1 = rowpermute(C, R);
1229 3938 : C2 = rowpermute(C, Rc);
1230 3938 : S = gen_lsolve_lower_unit(C1, C2, E, ff, mul);
1231 3938 : K = vecpermute(shallowconcat(gen_matneg(S, E, ff), gen_matid(n - r, E, ff)),
1232 : perm_inv(vecsmall_concat(R, Rc)));
1233 3938 : K = shallowtrans(K);
1234 3938 : return gerepilecopy(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 gerepilecopy(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 5858 : gen_ker_i(GEN x, long deplin, void *E, const struct bb_field *ff,
1280 : GEN (*mul)(void*, GEN, GEN)) {
1281 5858 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1282 4043 : return deplin? gen_deplin_echelon(x, E, ff, mul): gen_ker_echelon(x, E, ff, mul);
1283 1815 : 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 gerepilecopy(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 2788 : gen_pivots(GEN x, long *rr, void *E, const struct bb_field *ff,
1334 : GEN (*mul)(void*, GEN, GEN))
1335 : {
1336 2788 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1337 904 : return gen_pivots_CUP(x, rr, E, ff, mul);
1338 1884 : 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 gauss_pivot, 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 157256 : FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
1408 : void *E;
1409 : const struct bb_field *ff;
1410 157256 : long n = lg(A) - 1;
1411 :
1412 157256 : if (n == 0)
1413 0 : return cgetg(1, t_MAT);
1414 157256 : if (n > 1)
1415 88394 : return FlxqM_mul_Kronecker(A, B, T, p);
1416 68862 : ff = get_Flxq_field(&E, T, p);
1417 68862 : return gen_matmul(A, B, E, ff);
1418 : }
1419 :
1420 : GEN
1421 86037 : FqM_mul(GEN A, GEN B, GEN T, GEN p) {
1422 : void *E;
1423 86037 : long n = lg(A) - 1;
1424 : const struct bb_field *ff;
1425 86037 : if (n == 0)
1426 0 : return cgetg(1, t_MAT);
1427 86037 : if (n > 1)
1428 81872 : 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 16072 : _FlxqM_mul(void *E, GEN A, GEN B)
1451 : {
1452 16072 : struct _Flxq *D = (struct _Flxq*)E;
1453 16072 : return FlxqM_mul(A, B, D->T, D->p);
1454 : }
1455 :
1456 : static GEN
1457 22741 : _FpM_mul(void *E, GEN A, GEN B)
1458 22741 : { 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 1277129 : FpM_init(GEN a, GEN p, ulong *pp)
1474 : {
1475 1277129 : if (lgefint(p) == 3)
1476 : {
1477 1272843 : *pp = uel(p,2);
1478 1272843 : return (*pp==2)? ZM_to_F2m(a): ZM_to_Flm(a, *pp);
1479 : }
1480 4286 : *pp = 0; return a;
1481 : }
1482 : static GEN
1483 1305477 : FpM_init3(GEN a, GEN p, ulong *pp)
1484 : {
1485 1305477 : if (lgefint(p) == 3)
1486 : {
1487 1302901 : *pp = uel(p,2);
1488 1302901 : switch(*pp)
1489 : {
1490 705938 : case 2: return ZM_to_F2m(a);
1491 156430 : case 3: return ZM_to_F3m(a);
1492 440533 : default:return ZM_to_Flm(a, *pp);
1493 : }
1494 : }
1495 2576 : *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 4676 : FpM_det(GEN a, GEN p)
1532 : {
1533 4676 : pari_sp av = avma;
1534 : ulong pp, d;
1535 4676 : a = FpM_init(a, p, &pp);
1536 4676 : switch(pp)
1537 : {
1538 315 : case 0: return FpM_det_gen(a, p);
1539 1750 : case 2: d = F2m_det_sp(a); break;
1540 2611 : default:d = Flm_det_sp(a,pp); break;
1541 : }
1542 4361 : 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 631626 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
1578 : {
1579 : ulong pp;
1580 631626 : if (lg(x)==1) { *rr = 0; return NULL; }
1581 627555 : x = FpM_init(x, p, &pp);
1582 627554 : switch(pp)
1583 : {
1584 1214 : case 0: return FpM_gauss_pivot_gen(x, p, rr);
1585 350084 : case 2: return F2m_gauss_pivot(x, rr);
1586 276256 : 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 1427 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr) {
1600 : void *E;
1601 1427 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1602 1427 : 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 1504 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
1614 : {
1615 1504 : if (lg(x)==1) { *rr = 0; return NULL; }
1616 1504 : if (!T) return FpM_gauss_pivot(x, p, rr);
1617 1504 : if (lgefint(p) == 3)
1618 : {
1619 1399 : pari_sp av = avma;
1620 1399 : ulong pp = uel(p,2);
1621 1399 : GEN Tp = ZXT_to_FlxT(T, pp);
1622 1399 : GEN d = FlxqM_gauss_pivot(ZXM_to_FlxM(x, pp, get_Flx_var(Tp)), Tp, pp, rr);
1623 1399 : return d ? gerepileuptoleaf(av, d): d;
1624 : }
1625 105 : return FqM_gauss_pivot_gen(x, T, p, rr);
1626 : }
1627 :
1628 : GEN
1629 324146 : FpM_image(GEN x, GEN p)
1630 : {
1631 : long r;
1632 324146 : GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
1633 324144 : 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 297120 : FpM_FpC_invimage(GEN A, GEN x, GEN p)
1774 : {
1775 297120 : pari_sp av = avma;
1776 : ulong pp;
1777 : GEN y;
1778 :
1779 297120 : A = FpM_init(A, p, &pp);
1780 297136 : switch(pp)
1781 : {
1782 8 : case 0: return FpM_FpC_invimage_gen(A, x, p);
1783 193033 : case 2:
1784 193033 : y = F2m_F2c_invimage(A, ZV_to_F2v(x));
1785 193029 : if (!y) return y;
1786 193029 : y = F2c_to_ZC(y);
1787 193029 : return gerepileupto(av, y);
1788 104095 : default:
1789 104095 : y = Flm_Flc_invimage(A, ZV_to_Flv(x, pp), pp);
1790 104094 : if (!y) return y;
1791 104094 : y = Flc_to_ZC(y);
1792 104095 : 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 2646 : FpM_ker_gen(GEN x, GEN p, long deplin)
1819 : {
1820 : void *E;
1821 2646 : const struct bb_field *S = get_Fp_field(&E,p);
1822 2646 : return gen_ker_i(x, deplin, E, S, _FpM_mul);
1823 : }
1824 : static GEN
1825 1305476 : FpM_ker_i(GEN x, GEN p, long deplin)
1826 : {
1827 1305476 : pari_sp av = avma;
1828 : ulong pp;
1829 : GEN y;
1830 :
1831 1305476 : if (lg(x)==1) return cgetg(1,t_MAT);
1832 1305476 : x = FpM_init3(x, p, &pp);
1833 1305499 : switch(pp)
1834 : {
1835 2576 : case 0: return FpM_ker_gen(x,p,deplin);
1836 705955 : case 2:
1837 705955 : y = F2m_ker_sp(x, deplin);
1838 705946 : if (!y) return gc_NULL(av);
1839 705960 : y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
1840 705957 : return gerepileupto(av, y);
1841 156431 : case 3:
1842 156431 : y = F3m_ker_sp(x, deplin);
1843 156431 : if (!y) return gc_NULL(av);
1844 156431 : y = deplin? F3c_to_ZC(y): F3m_to_ZM(y);
1845 156431 : return gerepileupto(av, y);
1846 440537 : default:
1847 440537 : y = Flm_ker_sp(x, pp, deplin);
1848 440539 : if (!y) return gc_NULL(av);
1849 440539 : y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
1850 440539 : return gerepileupto(av, y);
1851 : }
1852 : }
1853 :
1854 : GEN
1855 841985 : 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 2994 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin) {
1876 : void *E;
1877 2994 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1878 2994 : 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 11310 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
1896 : {
1897 11310 : if (!T) return FpM_ker_i(x,p,deplin);
1898 3057 : if (lg(x)==1) return cgetg(1,t_MAT);
1899 :
1900 3057 : if (lgefint(p)==3)
1901 : {
1902 2931 : pari_sp av = avma;
1903 2931 : ulong l = p[2];
1904 2931 : GEN Tl = ZXT_to_FlxT(T,l);
1905 2931 : GEN Ml = ZXM_to_FlxM(x, l, get_Flx_var(Tl));
1906 2931 : GEN K = FlxqM_ker_i(Ml, Tl, l, deplin);
1907 2931 : 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 2924 : return gerepileupto(av, K);
1911 : }
1912 126 : return FqM_ker_gen(x, T, p, deplin);
1913 : }
1914 :
1915 : GEN
1916 11226 : FqM_ker(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,0); }
1917 :
1918 : GEN
1919 455251 : 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 347800 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
1946 : {
1947 347800 : long n = nbrows(a);
1948 347799 : a = FpM_init(a,p,pp);
1949 347798 : switch(*pp)
1950 : {
1951 2749 : case 0:
1952 2749 : if (!b) b = matid(n);
1953 2749 : return FpM_gauss_gen(a,b,p);
1954 226963 : case 2:
1955 226963 : if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
1956 226961 : return F2m_gauss_sp(a,b);
1957 118086 : default:
1958 118086 : if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
1959 118086 : 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 gerepilecopy(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) { set_avma(av); return cgetg(1, t_MAT); }
1995 14 : u = F2xqM_gauss_gen(a, b, T);
1996 14 : if (!u) return gc_NULL(av);
1997 14 : return gerepilecopy(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) { set_avma(av); return 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 gerepilecopy(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 gerepilecopy(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 gerepilecopy(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 gerepilecopy(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 gerepilecopy(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 gerepilecopy(av, gel(u,1));
2089 : }
2090 :
2091 : GEN
2092 347751 : FpM_inv(GEN a, GEN p)
2093 : {
2094 347751 : pari_sp av = avma;
2095 : ulong pp;
2096 : GEN u;
2097 347751 : if (lg(a) == 1) return cgetg(1, t_MAT);
2098 347751 : u = FpM_gauss_i(a, NULL, p, &pp);
2099 347746 : if (!u) return gc_NULL(av);
2100 347732 : switch(pp)
2101 : {
2102 2693 : case 0: return gerepilecopy(av, u);
2103 226953 : case 2: u = F2m_to_ZM(u); break;
2104 118086 : default: u = Flm_to_ZM(u); break;
2105 : }
2106 345039 : 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) { set_avma(av); return 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 gerepilecopy(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) { set_avma(av); return 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 gerepilecopy(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 gerepilecopy(av, u);
2141 : }
2142 :
2143 : GEN
2144 262923 : FpM_intersect_i(GEN x, GEN y, GEN p)
2145 : {
2146 262923 : long j, lx = lg(x);
2147 : GEN z;
2148 :
2149 262923 : if (lx == 1 || lg(y) == 1) return cgetg(1,t_MAT);
2150 262923 : if (lgefint(p) == 3)
2151 : {
2152 262923 : ulong pp = p[2];
2153 262923 : 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 : static void
2175 269722 : init_suppl(GEN x)
2176 : {
2177 269722 : if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
2178 : /* HACK: avoid overwriting d from gauss_pivot() after set_avma(av) */
2179 269722 : (void)new_chunk(lgcols(x) * 2);
2180 269720 : }
2181 :
2182 : GEN
2183 268154 : FpM_suppl(GEN x, GEN p)
2184 : {
2185 : GEN d;
2186 : long r;
2187 268154 : init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
2188 268153 : return get_suppl(x,d,nbrows(x),r,&col_ei);
2189 : }
2190 :
2191 : GEN
2192 14 : F2m_suppl(GEN x)
2193 : {
2194 : GEN d;
2195 : long r;
2196 14 : init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
2197 14 : return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
2198 : }
2199 :
2200 : GEN
2201 105 : Flm_suppl(GEN x, ulong p)
2202 : {
2203 : GEN d;
2204 : long r;
2205 105 : init_suppl(x); d = Flm_pivots(x, p, &r, 0);
2206 105 : return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
2207 : }
2208 :
2209 : GEN
2210 7 : F2xqM_suppl(GEN x, GEN T)
2211 : {
2212 : void *E;
2213 7 : const struct bb_field *S = get_F2xq_field(&E, T);
2214 7 : return gen_suppl(x, E, S, _F2xqM_mul);
2215 : }
2216 :
2217 : GEN
2218 14 : FlxqM_suppl(GEN x, GEN T, ulong p)
2219 : {
2220 : void *E;
2221 14 : const struct bb_field *S = get_Flxq_field(&E, T, p);
2222 14 : return gen_suppl(x, E, S, _FlxqM_mul);
2223 : }
2224 :
2225 : GEN
2226 5319 : FqM_suppl(GEN x, GEN T, GEN p)
2227 : {
2228 5319 : pari_sp av = avma;
2229 : GEN d;
2230 : long r;
2231 :
2232 5319 : if (!T) return FpM_suppl(x,p);
2233 1378 : init_suppl(x);
2234 1378 : d = FqM_gauss_pivot(x,T,p,&r);
2235 1378 : set_avma(av); return get_suppl(x,d,nbrows(x),r,&col_ei);
2236 : }
2237 :
2238 : static void
2239 1951258 : init_indexrank(GEN x) {
2240 1951258 : (void)new_chunk(3 + 2*lg(x)); /* HACK */
2241 1951257 : }
2242 :
2243 : GEN
2244 39275 : FpM_indexrank(GEN x, GEN p) {
2245 39275 : pari_sp av = avma;
2246 : long r;
2247 : GEN d;
2248 39275 : init_indexrank(x);
2249 39275 : d = FpM_gauss_pivot(x,p,&r);
2250 39275 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2251 : }
2252 :
2253 : GEN
2254 56448 : Flm_indexrank(GEN x, ulong p) {
2255 56448 : pari_sp av = avma;
2256 : long r;
2257 : GEN d;
2258 56448 : init_indexrank(x);
2259 56448 : d = Flm_pivots(x, p, &r, 0);
2260 56448 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2261 : }
2262 :
2263 : GEN
2264 53 : F2m_indexrank(GEN x) {
2265 53 : pari_sp av = avma;
2266 : long r;
2267 : GEN d;
2268 53 : init_indexrank(x);
2269 53 : d = F2m_gauss_pivot(F2m_copy(x),&r);
2270 53 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2271 : }
2272 :
2273 : GEN
2274 7 : F2xqM_indexrank(GEN x, GEN T) {
2275 7 : pari_sp av = avma;
2276 : long r;
2277 : GEN d;
2278 7 : init_indexrank(x);
2279 7 : d = F2xqM_gauss_pivot(x, T, &r);
2280 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2281 : }
2282 :
2283 : GEN
2284 7 : FlxqM_indexrank(GEN x, GEN T, ulong p) {
2285 7 : pari_sp av = avma;
2286 : long r;
2287 : GEN d;
2288 7 : init_indexrank(x);
2289 7 : d = FlxqM_gauss_pivot(x, T, p, &r);
2290 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2291 : }
2292 :
2293 : GEN
2294 7 : FqM_indexrank(GEN x, GEN T, GEN p) {
2295 7 : pari_sp av = avma;
2296 : long r;
2297 : GEN d;
2298 7 : init_indexrank(x);
2299 7 : d = FqM_gauss_pivot(x, T, p, &r);
2300 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2301 : }
2302 :
2303 : /*******************************************************************/
2304 : /* */
2305 : /* Solve A*X=B (Gauss pivot) */
2306 : /* */
2307 : /*******************************************************************/
2308 : /* x a column, x0 same column in the original input matrix (for reference),
2309 : * c list of pivots so far */
2310 : static long
2311 2592459 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
2312 : {
2313 2592459 : GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
2314 2592459 : long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
2315 2592459 : if (c)
2316 : {
2317 580234 : for (i=1; i<lx; i++)
2318 361106 : if (!c[i])
2319 : {
2320 146270 : long e = gexpo(gel(x,i));
2321 146270 : if (e > ex) { ex = e; k = i; }
2322 : }
2323 : }
2324 : else
2325 : {
2326 8418557 : for (i=ix; i<lx; i++)
2327 : {
2328 6045218 : long e = gexpo(gel(x,i));
2329 6045226 : if (e > ex) { ex = e; k = i; }
2330 : }
2331 : }
2332 2592467 : if (!k) return lx;
2333 2477484 : p = gel(x,k);
2334 2477484 : r = gel(x0,k); if (isrationalzero(r)) r = x0;
2335 2477496 : return cx_approx0(p, r)? lx: k;
2336 : }
2337 : static long
2338 201792 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
2339 : {
2340 201792 : GEN x = gel(X, ix);
2341 201792 : long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
2342 201792 : if (c)
2343 : {
2344 504 : for (i=1; i<lx; i++)
2345 378 : if (!c[i] && !gequal0(gel(x,i)))
2346 : {
2347 245 : long e = gvaluation(gel(x,i), p);
2348 245 : if (e < ex) { ex = e; k = i; }
2349 : }
2350 : }
2351 : else
2352 : {
2353 1721226 : for (i=ix; i<lx; i++)
2354 1519560 : if (!gequal0(gel(x,i)))
2355 : {
2356 1146977 : long e = gvaluation(gel(x,i), p);
2357 1146977 : if (e < ex) { ex = e; k = i; }
2358 : }
2359 : }
2360 201792 : return k? k: lx;
2361 : }
2362 : static long
2363 4501 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
2364 : {
2365 4501 : GEN x = gel(X, ix);
2366 4501 : long i, lx = lg(x);
2367 : (void)x0;
2368 4501 : if (c)
2369 : {
2370 12775 : for (i=1; i<lx; i++)
2371 11795 : if (!c[i] && !gequal0(gel(x,i))) return i;
2372 : }
2373 : else
2374 : {
2375 2380 : for (i=ix; i<lx; i++)
2376 2366 : if (!gequal0(gel(x,i))) return i;
2377 : }
2378 994 : return lx;
2379 : }
2380 :
2381 : /* Return pivot seeking function appropriate for the domain of the RgM x
2382 : * (first non zero pivot, maximal pivot...)
2383 : * x0 is a reference point used when guessing whether x[i,j] ~ 0
2384 : * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
2385 : * but use original x when deciding whether a prospective pivot is nonzero */
2386 : static pivot_fun
2387 801731 : get_pivot_fun(GEN x, GEN x0, GEN *data)
2388 : {
2389 801731 : long i, j, hx, lx = lg(x);
2390 801731 : int res = t_INT;
2391 801731 : GEN p = NULL;
2392 :
2393 801731 : *data = NULL;
2394 801731 : if (lx == 1) return &gauss_get_pivot_NZ;
2395 801696 : hx = lgcols(x);
2396 3659768 : for (j=1; j<lx; j++)
2397 : {
2398 2858115 : GEN xj = gel(x,j);
2399 15866560 : for (i=1; i<hx; i++)
2400 : {
2401 13008487 : GEN c = gel(xj,i);
2402 13008487 : switch(typ(c))
2403 : {
2404 7423929 : case t_REAL:
2405 7423929 : res = t_REAL;
2406 7423929 : break;
2407 3640 : case t_COMPLEX:
2408 3640 : if (typ(gel(c,1)) == t_REAL || typ(gel(c,2)) == t_REAL) res = t_REAL;
2409 3640 : break;
2410 4399018 : case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT: case t_QUAD:
2411 : case t_POLMOD: /* exact types */
2412 4399018 : break;
2413 1181858 : case t_PADIC:
2414 1181858 : p = gel(c,2);
2415 1181858 : res = t_PADIC;
2416 1181858 : break;
2417 42 : default: return &gauss_get_pivot_NZ;
2418 : }
2419 : }
2420 : }
2421 801653 : switch(res)
2422 : {
2423 773319 : case t_REAL: *data = x0; return &gauss_get_pivot_max;
2424 26991 : case t_PADIC: *data = p; return &gauss_get_pivot_padic;
2425 1343 : default: return &gauss_get_pivot_NZ;
2426 : }
2427 : }
2428 :
2429 : static GEN
2430 1265554 : get_col(GEN a, GEN b, GEN p, long li)
2431 : {
2432 1265554 : GEN u = cgetg(li+1,t_COL);
2433 : long i, j;
2434 :
2435 1265556 : gel(u,li) = gdiv(gel(b,li), p);
2436 5149499 : for (i=li-1; i>0; i--)
2437 : {
2438 3883956 : pari_sp av = avma;
2439 3883956 : GEN m = gel(b,i);
2440 17086808 : for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
2441 3883924 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
2442 : }
2443 1265543 : return u;
2444 : }
2445 :
2446 : /* bk -= m * bi */
2447 : static void
2448 18289033 : _submul(GEN b, long k, long i, GEN m)
2449 : {
2450 18289033 : gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
2451 18288945 : }
2452 : static int
2453 2408669 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
2454 : {
2455 2408669 : *iscol = *b ? (typ(*b) == t_COL): 0;
2456 2408669 : *aco = lg(a) - 1;
2457 2408669 : if (!*aco) /* a empty */
2458 : {
2459 70 : if (*b && lg(*b) != 1) pari_err_DIM("gauss");
2460 70 : *li = 0; return 0;
2461 : }
2462 2408599 : *li = nbrows(a);
2463 2408600 : if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
2464 2408601 : if (*b)
2465 : {
2466 2144345 : switch(typ(*b))
2467 : {
2468 121613 : case t_MAT:
2469 121613 : if (lg(*b) == 1) return 0;
2470 121613 : *b = RgM_shallowcopy(*b);
2471 121613 : break;
2472 2022732 : case t_COL:
2473 2022732 : *b = mkmat( leafcopy(*b) );
2474 2022733 : break;
2475 0 : default: pari_err_TYPE("gauss",*b);
2476 : }
2477 2144346 : if (nbrows(*b) != *li) pari_err_DIM("gauss");
2478 : }
2479 : else
2480 264256 : *b = matid(*li);
2481 2408601 : return 1;
2482 : }
2483 :
2484 : static GEN
2485 2051 : RgM_inv_FpM(GEN a, GEN p)
2486 : {
2487 : ulong pp;
2488 2051 : a = RgM_Fp_init(a, p, &pp);
2489 2051 : switch(pp)
2490 : {
2491 35 : case 0:
2492 35 : a = FpM_inv(a,p);
2493 35 : if (a) a = FpM_to_mod(a, p);
2494 35 : break;
2495 189 : case 2:
2496 189 : a = F2m_inv(a);
2497 189 : if (a) a = F2m_to_mod(a);
2498 189 : break;
2499 1827 : default:
2500 1827 : a = Flm_inv_sp(a, NULL, pp);
2501 1827 : if (a) a = Flm_to_mod(a, pp);
2502 : }
2503 2051 : return a;
2504 : }
2505 :
2506 : static GEN
2507 42 : RgM_inv_FqM(GEN x, GEN pol, GEN p)
2508 : {
2509 42 : pari_sp av = avma;
2510 42 : GEN b, T = RgX_to_FpX(pol, p);
2511 42 : if (signe(T) == 0) pari_err_OP("^",x,gen_m1);
2512 42 : b = FqM_inv(RgM_to_FqM(x, T, p), T, p);
2513 42 : if (!b) return gc_NULL(av);
2514 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
2515 : }
2516 :
2517 : #define code(t1,t2) ((t1 << 6) | t2)
2518 : static GEN
2519 526944 : RgM_inv_fast(GEN x)
2520 : {
2521 : GEN p, pol;
2522 : long pa;
2523 526944 : long t = RgM_type(x, &p,&pol,&pa);
2524 526948 : switch(t)
2525 : {
2526 48314 : case t_INT: /* Fall back */
2527 48314 : case t_FRAC: return QM_inv(x);
2528 147 : case t_FFELT: return FFM_inv(x, pol);
2529 2051 : case t_INTMOD: return RgM_inv_FpM(x, p);
2530 42 : case code(t_POLMOD, t_INTMOD):
2531 42 : return RgM_inv_FqM(x, pol, p);
2532 476394 : default: return gen_0;
2533 : }
2534 : }
2535 : #undef code
2536 :
2537 : static GEN
2538 63 : RgM_RgC_solve_FpC(GEN a, GEN b, GEN p)
2539 : {
2540 63 : pari_sp av = avma;
2541 : ulong pp;
2542 63 : a = RgM_Fp_init(a, p, &pp);
2543 63 : switch(pp)
2544 : {
2545 14 : case 0:
2546 14 : b = RgC_to_FpC(b, p);
2547 14 : a = FpM_FpC_gauss(a,b,p);
2548 14 : return a ? gerepileupto(av, FpC_to_mod(a, p)): NULL;
2549 28 : case 2:
2550 28 : b = RgV_to_F2v(b);
2551 28 : a = F2m_F2c_gauss(a,b);
2552 28 : return a ? gerepileupto(av, F2c_to_mod(a)): NULL;
2553 21 : default:
2554 21 : b = RgV_to_Flv(b, pp);
2555 21 : a = Flm_Flc_gauss(a, b, pp);
2556 21 : return a ? gerepileupto(av, Flc_to_mod(a, pp)): NULL;
2557 : }
2558 : }
2559 :
2560 : static GEN
2561 105 : RgM_solve_FpM(GEN a, GEN b, GEN p)
2562 : {
2563 105 : pari_sp av = avma;
2564 : ulong pp;
2565 105 : a = RgM_Fp_init(a, p, &pp);
2566 105 : switch(pp)
2567 : {
2568 35 : case 0:
2569 35 : b = RgM_to_FpM(b, p);
2570 35 : a = FpM_gauss(a,b,p);
2571 35 : return a ? gerepileupto(av, FpM_to_mod(a, p)): NULL;
2572 28 : case 2:
2573 28 : b = RgM_to_F2m(b);
2574 28 : a = F2m_gauss(a,b);
2575 28 : return a ? gerepileupto(av, F2m_to_mod(a)): NULL;
2576 42 : default:
2577 42 : b = RgM_to_Flm(b, pp);
2578 42 : a = Flm_gauss(a,b,pp);
2579 42 : return a ? gerepileupto(av, Flm_to_mod(a, pp)): NULL;
2580 : }
2581 : }
2582 :
2583 : /* Gaussan Elimination. If a is square, return a^(-1)*b;
2584 : * if a has more rows than columns and b is NULL, return c such that c a = Id.
2585 : * a is a (not necessarily square) matrix
2586 : * b is a matrix or column vector, NULL meaning: take the identity matrix,
2587 : * effectively returning the inverse of a
2588 : * If a and b are empty, the result is the empty matrix.
2589 : *
2590 : * li: number of rows of a and b
2591 : * aco: number of columns of a
2592 : * bco: number of columns of b (if matrix)
2593 : */
2594 : static GEN
2595 1693154 : RgM_solve_basecase(GEN a, GEN b)
2596 : {
2597 1693154 : pari_sp av = avma;
2598 : long i, j, k, li, bco, aco;
2599 : int iscol;
2600 : pivot_fun pivot;
2601 : GEN p, u, data;
2602 :
2603 1693154 : set_avma(av);
2604 :
2605 1693153 : if (lg(a)-1 == 2 && nbrows(a) == 2) {
2606 : /* 2x2 matrix, start by inverting a */
2607 1027057 : GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
2608 1027057 : GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
2609 1027057 : GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
2610 1027056 : if (gequal0(D)) return NULL;
2611 1027054 : ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
2612 1027052 : ainv = gmul(ainv, ginv(D));
2613 1027044 : if (b) ainv = gmul(ainv, b);
2614 1027045 : return gerepileupto(av, ainv);
2615 : }
2616 :
2617 666096 : if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
2618 666097 : pivot = get_pivot_fun(a, a, &data);
2619 666097 : a = RgM_shallowcopy(a);
2620 666096 : bco = lg(b)-1;
2621 666096 : if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
2622 :
2623 666096 : p = NULL; /* gcc -Wall */
2624 2292190 : for (i=1; i<=aco; i++)
2625 : {
2626 : /* k is the line where we find the pivot */
2627 2292183 : k = pivot(a, data, i, NULL);
2628 2292194 : if (k > li) return NULL;
2629 2292179 : if (k != i)
2630 : { /* exchange the lines s.t. k = i */
2631 1795602 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
2632 1739854 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
2633 : }
2634 2292179 : p = gcoeff(a,i,i);
2635 2292179 : if (i == aco) break;
2636 :
2637 5118268 : for (k=i+1; k<=li; k++)
2638 : {
2639 3492195 : GEN m = gcoeff(a,k,i);
2640 3492195 : if (!gequal0(m))
2641 : {
2642 2837097 : m = gdiv(m,p);
2643 12111213 : for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
2644 11852294 : for (j=1; j<=bco; j++) _submul(gel(b,j),k,i,m);
2645 : }
2646 : }
2647 1626073 : if (gc_needed(av,1))
2648 : {
2649 12 : if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
2650 12 : gerepileall(av,2, &a,&b);
2651 : }
2652 : }
2653 :
2654 666087 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
2655 666087 : u = cgetg(bco+1,t_MAT);
2656 1931623 : for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
2657 666070 : return gerepilecopy(av, iscol? gel(u,1): u);
2658 : }
2659 :
2660 : static GEN
2661 1183437 : RgM_RgC_solve_fast(GEN x, GEN y)
2662 : {
2663 : GEN p, pol;
2664 : long pa;
2665 1183437 : long t = RgM_RgC_type(x, y, &p,&pol,&pa);
2666 1183437 : switch(t)
2667 : {
2668 14930 : case t_INT: return ZM_gauss(x, y);
2669 175 : case t_FRAC: return QM_gauss(x, y);
2670 63 : case t_INTMOD: return RgM_RgC_solve_FpC(x, y, p);
2671 42 : case t_FFELT: return FFM_FFC_gauss(x, y, pol);
2672 1168227 : default: return gen_0;
2673 : }
2674 : }
2675 :
2676 : static GEN
2677 48790 : RgM_solve_fast(GEN x, GEN y)
2678 : {
2679 : GEN p, pol;
2680 : long pa;
2681 48790 : long t = RgM_type2(x, y, &p,&pol,&pa);
2682 48790 : switch(t)
2683 : {
2684 77 : case t_INT: return ZM_gauss(x, y);
2685 14 : case t_FRAC: return QM_gauss(x, y);
2686 105 : case t_INTMOD: return RgM_solve_FpM(x, y, p);
2687 56 : case t_FFELT: return FFM_gauss(x, y, pol);
2688 48538 : default: return gen_0;
2689 : }
2690 : }
2691 :
2692 : GEN
2693 1232227 : RgM_solve(GEN a, GEN b)
2694 : {
2695 1232227 : pari_sp av = avma;
2696 : GEN u;
2697 1232227 : if (!b) return RgM_inv(a);
2698 1232227 : u = typ(b)==t_MAT ? RgM_solve_fast(a, b): RgM_RgC_solve_fast(a, b);
2699 1232227 : if (!u) return gc_NULL(av);
2700 1232122 : if (u != gen_0) return u;
2701 1216765 : return RgM_solve_basecase(a, b);
2702 : }
2703 :
2704 : GEN
2705 28 : RgM_div(GEN a, GEN b)
2706 : {
2707 28 : pari_sp av = avma;
2708 28 : GEN u = RgM_solve(shallowtrans(b), shallowtrans(a));
2709 28 : if (!u) return gc_NULL(av);
2710 21 : return gerepilecopy(av, shallowtrans(u));
2711 : }
2712 :
2713 : GEN
2714 526944 : RgM_inv(GEN a)
2715 : {
2716 526944 : GEN b = RgM_inv_fast(a);
2717 526934 : return b==gen_0? RgM_solve_basecase(a, NULL): b;
2718 : }
2719 :
2720 : /* assume dim A >= 1, A invertible + upper triangular */
2721 : static GEN
2722 3229723 : RgM_inv_upper_ind(GEN A, long index)
2723 : {
2724 3229723 : long n = lg(A)-1, i = index, j;
2725 3229723 : GEN u = zerocol(n);
2726 3229714 : gel(u,i) = ginv(gcoeff(A,i,i));
2727 6532867 : for (i--; i>0; i--)
2728 : {
2729 3303146 : pari_sp av = avma;
2730 3303146 : GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
2731 14638164 : for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
2732 3303125 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
2733 : }
2734 3229721 : return u;
2735 : }
2736 : GEN
2737 1617027 : RgM_inv_upper(GEN A)
2738 : {
2739 : long i, l;
2740 1617027 : GEN B = cgetg_copy(A, &l);
2741 4846740 : for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
2742 1617017 : return B;
2743 : }
2744 :
2745 : static GEN
2746 4517678 : split_realimag_col(GEN z, long r1, long r2)
2747 : {
2748 4517678 : long i, ru = r1+r2;
2749 4517678 : GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
2750 12541072 : for (i=1; i<=r1; i++) {
2751 8023394 : GEN a = gel(z,i);
2752 8023394 : if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
2753 8023394 : gel(x,i) = a;
2754 : }
2755 7225327 : for ( ; i<=ru; i++) {
2756 2707649 : GEN b, a = gel(z,i);
2757 2707649 : if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
2758 2707649 : gel(x,i) = a;
2759 2707649 : gel(y,i) = b;
2760 : }
2761 4517678 : return x;
2762 : }
2763 : GEN
2764 2570416 : split_realimag(GEN x, long r1, long r2)
2765 : {
2766 2570416 : if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
2767 4503457 : pari_APPLY_same(split_realimag_col(gel(x,i), r1, r2));
2768 : }
2769 :
2770 : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
2771 : * r1 first lines of M,y are real. Solve the system obtained by splitting
2772 : * real and imaginary parts. */
2773 : GEN
2774 1215627 : RgM_solve_realimag(GEN M, GEN y)
2775 : {
2776 1215627 : long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
2777 1215627 : return RgM_solve(split_realimag(M, r1,r2),
2778 : split_realimag(y, r1,r2));
2779 : }
2780 :
2781 : GEN
2782 434 : gauss(GEN a, GEN b)
2783 : {
2784 : GEN z;
2785 434 : long t = typ(b);
2786 434 : if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
2787 434 : if (t!=t_COL && t!=t_MAT) pari_err_TYPE("gauss",b);
2788 434 : z = RgM_solve(a,b);
2789 434 : if (!z) pari_err_INV("gauss",a);
2790 329 : return z;
2791 : }
2792 :
2793 : /* #C = n, C[z[i]] = K[i], complete by 0s */
2794 : static GEN
2795 14 : RgC_inflate(GEN K, GEN z, long n)
2796 : {
2797 14 : GEN c = zerocol(n);
2798 14 : long j, l = lg(K);
2799 28 : for (j = 1; j < l; j++) gel(c, z[j]) = gel(K, j);
2800 14 : return c;
2801 : }
2802 : /* in place: C[i] *= cB / v[i] */
2803 : static void
2804 6482 : QC_normalize(GEN C, GEN v, GEN cB)
2805 : {
2806 6482 : long l = lg(C), i;
2807 48531 : for (i = 1; i < l; i++)
2808 : {
2809 42049 : GEN c = cB, k = gel(C,i), d = gel(v,i);
2810 42049 : if (d)
2811 : {
2812 24958 : if (isintzero(d)) { gel(C,i) = gen_0; continue; }
2813 24958 : c = div_content(c, d);
2814 : }
2815 42049 : gel(C,i) = c? gmul(k,c): k;
2816 : }
2817 6482 : }
2818 :
2819 : /* same as above, M rational; if flag = 1, call indexrank and return 1 sol */
2820 : GEN
2821 6475 : QM_gauss_i(GEN M, GEN B, long flag)
2822 : {
2823 6475 : pari_sp av = avma;
2824 : long i, l, n;
2825 6475 : int col = typ(B) == t_COL;
2826 6475 : GEN K, cB, N = cgetg_copy(M, &l), v = cgetg(l, t_VEC), z2 = NULL;
2827 :
2828 48552 : for (i = 1; i < l; i++)
2829 42077 : gel(N,i) = Q_primitive_part(gel(M,i), &gel(v,i));
2830 6475 : if (flag)
2831 : {
2832 329 : GEN z = ZM_indexrank(N), z1 = gel(z,1);
2833 329 : z2 = gel(z,2);
2834 329 : N = shallowmatextract(N, z1, z2);
2835 329 : B = col? vecpermute(B,z1): rowpermute(B,z1);
2836 329 : if (lg(z2) == l) z2 = NULL; else v = vecpermute(v, z2);
2837 : }
2838 6475 : B = Q_primitive_part(B, &cB);
2839 6475 : K = ZM_gauss(N, B); if (!K) return gc_NULL(av);
2840 6475 : n = l - 1;
2841 6475 : if (col)
2842 : {
2843 6447 : QC_normalize(K, v, cB);
2844 6447 : if (z2) K = RgC_inflate(K, z2, n);
2845 : }
2846 : else
2847 : {
2848 28 : long lK = lg(K);
2849 63 : for (i = 1; i < lK; i++)
2850 : {
2851 35 : QC_normalize(gel(K,i), v, cB);
2852 35 : if (z2) gel(K,i) = RgC_inflate(gel(K,i), z2, n);
2853 : }
2854 : }
2855 6475 : return gerepilecopy(av, K);
2856 : }
2857 : GEN
2858 6146 : QM_gauss(GEN M, GEN B) { return QM_gauss_i(M, B, 0); }
2859 :
2860 : static GEN
2861 790318 : ZM_inv_slice(GEN A, GEN P, GEN *mod)
2862 : {
2863 790318 : pari_sp av = avma;
2864 790318 : long i, n = lg(P)-1;
2865 : GEN H, T;
2866 790318 : if (n == 1)
2867 : {
2868 757852 : ulong p = uel(P,1);
2869 757852 : GEN Hp, a = ZM_to_Flm(A, p);
2870 757854 : Hp = Flm_adjoint(a, p);
2871 757854 : Hp = gerepileupto(av, Flm_to_ZM(Hp));
2872 757852 : *mod = utoipos(p); return Hp;
2873 : }
2874 32466 : T = ZV_producttree(P);
2875 32467 : A = ZM_nv_mod_tree(A, P, T);
2876 32467 : H = cgetg(n+1, t_VEC);
2877 180538 : for(i=1; i <= n; i++)
2878 148071 : gel(H,i) = Flm_adjoint(gel(A, i), uel(P,i));
2879 32467 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
2880 32467 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2881 : }
2882 :
2883 : static GEN
2884 717455 : RgM_true_Hadamard(GEN a)
2885 : {
2886 717455 : pari_sp av = avma;
2887 717455 : long n = lg(a)-1, i;
2888 : GEN B;
2889 717455 : if (n == 0) return gen_1;
2890 717455 : a = RgM_gtofp(a, LOWDEFAULTPREC);
2891 717455 : B = gnorml2(gel(a,1));
2892 2956277 : for (i = 2; i <= n; i++) B = gmul(B, gnorml2(gel(a,i)));
2893 717454 : return gerepileuptoint(av, ceil_safe(sqrtr(B)));
2894 : }
2895 :
2896 : GEN
2897 790318 : ZM_inv_worker(GEN P, GEN A)
2898 : {
2899 790318 : GEN V = cgetg(3, t_VEC);
2900 790318 : gel(V,1) = ZM_inv_slice(A, P, &gel(V,2));
2901 790319 : return V;
2902 : }
2903 :
2904 : static GEN
2905 43504 : ZM_inv0(GEN A, GEN *pden)
2906 : {
2907 43504 : if (pden) *pden = gen_1;
2908 43504 : (void)A; return cgetg(1, t_MAT);
2909 : }
2910 : static GEN
2911 644112 : ZM_inv1(GEN A, GEN *pden)
2912 : {
2913 644112 : GEN a = gcoeff(A,1,1);
2914 644112 : long s = signe(a);
2915 644112 : if (!s) return NULL;
2916 644112 : if (pden) *pden = absi(a);
2917 644112 : retmkmat(mkcol(s == 1? gen_1: gen_m1));
2918 : }
2919 : static GEN
2920 724879 : ZM_inv2(GEN A, GEN *pden)
2921 : {
2922 : GEN a, b, c, d, D, cA;
2923 : long s;
2924 724879 : A = Q_primitive_part(A, &cA);
2925 724877 : a = gcoeff(A,1,1); b = gcoeff(A,1,2);
2926 724877 : c = gcoeff(A,2,1); d = gcoeff(A,2,2);
2927 724877 : D = subii(mulii(a,d), mulii(b,c)); /* left on stack */
2928 724874 : s = signe(D);
2929 724874 : if (!s) return NULL;
2930 724860 : if (s < 0) D = negi(D);
2931 724861 : if (pden) *pden = mul_denom(D, cA);
2932 724861 : if (s > 0)
2933 683363 : retmkmat2(mkcol2(icopy(d), negi(c)), mkcol2(negi(b), icopy(a)));
2934 : else
2935 41498 : retmkmat2(mkcol2(negi(d), icopy(c)), mkcol2(icopy(b), negi(a)));
2936 : }
2937 :
2938 : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
2939 : * not available. Return H primitive such that M*H = den*Id */
2940 : GEN
2941 0 : ZM_inv_ratlift(GEN M, GEN *pden)
2942 : {
2943 0 : pari_sp av2, av = avma;
2944 : GEN Hp, q, H;
2945 : ulong p;
2946 0 : long m = lg(M)-1;
2947 : forprime_t S;
2948 : pari_timer ti;
2949 :
2950 0 : if (m == 0) return ZM_inv0(M,pden);
2951 0 : if (m == 1 && nbrows(M)==1) return ZM_inv1(M,pden);
2952 0 : if (m == 2 && nbrows(M)==2) return ZM_inv2(M,pden);
2953 :
2954 0 : if (DEBUGLEVEL>5) timer_start(&ti);
2955 0 : init_modular_big(&S);
2956 0 : av2 = avma;
2957 0 : H = NULL;
2958 0 : while ((p = u_forprime_next(&S)))
2959 : {
2960 : GEN Mp, B, Hr;
2961 0 : Mp = ZM_to_Flm(M,p);
2962 0 : Hp = Flm_inv_sp(Mp, NULL, p);
2963 0 : if (!Hp) continue;
2964 0 : if (!H)
2965 : {
2966 0 : H = ZM_init_CRT(Hp, p);
2967 0 : q = utoipos(p);
2968 : }
2969 : else
2970 0 : ZM_incremental_CRT(&H, Hp, &q, p);
2971 0 : B = sqrti(shifti(q,-1));
2972 0 : Hr = FpM_ratlift(H,q,B,B,NULL);
2973 0 : if (DEBUGLEVEL>5)
2974 0 : timer_printf(&ti,"ZM_inv mod %lu (ratlift=%ld)", p,!!Hr);
2975 0 : if (Hr) {/* DONE ? */
2976 0 : GEN Hl = Q_remove_denom(Hr, pden);
2977 0 : if (ZM_isscalar(ZM_mul(Hl, M), *pden)) { H = Hl; break; }
2978 : }
2979 :
2980 0 : if (gc_needed(av,2))
2981 : {
2982 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
2983 0 : gerepileall(av2, 2, &H, &q);
2984 : }
2985 : }
2986 0 : if (!*pden) *pden = gen_1;
2987 0 : return gc_all(av, 2, &H, pden);
2988 : }
2989 :
2990 : GEN
2991 74464 : FpM_ratlift_worker(GEN A, GEN mod, GEN B)
2992 : {
2993 : long l, i;
2994 74464 : GEN H = cgetg_copy(A, &l);
2995 156576 : for (i = 1; i < l; i++)
2996 : {
2997 82119 : GEN c = FpC_ratlift(gel(A,i), mod, B, B, NULL);
2998 82114 : gel(H,i) = c? c: gen_0;
2999 : }
3000 74457 : return H;
3001 : }
3002 : static int
3003 762186 : can_ratlift(GEN x, GEN mod, GEN B)
3004 : {
3005 762186 : pari_sp av = avma;
3006 : GEN a, b;
3007 762186 : return gc_bool(av, Fp_ratlift(x, mod, B, B, &a,&b));
3008 : }
3009 : static GEN
3010 2739507 : FpM_ratlift_parallel(GEN A, GEN mod, GEN B)
3011 : {
3012 2739507 : pari_sp av = avma;
3013 : GEN worker;
3014 2739507 : long i, l = lg(A), m = mt_nbthreads();
3015 2739507 : int test = !!B;
3016 :
3017 2739507 : if (l == 1 || lgcols(A) == 1) return gcopy(A);
3018 2739508 : if (!B) B = sqrti(shifti(mod,-1));
3019 2739454 : if (m == 1 || l == 2 || lgcols(A) < 10)
3020 : {
3021 2732005 : A = FpM_ratlift(A, mod, B, B, NULL);
3022 2732053 : return A? A: gc_NULL(av);
3023 : }
3024 : /* test one coefficient first */
3025 7449 : if (test && !can_ratlift(gcoeff(A,1,1), mod, B)) return gc_NULL(av);
3026 7331 : worker = snm_closure(is_entry("_FpM_ratlift_worker"), mkvec2(mod,B));
3027 7331 : A = gen_parapply_slice(worker, A, m);
3028 81273 : for (i = 1; i < l; i++) if (typ(gel(A,i)) != t_COL) return gc_NULL(av);
3029 6337 : return A;
3030 : }
3031 :
3032 : static GEN
3033 755057 : ZM_adj_ratlift(GEN A, GEN H, GEN mod, GEN T)
3034 : {
3035 755057 : pari_sp av = avma;
3036 : GEN B, D, g;
3037 755057 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3038 755056 : if (T) D = mulii(T, D);
3039 755056 : g = gcdii(D, mod);
3040 755051 : if (!equali1(g))
3041 : {
3042 14 : mod = diviiexact(mod, g);
3043 14 : H = FpM_red(H, mod);
3044 : }
3045 755049 : D = Fp_inv(Fp_red(D, mod), mod);
3046 : /* test 1 coeff first */
3047 755049 : B = sqrti(shifti(mod,-1));
3048 755044 : if (!can_ratlift(Fp_mul(D, gcoeff(A,1,1), mod), mod, B)) return gc_NULL(av);
3049 734200 : H = FpM_Fp_mul(H, D, mod);
3050 734193 : H = FpM_ratlift_parallel(H, mod, B);
3051 734200 : return H? H: gc_NULL(av);
3052 : }
3053 :
3054 : /* if (T) return T A^(-1) in Mn(Q), else B in Mn(Z) such that A B = den*Id */
3055 : static GEN
3056 2129955 : ZM_inv_i(GEN A, GEN *pden, GEN T)
3057 : {
3058 2129955 : pari_sp av = avma;
3059 2129955 : long m = lg(A)-1, n, k1 = 1, k2;
3060 2129955 : GEN H = NULL, D, H1 = NULL, mod1 = NULL, worker;
3061 : ulong bnd, mask;
3062 : forprime_t S;
3063 : pari_timer ti;
3064 :
3065 2129955 : if (m == 0) return ZM_inv0(A,pden);
3066 2086451 : if (pden) *pden = gen_1;
3067 2086451 : if (nbrows(A) < m) return NULL;
3068 2086445 : if (m == 1 && nbrows(A)==1 && !T) return ZM_inv1(A,pden);
3069 1442333 : if (m == 2 && nbrows(A)==2 && !T) return ZM_inv2(A,pden);
3070 :
3071 717454 : if (DEBUGLEVEL>=5) timer_start(&ti);
3072 717454 : init_modular_big(&S);
3073 717455 : bnd = expi(RgM_true_Hadamard(A));
3074 717455 : worker = snm_closure(is_entry("_ZM_inv_worker"), mkvec(A));
3075 717455 : gen_inccrt("ZM_inv_r", worker, NULL, k1, 0, &S, &H1, &mod1, nmV_chinese_center, FpM_center);
3076 717456 : n = (bnd+1)/expu(S.p)+1;
3077 717456 : if (DEBUGLEVEL>=5) timer_printf(&ti,"inv (%ld/%ld primes)", k1, n);
3078 717456 : mask = quadratic_prec_mask(n);
3079 717456 : for (k2 = 0;;)
3080 65699 : {
3081 : GEN Hr;
3082 783155 : if (k2 > 0)
3083 : {
3084 58411 : gen_inccrt("ZM_inv_r", worker, NULL, k2, 0, &S, &H1, &mod1,nmV_chinese_center,FpM_center);
3085 58411 : k1 += k2;
3086 58411 : if (DEBUGLEVEL>=5) timer_printf(&ti,"CRT (%ld/%ld primes)", k1, n);
3087 : }
3088 783155 : if (mask == 1) break;
3089 755057 : k2 = (mask&1UL) ? k1-1: k1;
3090 755057 : mask >>= 1;
3091 :
3092 755057 : Hr = ZM_adj_ratlift(A, H1, mod1, T);
3093 755055 : if (DEBUGLEVEL>=5) timer_printf(&ti,"ratlift (%ld/%ld primes)", k1, n);
3094 755055 : if (Hr) {/* DONE ? */
3095 693204 : GEN Hl = Q_primpart(Hr), R = ZM_mul(Hl, A), d = gcoeff(R,1,1);
3096 693202 : if (gsigne(d) < 0) { d = gneg(d); Hl = ZM_neg(Hl); }
3097 693203 : if (DEBUGLEVEL>=5) timer_printf(&ti,"mult (%ld/%ld primes)", k1, n);
3098 693203 : if (equali1(d))
3099 : {
3100 595725 : if (ZM_isidentity(R)) { H = Hl; break; }
3101 : }
3102 97478 : else if (ZM_isscalar(R, d))
3103 : {
3104 93633 : if (T) T = gdiv(T,d);
3105 87334 : else if (pden) *pden = d;
3106 93633 : H = Hl; break;
3107 : }
3108 : }
3109 : }
3110 717455 : if (!H)
3111 : {
3112 : GEN d;
3113 28098 : H = H1;
3114 28098 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3115 28098 : if (signe(D)==0) pari_err_INV("ZM_inv", A);
3116 28098 : if (T) T = gdiv(T, D);
3117 : else
3118 : {
3119 27052 : d = gcdii(Q_content_safe(H), D);
3120 27052 : if (signe(D) < 0) d = negi(d);
3121 27052 : if (!equali1(d))
3122 : {
3123 15332 : H = ZM_Z_divexact(H, d);
3124 15332 : D = diviiexact(D, d);
3125 : }
3126 27052 : if (pden) *pden = D;
3127 : }
3128 : }
3129 717455 : if (T && !isint1(T)) H = ZM_Q_mul(H, T);
3130 717455 : return gc_all(av, pden? 2: 1, &H, pden);
3131 : }
3132 : GEN
3133 2065671 : ZM_inv(GEN A, GEN *pden) { return ZM_inv_i(A, pden, NULL); }
3134 :
3135 : /* same as above, M rational */
3136 : GEN
3137 64283 : QM_inv(GEN M)
3138 : {
3139 64283 : pari_sp av = avma;
3140 : GEN den, dM, K;
3141 64283 : M = Q_remove_denom(M, &dM);
3142 64283 : K = ZM_inv_i(M, &den, dM);
3143 64283 : if (!K) return gc_NULL(av);
3144 64262 : if (den && !equali1(den)) K = ZM_Q_mul(K, ginv(den));
3145 64248 : return gerepileupto(av, K);
3146 : }
3147 :
3148 : static GEN
3149 105226 : ZM_ker_filter(GEN A, GEN P)
3150 : {
3151 105226 : long i, j, l = lg(A), n = 1, d = lg(gmael(A,1,1));
3152 105226 : GEN B, Q, D = gmael(A,1,2);
3153 215193 : for (i=2; i<l; i++)
3154 : {
3155 109967 : GEN Di = gmael(A,i,2);
3156 109967 : long di = lg(gmael(A,i,1));
3157 109967 : int c = vecsmall_lexcmp(D, Di);
3158 109967 : if (di==d && c==0) n++;
3159 45588 : else if (d > di || (di==d && c>0))
3160 37680 : { n = 1; d = di; D = Di; }
3161 : }
3162 105226 : B = cgetg(n+1, t_VEC);
3163 105226 : Q = cgetg(n+1, typ(P));
3164 320419 : for (i=1, j=1; i<l; i++)
3165 : {
3166 215193 : if (lg(gmael(A,i,1))==d && vecsmall_lexcmp(D, gmael(A,i,2))==0)
3167 : {
3168 169605 : gel(B,j) = gmael(A,i,1);
3169 169605 : Q[j] = P[i];
3170 169605 : j++;
3171 : }
3172 : }
3173 105226 : return mkvec3(B,Q,D);
3174 : }
3175 :
3176 : static GEN
3177 69559 : ZM_ker_chinese(GEN A, GEN P, GEN *mod)
3178 : {
3179 69559 : GEN BQD = ZM_ker_filter(A, P);
3180 69559 : return mkvec2(nmV_chinese_center(gel(BQD,1), gel(BQD,2), mod), gel(BQD,3));
3181 : }
3182 :
3183 : static GEN
3184 131838 : ZM_ker_slice(GEN A, GEN P, GEN *mod)
3185 : {
3186 131838 : pari_sp av = avma;
3187 131838 : long i, n = lg(P)-1;
3188 : GEN BQD, B, Q, D, H, HD, T;
3189 131838 : if (n == 1)
3190 : {
3191 96171 : ulong p = uel(P,1);
3192 96171 : GEN K = Flm_ker_sp(ZM_to_Flm(A, p), p, 2);
3193 96171 : *mod = utoipos(p); return mkvec2(Flm_to_ZM(gel(K,1)), gel(K,2));
3194 : }
3195 35667 : T = ZV_producttree(P);
3196 35667 : A = ZM_nv_mod_tree(A, P, T);
3197 35667 : H = cgetg(n+1, t_VEC);
3198 111507 : for(i=1 ; i <= n; i++)
3199 75840 : gel(H,i) = Flm_ker_sp(gel(A, i), P[i], 2);
3200 35667 : BQD = ZM_ker_filter(H, P);
3201 35667 : B = gel(BQD,1); Q = gel(BQD,2); D = gel(BQD, 3);
3202 35667 : if (lg(Q) != lg(P)) T = ZV_producttree(Q);
3203 35667 : H = nmV_chinese_center_tree_seq(B, Q, T, ZV_chinesetree(Q,T));
3204 35667 : *mod = gmael(T, lg(T)-1, 1);
3205 35667 : HD = mkvec2(H, D);
3206 35667 : return gc_all(av, 2, &HD, mod);
3207 : }
3208 :
3209 : GEN
3210 131838 : ZM_ker_worker(GEN P, GEN A)
3211 : {
3212 131838 : GEN V = cgetg(3, t_VEC);
3213 131838 : gel(V,1) = ZM_ker_slice(A, P, &gel(V,2));
3214 131838 : return V;
3215 : }
3216 :
3217 : /* assume lg(A) > 1 */
3218 : static GEN
3219 65103 : ZM_ker_i(GEN A)
3220 : {
3221 : pari_sp av;
3222 65103 : long k, m = lg(A)-1;
3223 65103 : GEN HD = NULL, mod = gen_1, worker;
3224 : forprime_t S;
3225 :
3226 65103 : if (m >= 2*nbrows(A))
3227 : {
3228 3059 : GEN v = ZM_indexrank(A), y = gel(v,2), z = indexcompl(y, m);
3229 : GEN B, A1, A1i, d;
3230 3059 : A = rowpermute(A, gel(v,1)); /* same kernel */
3231 3059 : A1 = vecpermute(A, y); /* maximal rank submatrix */
3232 3059 : B = vecpermute(A, z);
3233 3059 : A1i = ZM_inv(A1, &d);
3234 3059 : if (!d) d = gen_1;
3235 3059 : B = vconcat(ZM_mul(ZM_neg(A1i), B), scalarmat_shallow(d, lg(B)-1));
3236 3059 : if (!gequal(y, identity_perm(lg(y)-1)))
3237 685 : B = rowpermute(B, perm_inv(shallowconcat(y,z)));
3238 3059 : return vec_Q_primpart(B);
3239 : }
3240 62044 : init_modular_big(&S);
3241 62044 : worker = snm_closure(is_entry("_ZM_ker_worker"), mkvec(A));
3242 62044 : av = avma;
3243 62044 : for (k = 1;; k <<= 1)
3244 65346 : {
3245 : pari_timer ti;
3246 : GEN H, Hr;
3247 127390 : gen_inccrt_i("ZM_ker", worker, NULL, (k+1)>>1, 0,
3248 : &S, &HD, &mod, ZM_ker_chinese, NULL);
3249 127390 : gerepileall(av, 2, &HD, &mod);
3250 143017 : H = gel(HD, 1); if (lg(H) == 1) return H;
3251 80973 : if (DEBUGLEVEL >= 4) timer_start(&ti);
3252 80973 : Hr = FpM_ratlift_parallel(H, mod, NULL);
3253 80973 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: ratlift (%ld)",!!Hr);
3254 80973 : if (Hr)
3255 : {
3256 : GEN MH;
3257 70223 : Hr = vec_Q_primpart(Hr);
3258 70223 : MH = ZM_mul(A, Hr);
3259 70223 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: QM_mul");
3260 70223 : if (ZM_equal0(MH)) return Hr;
3261 : }
3262 : }
3263 : }
3264 :
3265 : GEN
3266 49269 : ZM_ker(GEN M)
3267 : {
3268 49269 : pari_sp av = avma;
3269 49269 : long l = lg(M)-1;
3270 49269 : if (l==0) return cgetg(1, t_MAT);
3271 49269 : if (lgcols(M)==1) return matid(l);
3272 49269 : return gerepilecopy(av, ZM_ker_i(M));
3273 : }
3274 :
3275 : static GEN
3276 1999262 : ZM_gauss_slice(GEN A, GEN B, GEN P, GEN *mod)
3277 : {
3278 1999262 : pari_sp av = avma;
3279 1999262 : long i, n = lg(P)-1;
3280 : GEN H, T;
3281 1999262 : if (n == 1)
3282 : {
3283 1959750 : ulong p = uel(P,1);
3284 1959750 : GEN Hp = Flm_gauss(ZM_to_Flm(A, p) , ZM_to_Flm(B, p) ,p);
3285 1959751 : if (!Hp) { *mod=gen_1; return zeromat(lg(A)-1,lg(B)-1); }
3286 1959751 : Hp = gerepileupto(av, Flm_to_ZM(Hp));
3287 1959755 : *mod = utoipos(p); return Hp;
3288 : }
3289 39512 : T = ZV_producttree(P);
3290 39512 : A = ZM_nv_mod_tree(A, P, T);
3291 39512 : B = ZM_nv_mod_tree(B, P, T);
3292 39512 : H = cgetg(n+1, t_VEC);
3293 223420 : for(i=1; i <= n; i++)
3294 : {
3295 183908 : GEN Hi = Flm_gauss(gel(A, i), gel(B,i), uel(P,i));
3296 183908 : gel(H,i) = Hi ? Hi: zero_Flm(lg(A)-1,lg(B)-1);
3297 183908 : if (!Hi) uel(P,i)=1;
3298 : }
3299 39512 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
3300 39512 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3301 : }
3302 :
3303 : GEN
3304 1999261 : ZM_gauss_worker(GEN P, GEN A, GEN B)
3305 : {
3306 1999261 : GEN V = cgetg(3, t_VEC);
3307 1999262 : gel(V,1) = ZM_gauss_slice(A, B, P, &gel(V,2));
3308 1999263 : return V;
3309 : }
3310 :
3311 : /* assume lg(A) > 1 */
3312 : static GEN
3313 1742573 : ZM_gauss_i(GEN A, GEN B)
3314 : {
3315 : pari_sp av;
3316 : long k, m, ncol;
3317 : int iscol;
3318 1742573 : GEN y, y1, y2, Hr, H = NULL, mod = gen_1, worker;
3319 : forprime_t S;
3320 1742573 : if (!init_gauss(A, &B, &m, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
3321 1742504 : init_modular_big(&S);
3322 1742503 : y = ZM_indexrank(A); y1 = gel(y,1); y2 = gel(y,2);
3323 1742506 : if (lg(y2)-1 != m) return NULL;
3324 1742478 : A = rowpermute(A, y1);
3325 1742474 : B = rowpermute(B, y1);
3326 : /* a is square and invertible */
3327 1742474 : ncol = lg(B);
3328 1742474 : worker = snm_closure(is_entry("_ZM_gauss_worker"), mkvec2(A,B));
3329 1742482 : av = avma;
3330 1742482 : for (k = 1;; k <<= 1)
3331 181868 : {
3332 : pari_timer ti;
3333 1924350 : gen_inccrt_i("ZM_gauss", worker, NULL, (k+1)>>1 , m,
3334 : &S, &H, &mod, nmV_chinese_center, FpM_center);
3335 1924331 : gerepileall(av, 2, &H, &mod);
3336 1924348 : if (DEBUGLEVEL >= 4) timer_start(&ti);
3337 1924348 : Hr = FpM_ratlift_parallel(H, mod, NULL);
3338 1924333 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_gauss: ratlift (%ld)",!!Hr);
3339 1924338 : if (Hr)
3340 : {
3341 : GEN MH, c;
3342 1788078 : MH = ZM_mul(A, Q_remove_denom(Hr, &c));
3343 1788059 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_gauss: QM_mul");
3344 1788071 : if (ZM_equal(MH, c ? ZM_Z_mul(B, c): B)) break;
3345 : }
3346 : }
3347 1742471 : return iscol ? gel(Hr, 1): Hr;
3348 : }
3349 :
3350 : GEN
3351 1742573 : ZM_gauss(GEN A, GEN B)
3352 : {
3353 1742573 : pari_sp av = avma;
3354 1742573 : GEN C = ZM_gauss_i(A,B);
3355 1742565 : return C ? gerepilecopy(av, C): NULL;
3356 : }
3357 :
3358 : GEN
3359 16709 : QM_ker(GEN M)
3360 : {
3361 16709 : pari_sp av = avma;
3362 16709 : long l = lg(M)-1;
3363 16709 : if (l==0) return cgetg(1, t_MAT);
3364 16667 : if (lgcols(M)==1) return matid(l);
3365 15750 : return gerepilecopy(av, ZM_ker_i(row_Q_primpart(M)));
3366 : }
3367 :
3368 : /* x a ZM. Return a multiple of the determinant of the lattice generated by
3369 : * the columns of x. From Algorithm 2.2.6 in GTM138 */
3370 : GEN
3371 49964 : detint(GEN A)
3372 : {
3373 49964 : if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
3374 49964 : RgM_check_ZM(A, "detint");
3375 49964 : return ZM_detmult(A);
3376 : }
3377 : GEN
3378 169278 : ZM_detmult(GEN A)
3379 : {
3380 169278 : pari_sp av1, av = avma;
3381 : GEN B, c, v, piv;
3382 169278 : long rg, i, j, k, m, n = lg(A) - 1;
3383 :
3384 169278 : if (!n) return gen_1;
3385 169278 : m = nbrows(A);
3386 169278 : if (n < m) return gen_0;
3387 169201 : c = zero_zv(m);
3388 169201 : av1 = avma;
3389 169201 : B = zeromatcopy(m,m);
3390 169202 : v = cgetg(m+1, t_COL);
3391 169202 : piv = gen_1; rg = 0;
3392 726471 : for (k=1; k<=n; k++)
3393 : {
3394 726457 : GEN pivprec = piv;
3395 726457 : long t = 0;
3396 5363810 : for (i=1; i<=m; i++)
3397 : {
3398 4637352 : pari_sp av2 = avma;
3399 : GEN vi;
3400 4637352 : if (c[i]) continue;
3401 :
3402 2682155 : vi = mulii(piv, gcoeff(A,i,k));
3403 28372491 : for (j=1; j<=m; j++)
3404 25690323 : if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
3405 2682168 : if (!t && signe(vi)) t = i;
3406 2682168 : gel(v,i) = gerepileuptoint(av2, vi);
3407 : }
3408 726458 : if (!t) continue;
3409 : /* at this point c[t] = 0 */
3410 :
3411 726367 : if (++rg >= m) { /* full rank; mostly done */
3412 169188 : GEN det = gel(v,t); /* last on stack */
3413 169188 : if (++k > n)
3414 169056 : det = absi(det);
3415 : else
3416 : {
3417 : /* improve further; at this point c[i] is set for all i != t */
3418 132 : gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
3419 418 : for ( ; k<=n; k++)
3420 286 : det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
3421 : }
3422 169188 : return gerepileuptoint(av, det);
3423 : }
3424 :
3425 557179 : piv = gel(v,t);
3426 4467671 : for (i=1; i<=m; i++)
3427 : {
3428 : GEN mvi;
3429 3910492 : if (c[i] || i == t) continue;
3430 :
3431 1955246 : gcoeff(B,t,i) = mvi = negi(gel(v,i));
3432 23004783 : for (j=1; j<=m; j++)
3433 21049537 : if (c[j]) /* implies j != t */
3434 : {
3435 5713015 : pari_sp av2 = avma;
3436 5713015 : GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
3437 5713015 : if (rg > 1) z = diviiexact(z, pivprec);
3438 5713015 : gcoeff(B,j,i) = gerepileuptoint(av2, z);
3439 : }
3440 : }
3441 557179 : c[t] = k;
3442 557179 : if (gc_needed(av,1))
3443 : {
3444 0 : if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
3445 0 : gerepileall(av1, 2, &piv,&B); v = zerovec(m);
3446 : }
3447 : }
3448 14 : return gc_const(av, gen_0);
3449 : }
3450 :
3451 : /* Reduce x modulo (invertible) y */
3452 : GEN
3453 14979 : closemodinvertible(GEN x, GEN y)
3454 : {
3455 14979 : return gmul(y, ground(RgM_solve(y,x)));
3456 : }
3457 : GEN
3458 7 : reducemodinvertible(GEN x, GEN y)
3459 : {
3460 7 : return gsub(x, closemodinvertible(x,y));
3461 : }
3462 : GEN
3463 0 : reducemodlll(GEN x,GEN y)
3464 : {
3465 0 : return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
3466 : }
3467 :
3468 : /*******************************************************************/
3469 : /* */
3470 : /* KERNEL of an m x n matrix */
3471 : /* return n - rk(x) linearly independent vectors */
3472 : /* */
3473 : /*******************************************************************/
3474 : static GEN
3475 28 : RgM_deplin_i(GEN x0)
3476 : {
3477 28 : pari_sp av = avma, av2;
3478 28 : long i, j, k, nl, nc = lg(x0)-1;
3479 : GEN D, x, y, c, l, d, ck;
3480 :
3481 28 : if (!nc) return NULL;
3482 28 : nl = nbrows(x0);
3483 28 : c = zero_zv(nl);
3484 28 : l = cgetg(nc+1, t_VECSMALL); /* not initialized */
3485 28 : av2 = avma;
3486 28 : x = RgM_shallowcopy(x0);
3487 28 : d = const_vec(nl, gen_1); /* pivot list */
3488 28 : ck = NULL; /* gcc -Wall */
3489 98 : for (k=1; k<=nc; k++)
3490 : {
3491 91 : ck = gel(x,k);
3492 196 : for (j=1; j<k; j++)
3493 : {
3494 105 : GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
3495 420 : for (i=1; i<=nl; i++)
3496 315 : if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
3497 : }
3498 :
3499 91 : i = gauss_get_pivot_NZ(x, NULL, k, c);
3500 91 : if (i > nl) break;
3501 70 : if (gc_needed(av,1))
3502 : {
3503 0 : if (DEBUGMEM>1) pari_warn(warnmem,"deplin k = %ld/%ld",k,nc);
3504 0 : gerepileall(av2, 2, &x, &d);
3505 0 : ck = gel(x,k);
3506 : }
3507 70 : gel(d,k) = gel(ck,i);
3508 70 : c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
3509 : }
3510 28 : if (k > nc) return gc_NULL(av);
3511 21 : if (k == 1) { set_avma(av); return scalarcol_shallow(gen_1,nc); }
3512 21 : y = cgetg(nc+1,t_COL);
3513 21 : gel(y,1) = gcopy(gel(ck, l[1]));
3514 49 : for (D=gel(d,1),j=2; j<k; j++)
3515 : {
3516 28 : gel(y,j) = gmul(gel(ck, l[j]), D);
3517 28 : D = gmul(D, gel(d,j));
3518 : }
3519 21 : gel(y,j) = gneg(D);
3520 21 : for (j++; j<=nc; j++) gel(y,j) = gen_0;
3521 21 : y = primitive_part(y, &c);
3522 21 : return c? gerepileupto(av, y): gerepilecopy(av, y);
3523 : }
3524 : static GEN
3525 0 : RgV_deplin(GEN v)
3526 : {
3527 0 : pari_sp av = avma;
3528 0 : long n = lg(v)-1;
3529 0 : GEN y, p = NULL;
3530 0 : if (n <= 1)
3531 : {
3532 0 : if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
3533 0 : return cgetg(1, t_COL);
3534 : }
3535 0 : if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
3536 0 : v = primpart(mkvec2(gel(v,1),gel(v,2)));
3537 0 : if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
3538 0 : y = zerocol(n);
3539 0 : gel(y,1) = gneg(gel(v,2));
3540 0 : gel(y,2) = gcopy(gel(v,1));
3541 0 : return gerepileupto(av, y);
3542 :
3543 : }
3544 :
3545 : static GEN
3546 105 : RgM_deplin_FpM(GEN x, GEN p)
3547 : {
3548 105 : pari_sp av = avma;
3549 : ulong pp;
3550 105 : x = RgM_Fp_init3(x, p, &pp);
3551 105 : switch(pp)
3552 : {
3553 35 : case 0:
3554 35 : x = FpM_ker_gen(x,p,1);
3555 35 : if (!x) return gc_NULL(av);
3556 21 : x = FpC_center(x,p,shifti(p,-1));
3557 21 : break;
3558 14 : case 2:
3559 14 : x = F2m_ker_sp(x,1);
3560 14 : if (!x) return gc_NULL(av);
3561 7 : x = F2c_to_ZC(x); break;
3562 0 : case 3:
3563 0 : x = F3m_ker_sp(x,1);
3564 0 : if (!x) return gc_NULL(av);
3565 0 : x = F3c_to_ZC(x); break;
3566 56 : default:
3567 56 : x = Flm_ker_sp(x,pp,1);
3568 56 : if (!x) return gc_NULL(av);
3569 35 : x = Flv_center(x, pp, pp>>1);
3570 35 : x = zc_to_ZC(x);
3571 35 : break;
3572 : }
3573 63 : return gerepileupto(av, x);
3574 : }
3575 :
3576 : /* FIXME: implement direct modular ZM_deplin ? */
3577 : static GEN
3578 119 : QM_deplin(GEN M)
3579 : {
3580 119 : pari_sp av = avma;
3581 119 : long l = lg(M)-1;
3582 : GEN k;
3583 119 : if (l==0) return NULL;
3584 84 : if (lgcols(M)==1) return col_ei(l, 1);
3585 84 : k = ZM_ker_i(row_Q_primpart(M));
3586 84 : if (lg(k)== 1) return gc_NULL(av);
3587 70 : return gerepilecopy(av, gel(k,1));
3588 : }
3589 :
3590 : static GEN
3591 49 : RgM_deplin_FqM(GEN x, GEN pol, GEN p)
3592 : {
3593 49 : pari_sp av = avma;
3594 49 : GEN b, T = RgX_to_FpX(pol, p);
3595 49 : if (signe(T) == 0) pari_err_OP("deplin",x,pol);
3596 49 : b = FqM_deplin(RgM_to_FqM(x, T, p), T, p);
3597 49 : if (!b) return gc_NULL(av);
3598 35 : return gerepileupto(av, b);
3599 : }
3600 :
3601 : #define code(t1,t2) ((t1 << 6) | t2)
3602 : static GEN
3603 385 : RgM_deplin_fast(GEN x)
3604 : {
3605 : GEN p, pol;
3606 : long pa;
3607 385 : long t = RgM_type(x, &p,&pol,&pa);
3608 385 : switch(t)
3609 : {
3610 119 : case t_INT: /* fall through */
3611 119 : case t_FRAC: return QM_deplin(x);
3612 84 : case t_FFELT: return FFM_deplin(x, pol);
3613 105 : case t_INTMOD: return RgM_deplin_FpM(x, p);
3614 49 : case code(t_POLMOD, t_INTMOD):
3615 49 : return RgM_deplin_FqM(x, pol, p);
3616 28 : default: return gen_0;
3617 : }
3618 : }
3619 : #undef code
3620 :
3621 : static GEN
3622 385 : RgM_deplin(GEN x)
3623 : {
3624 385 : GEN z = RgM_deplin_fast(x);
3625 385 : if (z!= gen_0) return z;
3626 28 : return RgM_deplin_i(x);
3627 : }
3628 :
3629 : GEN
3630 385 : deplin(GEN x)
3631 : {
3632 385 : switch(typ(x))
3633 : {
3634 385 : case t_MAT:
3635 : {
3636 385 : GEN z = RgM_deplin(x);
3637 385 : if (z) return z;
3638 147 : return cgetg(1, t_COL);
3639 : }
3640 0 : case t_VEC: return RgV_deplin(x);
3641 0 : default: pari_err_TYPE("deplin",x);
3642 : }
3643 : return NULL;/*LCOV_EXCL_LINE*/
3644 : }
3645 :
3646 : /*******************************************************************/
3647 : /* */
3648 : /* GAUSS REDUCTION OF MATRICES (m lines x n cols) */
3649 : /* (kernel, image, complementary image, rank) */
3650 : /* */
3651 : /*******************************************************************/
3652 : /* return the transform of x under a standard Gauss pivot.
3653 : * x0 is a reference point when guessing whether x[i,j] ~ 0
3654 : * (iff x[i,j] << x0[i,j])
3655 : * Set r = dim ker(x). d[k] contains the index of the first nonzero pivot
3656 : * in column k */
3657 : static GEN
3658 1271 : gauss_pivot_ker(GEN x, GEN x0, GEN *dd, long *rr)
3659 : {
3660 : GEN c, d, p, data;
3661 : pari_sp av;
3662 : long i, j, k, r, t, n, m;
3663 : pivot_fun pivot;
3664 :
3665 1271 : n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
3666 1271 : m=nbrows(x); r=0;
3667 1271 : pivot = get_pivot_fun(x, x0, &data);
3668 1271 : x = RgM_shallowcopy(x);
3669 1271 : c = zero_zv(m);
3670 1271 : d = cgetg(n+1,t_VECSMALL);
3671 1271 : av=avma;
3672 7475 : for (k=1; k<=n; k++)
3673 : {
3674 6204 : j = pivot(x, data, k, c);
3675 6204 : if (j > m)
3676 : {
3677 1463 : r++; d[k]=0;
3678 6496 : for(j=1; j<k; j++)
3679 5033 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
3680 : }
3681 : else
3682 : { /* pivot for column k on row j */
3683 4741 : c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
3684 4741 : gcoeff(x,j,k) = gen_m1;
3685 : /* x[j,] /= - x[j,k] */
3686 24169 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3687 42136 : for (t=1; t<=m; t++)
3688 37395 : if (t!=j)
3689 : { /* x[t,] -= 1 / x[j,k] x[j,] */
3690 32654 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3691 32654 : if (gequal0(p)) continue;
3692 86920 : for (i=k+1; i<=n; i++)
3693 69463 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
3694 17457 : if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
3695 : }
3696 : }
3697 : }
3698 1271 : *dd=d; *rr=r; return x;
3699 : }
3700 :
3701 : /* r = dim ker(x).
3702 : * Returns d:
3703 : * d[k] != 0 contains the index of a nonzero pivot in column k
3704 : * d[k] == 0 if column k is a linear combination of the (k-1) first ones */
3705 : GEN
3706 167558 : RgM_pivots(GEN x0, GEN data, long *rr, pivot_fun pivot)
3707 : {
3708 : GEN x, c, d, p;
3709 167558 : long i, j, k, r, t, m, n = lg(x0)-1;
3710 : pari_sp av;
3711 :
3712 167558 : if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
3713 152266 : if (!n) { *rr = 0; return NULL; }
3714 :
3715 152266 : d = cgetg(n+1, t_VECSMALL);
3716 152266 : x = RgM_shallowcopy(x0);
3717 152265 : m = nbrows(x); r = 0;
3718 152265 : c = zero_zv(m);
3719 152270 : av = avma;
3720 927602 : for (k=1; k<=n; k++)
3721 : {
3722 775336 : j = pivot(x, data, k, c);
3723 775348 : if (j > m) { r++; d[k] = 0; }
3724 : else
3725 : {
3726 291495 : c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
3727 1882148 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3728 :
3729 1054025 : for (t=1; t<=m; t++)
3730 762546 : if (!c[t]) /* no pivot on that line yet */
3731 : {
3732 257071 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3733 4126159 : for (i=k+1; i<=n; i++)
3734 3869091 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
3735 257068 : if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
3736 : }
3737 2173723 : for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
3738 : }
3739 : }
3740 152266 : *rr = r; return gc_const((pari_sp)d, d);
3741 : }
3742 :
3743 : static long
3744 4244542 : ZM_count_0_cols(GEN M)
3745 : {
3746 4244542 : long i, l = lg(M), n = 0;
3747 18206982 : for (i = 1; i < l; i++)
3748 13962447 : if (ZV_equal0(gel(M,i))) n++;
3749 4244535 : return n;
3750 : }
3751 :
3752 : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
3753 : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
3754 : GEN
3755 4258196 : ZM_pivots(GEN M0, long *rr)
3756 : {
3757 4258196 : GEN d, dbest = NULL;
3758 : long m, mm, n, nn, i, imax, rmin, rbest, zc;
3759 4258196 : int beenthere = 0;
3760 4258196 : pari_sp av, av0 = avma;
3761 : forprime_t S;
3762 :
3763 4258196 : rbest = n = lg(M0)-1;
3764 4258196 : if (n == 0) { *rr = 0; return NULL; }
3765 4244545 : zc = ZM_count_0_cols(M0);
3766 4244537 : if (n == zc) { *rr = zc; return zero_zv(n); }
3767 :
3768 4244408 : m = nbrows(M0);
3769 4244409 : rmin = maxss(zc, n-m);
3770 4244404 : init_modular_small(&S);
3771 4244413 : if (n <= m) { nn = n; mm = m; } else { nn = m; mm = n; }
3772 4244413 : imax = (nn < 16)? 1: (nn < 64)? 2: 3; /* heuristic */
3773 :
3774 : for(;;)
3775 0 : {
3776 : GEN row, col, M, KM, IM, RHS, X, cX;
3777 : long rk;
3778 4267581 : for (av = avma, i = 0;; set_avma(av), i++)
3779 23174 : {
3780 4267581 : ulong p = u_forprime_next(&S);
3781 : long rp;
3782 4267558 : if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
3783 4267558 : d = Flm_pivots(ZM_to_Flm(M0, p), p, &rp, 1);
3784 4267584 : if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
3785 44961 : if (rp < rbest) { /* save best r so far */
3786 21812 : rbest = rp;
3787 21812 : guncloneNULL(dbest);
3788 21812 : dbest = gclone(d);
3789 21812 : if (beenthere) break;
3790 : }
3791 44961 : if (!beenthere && i >= imax) break;
3792 : }
3793 21787 : beenthere = 1;
3794 : /* Dubious case: there is (probably) a non trivial kernel */
3795 21787 : indexrank_all(m,n, rbest, dbest, &row, &col);
3796 21787 : M = rowpermute(vecpermute(M0, col), row);
3797 21787 : rk = n - rbest; /* (probable) dimension of image */
3798 21787 : if (n > m) M = shallowtrans(M);
3799 21787 : IM = vecslice(M,1,rk);
3800 21787 : KM = vecslice(M,rk+1, nn);
3801 21787 : M = rowslice(IM, 1,rk); /* square maximal rank */
3802 21787 : X = ZM_gauss(M, rowslice(KM, 1,rk));
3803 21787 : RHS = rowslice(KM,rk+1,mm);
3804 21787 : M = rowslice(IM,rk+1,mm);
3805 21787 : X = Q_remove_denom(X, &cX);
3806 21787 : if (cX) RHS = ZM_Z_mul(RHS, cX);
3807 21787 : if (ZM_equal(ZM_mul(M, X), RHS)) { d = vecsmall_copy(dbest); goto END; }
3808 0 : set_avma(av);
3809 : }
3810 4244410 : END:
3811 4244410 : *rr = rbest; guncloneNULL(dbest);
3812 4244409 : return gerepileuptoleaf(av0, d);
3813 : }
3814 :
3815 : /* set *pr = dim Ker x */
3816 : static GEN
3817 75888 : gauss_pivot(GEN x, long *pr) {
3818 : GEN data;
3819 75888 : pivot_fun pivot = get_pivot_fun(x, x, &data);
3820 75889 : return RgM_pivots(x, data, pr, pivot);
3821 : }
3822 :
3823 : /* compute ker(x), x0 is a reference point when guessing whether x[i,j] ~ 0
3824 : * (iff x[i,j] << x0[i,j]) */
3825 : static GEN
3826 1271 : ker_aux(GEN x, GEN x0)
3827 : {
3828 1271 : pari_sp av = avma;
3829 : GEN d,y;
3830 : long i,j,k,r,n;
3831 :
3832 1271 : x = gauss_pivot_ker(x,x0,&d,&r);
3833 1271 : if (!r) { set_avma(av); return cgetg(1,t_MAT); }
3834 1211 : n = lg(x)-1; y=cgetg(r+1,t_MAT);
3835 2674 : for (j=k=1; j<=r; j++,k++)
3836 : {
3837 1463 : GEN p = cgetg(n+1,t_COL);
3838 :
3839 5586 : gel(y,j) = p; while (d[k]) k++;
3840 6496 : for (i=1; i<k; i++)
3841 5033 : if (d[i])
3842 : {
3843 4641 : GEN p1=gcoeff(x,d[i],k);
3844 4641 : gel(p,i) = gcopy(p1); gunclone(p1);
3845 : }
3846 : else
3847 392 : gel(p,i) = gen_0;
3848 2541 : gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
3849 : }
3850 1211 : return gerepileupto(av,y);
3851 : }
3852 :
3853 : static GEN
3854 553 : RgM_ker_FpM(GEN x, GEN p)
3855 : {
3856 553 : pari_sp av = avma;
3857 : ulong pp;
3858 553 : x = RgM_Fp_init3(x, p, &pp);
3859 553 : switch(pp)
3860 : {
3861 35 : case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
3862 21 : case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
3863 77 : case 3: x = F3m_to_mod(F3m_ker_sp(x,0)); break;
3864 420 : default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
3865 : }
3866 553 : return gerepileupto(av, x);
3867 : }
3868 :
3869 : static GEN
3870 91 : RgM_ker_FqM(GEN x, GEN pol, GEN p)
3871 : {
3872 91 : pari_sp av = avma;
3873 91 : GEN b, T = RgX_to_FpX(pol, p);
3874 91 : if (signe(T) == 0) pari_err_OP("ker",x,pol);
3875 84 : b = FqM_ker(RgM_to_FqM(x, T, p), T, p);
3876 84 : return gerepileupto(av, FqM_to_mod(b, T, p));
3877 : }
3878 :
3879 : #define code(t1,t2) ((t1 << 6) | t2)
3880 : static GEN
3881 9198 : RgM_ker_fast(GEN x)
3882 : {
3883 : GEN p, pol;
3884 : long pa;
3885 9198 : long t = RgM_type(x, &p,&pol,&pa);
3886 9198 : switch(t)
3887 : {
3888 7609 : case t_INT: /* fall through */
3889 7609 : case t_FRAC: return QM_ker(x);
3890 63 : case t_FFELT: return FFM_ker(x, pol);
3891 553 : case t_INTMOD: return RgM_ker_FpM(x, p);
3892 91 : case code(t_POLMOD, t_INTMOD):
3893 91 : return RgM_ker_FqM(x, pol, p);
3894 882 : default: return NULL;
3895 : }
3896 : }
3897 : #undef code
3898 :
3899 : GEN
3900 9198 : ker(GEN x)
3901 : {
3902 9198 : GEN b = RgM_ker_fast(x);
3903 9191 : if (b) return b;
3904 882 : return ker_aux(x,x);
3905 : }
3906 :
3907 : GEN
3908 46221 : matker0(GEN x,long flag)
3909 : {
3910 46221 : if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
3911 46221 : if (!flag) return ker(x);
3912 45934 : RgM_check_ZM(x, "matker");
3913 45934 : return ZM_ker(x);
3914 : }
3915 :
3916 : static GEN
3917 525 : RgM_image_FpM(GEN x, GEN p)
3918 : {
3919 525 : pari_sp av = avma;
3920 : ulong pp;
3921 525 : x = RgM_Fp_init(x, p, &pp);
3922 525 : switch(pp)
3923 : {
3924 28 : case 0: x = FpM_to_mod(FpM_image(x,p),p); break;
3925 7 : case 2: x = F2m_to_mod(F2m_image(x)); break;
3926 490 : default:x = Flm_to_mod(Flm_image(x,pp), pp); break;
3927 : }
3928 525 : return gerepileupto(av, x);
3929 : }
3930 :
3931 : static GEN
3932 35 : RgM_image_FqM(GEN x, GEN pol, GEN p)
3933 : {
3934 35 : pari_sp av = avma;
3935 35 : GEN b, T = RgX_to_FpX(pol, p);
3936 35 : if (signe(T) == 0) pari_err_OP("image",x,pol);
3937 28 : b = FqM_image(RgM_to_FqM(x, T, p), T, p);
3938 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
3939 : }
3940 :
3941 : GEN
3942 6181 : QM_image_shallow(GEN A)
3943 : {
3944 6181 : A = vec_Q_primpart(A);
3945 6181 : return vecpermute(A, ZM_indeximage(A));
3946 : }
3947 : GEN
3948 5411 : QM_image(GEN A)
3949 : {
3950 5411 : pari_sp av = avma;
3951 5411 : return gerepilecopy(av, QM_image_shallow(A));
3952 : }
3953 :
3954 : #define code(t1,t2) ((t1 << 6) | t2)
3955 : static GEN
3956 6034 : RgM_image_fast(GEN x)
3957 : {
3958 : GEN p, pol;
3959 : long pa;
3960 6034 : long t = RgM_type(x, &p,&pol,&pa);
3961 6034 : switch(t)
3962 : {
3963 5411 : case t_INT: /* fall through */
3964 5411 : case t_FRAC: return QM_image(x);
3965 49 : case t_FFELT: return FFM_image(x, pol);
3966 525 : case t_INTMOD: return RgM_image_FpM(x, p);
3967 35 : case code(t_POLMOD, t_INTMOD):
3968 35 : return RgM_image_FqM(x, pol, p);
3969 14 : default: return NULL;
3970 : }
3971 : }
3972 : #undef code
3973 :
3974 : GEN
3975 6034 : image(GEN x)
3976 : {
3977 : GEN d, M;
3978 : long r;
3979 :
3980 6034 : if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
3981 6034 : M = RgM_image_fast(x);
3982 6027 : if (M) return M;
3983 14 : d = gauss_pivot(x,&r); /* d left on stack for efficiency */
3984 14 : return image_from_pivot(x,d,r);
3985 : }
3986 :
3987 : static GEN
3988 84 : imagecompl_aux(GEN x, GEN(*PIVOT)(GEN,long*))
3989 : {
3990 84 : pari_sp av = avma;
3991 : GEN d,y;
3992 : long j,i,r;
3993 :
3994 84 : if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
3995 84 : (void)new_chunk(lg(x) * 4 + 1); /* HACK */
3996 84 : d = PIVOT(x,&r); /* if (!d) then r = 0 */
3997 84 : set_avma(av); y = cgetg(r+1,t_VECSMALL);
3998 126 : for (i=j=1; j<=r; i++)
3999 42 : if (!d[i]) y[j++] = i;
4000 84 : return y;
4001 : }
4002 : GEN
4003 84 : imagecompl(GEN x) { return imagecompl_aux(x, &gauss_pivot); }
4004 : GEN
4005 0 : ZM_imagecompl(GEN x) { return imagecompl_aux(x, &ZM_pivots); }
4006 :
4007 : static GEN
4008 28 : RgM_RgC_invimage_FpC(GEN A, GEN y, GEN p)
4009 : {
4010 28 : pari_sp av = avma;
4011 : ulong pp;
4012 : GEN x;
4013 28 : A = RgM_Fp_init(A,p,&pp);
4014 28 : switch(pp)
4015 : {
4016 7 : case 0:
4017 7 : y = RgC_to_FpC(y,p);
4018 7 : x = FpM_FpC_invimage(A, y, p);
4019 7 : return x ? gerepileupto(av, FpC_to_mod(x,p)): NULL;
4020 7 : case 2:
4021 7 : y = RgV_to_F2v(y);
4022 7 : x = F2m_F2c_invimage(A, y);
4023 7 : return x ? gerepileupto(av, F2c_to_mod(x)): NULL;
4024 14 : default:
4025 14 : y = RgV_to_Flv(y,pp);
4026 14 : x = Flm_Flc_invimage(A, y, pp);
4027 14 : return x ? gerepileupto(av, Flc_to_mod(x,pp)): NULL;
4028 : }
4029 : }
4030 :
4031 : static GEN
4032 2184 : RgM_RgC_invimage_fast(GEN x, GEN y)
4033 : {
4034 : GEN p, pol;
4035 : long pa;
4036 2184 : long t = RgM_RgC_type(x, y, &p,&pol,&pa);
4037 2184 : switch(t)
4038 : {
4039 28 : case t_INTMOD: return RgM_RgC_invimage_FpC(x, y, p);
4040 63 : case t_FFELT: return FFM_FFC_invimage(x, y, pol);
4041 2093 : default: return gen_0;
4042 : }
4043 : }
4044 :
4045 : GEN
4046 2289 : RgM_RgC_invimage(GEN A, GEN y)
4047 : {
4048 2289 : pari_sp av = avma;
4049 2289 : long i, l = lg(A);
4050 : GEN M, x, t;
4051 2289 : if (l==1) return NULL;
4052 2184 : if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
4053 2184 : M = RgM_RgC_invimage_fast(A, y);
4054 2184 : if (!M) return gc_NULL(av);
4055 2163 : if (M != gen_0) return M;
4056 2093 : M = ker(shallowconcat(A, y));
4057 2093 : i = lg(M)-1;
4058 2093 : if (!i) return gc_NULL(av);
4059 :
4060 1834 : x = gel(M,i); t = gel(x,l);
4061 1834 : if (gequal0(t)) return gc_NULL(av);
4062 :
4063 1799 : t = gneg_i(t); setlg(x,l);
4064 1799 : return gerepileupto(av, RgC_Rg_div(x, t));
4065 : }
4066 :
4067 : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
4068 : * if no solution exist */
4069 : GEN
4070 2450 : inverseimage(GEN m, GEN v)
4071 : {
4072 : GEN y;
4073 2450 : if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
4074 2450 : switch(typ(v))
4075 : {
4076 2212 : case t_COL:
4077 2212 : y = RgM_RgC_invimage(m,v);
4078 2212 : return y? y: cgetg(1,t_COL);
4079 238 : case t_MAT:
4080 238 : y = RgM_invimage(m, v);
4081 238 : return y? y: cgetg(1,t_MAT);
4082 : }
4083 0 : pari_err_TYPE("inverseimage",v);
4084 : return NULL;/*LCOV_EXCL_LINE*/
4085 : }
4086 :
4087 : static GEN
4088 84 : RgM_invimage_FpM(GEN A, GEN B, GEN p)
4089 : {
4090 84 : pari_sp av = avma;
4091 : ulong pp;
4092 : GEN x;
4093 84 : A = RgM_Fp_init(A,p,&pp);
4094 84 : switch(pp)
4095 : {
4096 35 : case 0:
4097 35 : B = RgM_to_FpM(B,p);
4098 35 : x = FpM_invimage_gen(A, B, p);
4099 35 : return x ? gerepileupto(av, FpM_to_mod(x, p)): x;
4100 7 : case 2:
4101 7 : B = RgM_to_F2m(B);
4102 7 : x = F2m_invimage_i(A, B);
4103 7 : return x ? gerepileupto(av, F2m_to_mod(x)): x;
4104 42 : default:
4105 42 : B = RgM_to_Flm(B,pp);
4106 42 : x = Flm_invimage_i(A, B, pp);
4107 42 : return x ? gerepileupto(av, Flm_to_mod(x, pp)): x;
4108 : }
4109 : }
4110 :
4111 : static GEN
4112 364 : RgM_invimage_fast(GEN x, GEN y)
4113 : {
4114 : GEN p, pol;
4115 : long pa;
4116 364 : long t = RgM_type2(x, y, &p,&pol,&pa);
4117 364 : switch(t)
4118 : {
4119 84 : case t_INTMOD: return RgM_invimage_FpM(x, y, p);
4120 105 : case t_FFELT: return FFM_invimage(x, y, pol);
4121 175 : default: return gen_0;
4122 : }
4123 : }
4124 :
4125 : /* find Z such that A Z = B. Return NULL if no solution */
4126 : GEN
4127 364 : RgM_invimage(GEN A, GEN B)
4128 : {
4129 364 : pari_sp av = avma;
4130 : GEN d, x, X, Y;
4131 364 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
4132 364 : X = RgM_invimage_fast(A, B);
4133 364 : if (!X) return gc_NULL(av);
4134 252 : if (X != gen_0) return X;
4135 175 : x = ker(shallowconcat(RgM_neg(A), B));
4136 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
4137 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
4138 : * Y has at least nB columns and full rank */
4139 175 : nY = lg(x)-1;
4140 175 : if (nY < nB) return gc_NULL(av);
4141 161 : Y = rowslice(x, nA+1, nA+nB); /* nB rows */
4142 161 : d = cgetg(nB+1, t_VECSMALL);
4143 721 : for (i = nB, j = nY; i >= 1; i--, j--)
4144 : {
4145 805 : for (; j>=1; j--)
4146 756 : if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
4147 609 : if (!j) return gc_NULL(av);
4148 : }
4149 : /* reduce to the case Y square, upper triangular with 1s on diagonal */
4150 112 : Y = vecpermute(Y, d);
4151 112 : x = vecpermute(x, d);
4152 112 : X = rowslice(x, 1, nA);
4153 112 : return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
4154 : }
4155 :
4156 : static GEN
4157 70 : RgM_suppl_FpM(GEN x, GEN p)
4158 : {
4159 70 : pari_sp av = avma;
4160 : ulong pp;
4161 70 : x = RgM_Fp_init(x, p, &pp);
4162 70 : switch(pp)
4163 : {
4164 21 : case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
4165 14 : case 2: x = F2m_to_mod(F2m_suppl(x)); break;
4166 35 : default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
4167 : }
4168 70 : return gerepileupto(av, x);
4169 : }
4170 :
4171 : static GEN
4172 175 : RgM_suppl_fast(GEN x)
4173 : {
4174 : GEN p, pol;
4175 : long pa;
4176 175 : long t = RgM_type(x,&p,&pol,&pa);
4177 175 : switch(t)
4178 : {
4179 70 : case t_INTMOD: return RgM_suppl_FpM(x, p);
4180 35 : case t_FFELT: return FFM_suppl(x, pol);
4181 70 : default: return NULL;
4182 : }
4183 : }
4184 :
4185 : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
4186 : * whose first k columns are given by x. If rank(x) < k, undefined result. */
4187 : GEN
4188 175 : suppl(GEN x)
4189 : {
4190 175 : pari_sp av = avma;
4191 : GEN d, M;
4192 : long r;
4193 175 : if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
4194 175 : M = RgM_suppl_fast(x);
4195 175 : if (M) return M;
4196 70 : init_suppl(x);
4197 70 : d = gauss_pivot(x,&r);
4198 70 : set_avma(av); return get_suppl(x,d,nbrows(x),r,&col_ei);
4199 : }
4200 :
4201 : GEN
4202 7 : image2(GEN x)
4203 : {
4204 7 : pari_sp av = avma;
4205 : long k, n, i;
4206 : GEN A, B;
4207 :
4208 7 : if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
4209 7 : if (lg(x) == 1) return cgetg(1,t_MAT);
4210 7 : A = ker(x); k = lg(A)-1;
4211 7 : if (!k) { set_avma(av); return gcopy(x); }
4212 7 : A = suppl(A); n = lg(A)-1;
4213 7 : B = cgetg(n-k+1, t_MAT);
4214 21 : for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
4215 7 : return gerepileupto(av, B);
4216 : }
4217 :
4218 : GEN
4219 217 : matimage0(GEN x,long flag)
4220 : {
4221 217 : switch(flag)
4222 : {
4223 210 : case 0: return image(x);
4224 7 : case 1: return image2(x);
4225 0 : default: pari_err_FLAG("matimage");
4226 : }
4227 : return NULL; /* LCOV_EXCL_LINE */
4228 : }
4229 :
4230 : static long
4231 126 : RgM_rank_FpM(GEN x, GEN p)
4232 : {
4233 126 : pari_sp av = avma;
4234 : ulong pp;
4235 : long r;
4236 126 : x = RgM_Fp_init(x,p,&pp);
4237 126 : switch(pp)
4238 : {
4239 28 : case 0: r = FpM_rank(x,p); break;
4240 63 : case 2: r = F2m_rank(x); break;
4241 35 : default:r = Flm_rank(x,pp); break;
4242 : }
4243 126 : return gc_long(av, r);
4244 : }
4245 :
4246 : static long
4247 49 : RgM_rank_FqM(GEN x, GEN pol, GEN p)
4248 : {
4249 49 : pari_sp av = avma;
4250 : long r;
4251 49 : GEN T = RgX_to_FpX(pol, p);
4252 49 : if (signe(T) == 0) pari_err_OP("rank",x,pol);
4253 42 : r = FqM_rank(RgM_to_FqM(x, T, p), T, p);
4254 42 : return gc_long(av,r);
4255 : }
4256 :
4257 : #define code(t1,t2) ((t1 << 6) | t2)
4258 : static long
4259 315 : RgM_rank_fast(GEN x)
4260 : {
4261 : GEN p, pol;
4262 : long pa;
4263 315 : long t = RgM_type(x,&p,&pol,&pa);
4264 315 : switch(t)
4265 : {
4266 42 : case t_INT: return ZM_rank(x);
4267 21 : case t_FRAC: return QM_rank(x);
4268 126 : case t_INTMOD: return RgM_rank_FpM(x, p);
4269 70 : case t_FFELT: return FFM_rank(x, pol);
4270 49 : case code(t_POLMOD, t_INTMOD):
4271 49 : return RgM_rank_FqM(x, pol, p);
4272 7 : default: return -1;
4273 : }
4274 : }
4275 : #undef code
4276 :
4277 : long
4278 315 : rank(GEN x)
4279 : {
4280 315 : pari_sp av = avma;
4281 : long r;
4282 :
4283 315 : if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
4284 315 : r = RgM_rank_fast(x);
4285 308 : if (r >= 0) return r;
4286 7 : (void)gauss_pivot(x, &r);
4287 7 : return gc_long(av, lg(x)-1 - r);
4288 : }
4289 :
4290 : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
4291 : * followed by the missing indices */
4292 : static GEN
4293 43574 : perm_complete(GEN d, long n)
4294 : {
4295 43574 : GEN y = cgetg(n+1, t_VECSMALL);
4296 43574 : long i, j = 1, k = n, l = lg(d);
4297 43574 : pari_sp av = avma;
4298 43574 : char *T = stack_calloc(n+1);
4299 214418 : for (i = 1; i < l; i++) T[d[i]] = 1;
4300 417572 : for (i = 1; i <= n; i++)
4301 373998 : if (T[i]) y[j++] = i; else y[k--] = i;
4302 43574 : return gc_const(av, y);
4303 : }
4304 :
4305 : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
4306 : static GEN
4307 6181 : indeximage0(long n, long r, GEN d)
4308 : {
4309 : long i, j;
4310 : GEN v;
4311 :
4312 6181 : r = n - r; /* now r = dim Im(x) */
4313 6181 : v = cgetg(r+1,t_VECSMALL);
4314 34419 : if (d) for (i=j=1; j<=n; j++)
4315 28238 : if (d[j]) v[i++] = j;
4316 6181 : return v;
4317 : }
4318 : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from gauss_pivot */
4319 : static void
4320 21787 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
4321 : {
4322 21787 : GEN IR = indexrank0(n, r, d);
4323 21787 : *prow = perm_complete(gel(IR,1), m);
4324 21787 : *pcol = perm_complete(gel(IR,2), n);
4325 21787 : }
4326 :
4327 : static GEN
4328 28 : RgM_indexrank_FpM(GEN x, GEN p)
4329 : {
4330 28 : pari_sp av = avma;
4331 : ulong pp;
4332 : GEN r;
4333 28 : x = RgM_Fp_init(x,p,&pp);
4334 28 : switch(pp)
4335 : {
4336 7 : case 0: r = FpM_indexrank(x,p); break;
4337 7 : case 2: r = F2m_indexrank(x); break;
4338 14 : default: r = Flm_indexrank(x,pp); break;
4339 : }
4340 28 : return gerepileupto(av, r);
4341 : }
4342 :
4343 : static GEN
4344 0 : RgM_indexrank_FqM(GEN x, GEN pol, GEN p)
4345 : {
4346 0 : pari_sp av = avma;
4347 0 : GEN r, T = RgX_to_FpX(pol, p);
4348 0 : if (signe(T) == 0) pari_err_OP("indexrank",x,pol);
4349 0 : r = FqM_indexrank(RgM_to_FqM(x, T, p), T, p);
4350 0 : return gerepileupto(av, r);
4351 : }
4352 :
4353 : #define code(t1,t2) ((t1 << 6) | t2)
4354 : static GEN
4355 77512 : RgM_indexrank_fast(GEN x)
4356 : {
4357 : GEN p, pol;
4358 : long pa;
4359 77512 : long t = RgM_type(x,&p,&pol,&pa);
4360 77512 : switch(t)
4361 : {
4362 406 : case t_INT: return ZM_indexrank(x);
4363 1344 : case t_FRAC: return QM_indexrank(x);
4364 28 : case t_INTMOD: return RgM_indexrank_FpM(x, p);
4365 21 : case t_FFELT: return FFM_indexrank(x, pol);
4366 0 : case code(t_POLMOD, t_INTMOD):
4367 0 : return RgM_indexrank_FqM(x, pol, p);
4368 75713 : default: return NULL;
4369 : }
4370 : }
4371 : #undef code
4372 :
4373 : GEN
4374 77512 : indexrank(GEN x)
4375 : {
4376 : pari_sp av;
4377 : long r;
4378 : GEN d;
4379 77512 : if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
4380 77512 : d = RgM_indexrank_fast(x);
4381 77512 : if (d) return d;
4382 75713 : av = avma;
4383 75713 : init_indexrank(x);
4384 75713 : d = gauss_pivot(x, &r);
4385 75714 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4386 : }
4387 :
4388 : GEN
4389 6181 : ZM_indeximage(GEN x) {
4390 6181 : pari_sp av = avma;
4391 : long r;
4392 : GEN d;
4393 6181 : init_indexrank(x);
4394 6181 : d = ZM_pivots(x,&r);
4395 6181 : set_avma(av); return indeximage0(lg(x)-1, r, d);
4396 : }
4397 : long
4398 2224635 : ZM_rank(GEN x) {
4399 2224635 : pari_sp av = avma;
4400 : long r;
4401 2224635 : (void)ZM_pivots(x,&r);
4402 2224633 : return gc_long(av, lg(x)-1-r);
4403 : }
4404 : GEN
4405 1773569 : ZM_indexrank(GEN x) {
4406 1773569 : pari_sp av = avma;
4407 : long r;
4408 : GEN d;
4409 1773569 : init_indexrank(x);
4410 1773568 : d = ZM_pivots(x,&r);
4411 1773571 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4412 : }
4413 :
4414 : long
4415 21 : QM_rank(GEN x)
4416 : {
4417 21 : pari_sp av = avma;
4418 21 : long r = ZM_rank(Q_primpart(x));
4419 21 : set_avma(av);
4420 21 : return r;
4421 : }
4422 :
4423 : GEN
4424 1344 : QM_indexrank(GEN x)
4425 : {
4426 1344 : pari_sp av = avma;
4427 1344 : GEN r = ZM_indexrank(Q_primpart(x));
4428 1344 : return gerepileupto(av, r);
4429 : }
4430 :
4431 : /*******************************************************************/
4432 : /* */
4433 : /* ZabM */
4434 : /* */
4435 : /*******************************************************************/
4436 :
4437 : static GEN
4438 1276 : FpXM_ratlift(GEN a, GEN q)
4439 : {
4440 : GEN B, y;
4441 1276 : long i, j, l = lg(a), n;
4442 1276 : B = sqrti(shifti(q,-1));
4443 1276 : y = cgetg(l, t_MAT);
4444 1276 : if (l==1) return y;
4445 1276 : n = lgcols(a);
4446 3059 : for (i=1; i<l; i++)
4447 : {
4448 2404 : GEN yi = cgetg(n, t_COL);
4449 32311 : for (j=1; j<n; j++)
4450 : {
4451 30528 : GEN v = FpX_ratlift(gmael(a,i,j), q, B, B, NULL);
4452 30528 : if (!v) return NULL;
4453 29907 : gel(yi, j) = RgX_renormalize(v);
4454 : }
4455 1783 : gel(y,i) = yi;
4456 : }
4457 655 : return y;
4458 : }
4459 :
4460 : static GEN
4461 4485 : FlmV_recover_pre(GEN a, GEN M, ulong p, ulong pi, long sv)
4462 : {
4463 4485 : GEN a1 = gel(a,1);
4464 4485 : long i, j, k, l = lg(a1), n, lM = lg(M);
4465 4485 : GEN v = cgetg(lM, t_VECSMALL);
4466 4485 : GEN y = cgetg(l, t_MAT);
4467 4485 : if (l==1) return y;
4468 4485 : n = lgcols(a1);
4469 22521 : for (i=1; i<l; i++)
4470 : {
4471 18036 : GEN yi = cgetg(n, t_COL);
4472 347607 : for (j=1; j<n; j++)
4473 : {
4474 4676397 : for (k=1; k<lM; k++) uel(v,k) = umael(gel(a,k),i,j);
4475 329571 : gel(yi, j) = Flm_Flc_mul_pre_Flx(M, v, p, pi, sv);
4476 : }
4477 18036 : gel(y,i) = yi;
4478 : }
4479 4485 : return y;
4480 : }
4481 :
4482 : static GEN
4483 0 : FlkM_inv(GEN M, GEN P, ulong p)
4484 : {
4485 0 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4486 0 : GEN R = Flx_roots_pre(P, p, pi);
4487 0 : long l = lg(R), i;
4488 0 : GEN W = Flv_invVandermonde(R, 1UL, p);
4489 0 : GEN V = cgetg(l, t_VEC);
4490 0 : for(i=1; i<l; i++)
4491 : {
4492 0 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, PI);
4493 0 : GEN H = Flm_inv_sp(FlxM_eval_powers_pre(M, pows, p, pi), NULL, p);
4494 0 : if (!H) return NULL;
4495 0 : gel(V, i) = H;
4496 : }
4497 0 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4498 : }
4499 :
4500 : static GEN
4501 3209 : FlkM_adjoint(GEN M, GEN P, ulong p)
4502 : {
4503 3209 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4504 3209 : GEN R = Flx_roots_pre(P, p, pi);
4505 3209 : long l = lg(R), i;
4506 3209 : GEN W = Flv_invVandermonde(R, 1UL, p);
4507 3209 : GEN V = cgetg(l, t_VEC);
4508 15577 : for(i=1; i<l; i++)
4509 : {
4510 12368 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, PI);
4511 12368 : gel(V, i) = Flm_adjoint(FlxM_eval_powers_pre(M, pows, p, pi), p);
4512 : }
4513 3209 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4514 : }
4515 :
4516 : static GEN
4517 1985 : ZabM_inv_slice(GEN A, GEN Q, GEN P, GEN *mod)
4518 : {
4519 1985 : pari_sp av = avma;
4520 1985 : long i, n = lg(P)-1, w = varn(Q);
4521 : GEN H, T;
4522 1985 : if (n == 1)
4523 : {
4524 1554 : ulong p = uel(P,1);
4525 1554 : GEN Qp = ZX_to_Flx(Q, p);
4526 1554 : GEN Ap = ZXM_to_FlxM(A, p, get_Flx_var(Qp));
4527 1554 : GEN Hp = FlkM_adjoint(Ap, Qp, p);
4528 1554 : Hp = gerepileupto(av, FlxM_to_ZXM(Hp));
4529 1554 : *mod = utoipos(p); return Hp;
4530 : }
4531 431 : T = ZV_producttree(P);
4532 431 : A = ZXM_nv_mod_tree(A, P, T, w);
4533 431 : Q = ZX_nv_mod_tree(Q, P, T);
4534 431 : H = cgetg(n+1, t_VEC);
4535 2086 : for(i=1; i <= n; i++)
4536 : {
4537 1655 : ulong p = P[i];
4538 1655 : GEN a = gel(A,i), q = gel(Q, i);
4539 1655 : gel(H,i) = FlkM_adjoint(a, q, p);
4540 : }
4541 431 : H = nxMV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
4542 431 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
4543 : }
4544 :
4545 : GEN
4546 1985 : ZabM_inv_worker(GEN P, GEN A, GEN Q)
4547 : {
4548 1985 : GEN V = cgetg(3, t_VEC);
4549 1985 : gel(V,1) = ZabM_inv_slice(A, Q, P, &gel(V,2));
4550 1985 : return V;
4551 : }
4552 :
4553 : static GEN
4554 5509 : vecnorml1(GEN x)
4555 60508 : { pari_APPLY_same(gnorml1_fake(gel(x,i))); }
4556 :
4557 : static GEN
4558 1827 : ZabM_true_Hadamard(GEN a)
4559 : {
4560 1827 : pari_sp av = avma;
4561 1827 : long n = lg(a)-1, i;
4562 : GEN B;
4563 1827 : if (n == 0) return gen_1;
4564 1827 : if (n == 1) return gnorml1_fake(gcoeff(a,1,1));
4565 1183 : B = gen_1;
4566 6692 : for (i = 1; i <= n; i++)
4567 5509 : B = gmul(B, gnorml2(RgC_gtofp(vecnorml1(gel(a,i)),DEFAULTPREC)));
4568 1183 : return gerepileuptoint(av, ceil_safe(sqrtr_abs(B)));
4569 : }
4570 :
4571 : GEN
4572 1827 : ZabM_inv(GEN A, GEN Q, long n, GEN *pt_den)
4573 : {
4574 1827 : pari_sp av = avma;
4575 : forprime_t S;
4576 : GEN bnd, H, D, d, mod, worker;
4577 1827 : if (lg(A) == 1)
4578 : {
4579 0 : if (pt_den) *pt_den = gen_1;
4580 0 : return cgetg(1, t_MAT);
4581 : }
4582 1827 : bnd = ZabM_true_Hadamard(A);
4583 1827 : worker = snm_closure(is_entry("_ZabM_inv_worker"), mkvec2(A, Q));
4584 1827 : u_forprime_arith_init(&S, HIGHBIT+1, ULONG_MAX, 1, n);
4585 1827 : H = gen_crt("ZabM_inv", worker, &S, NULL, expi(bnd), 0, &mod,
4586 : nxMV_chinese_center, FpXM_center);
4587 1827 : D = RgMrow_RgC_mul(H, gel(A,1), 1);
4588 1827 : D = ZX_rem(D, Q);
4589 1827 : d = Z_content(mkvec2(H, D));
4590 1827 : if (d)
4591 : {
4592 518 : D = ZX_Z_divexact(D, d);
4593 518 : H = Q_div_to_int(H, d);
4594 : }
4595 1827 : if (!pt_den) return gerepileupto(av, H);
4596 1827 : *pt_den = D; return gc_all(av, 2, &H, pt_den);
4597 : }
4598 :
4599 : GEN
4600 0 : ZabM_inv_ratlift(GEN M, GEN P, long n, GEN *pden)
4601 : {
4602 0 : pari_sp av2, av = avma;
4603 : GEN q, H;
4604 0 : ulong m = LONG_MAX>>1;
4605 0 : ulong p= 1 + m - (m % n);
4606 0 : long lM = lg(M);
4607 0 : if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
4608 :
4609 0 : av2 = avma;
4610 0 : H = NULL;
4611 : for(;;)
4612 0 : {
4613 : GEN Hp, Pp, Mp, Hr;
4614 0 : do p += n; while(!uisprime(p));
4615 0 : Pp = ZX_to_Flx(P, p);
4616 0 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4617 0 : Hp = FlkM_inv(Mp, Pp, p);
4618 0 : if (!Hp) continue;
4619 0 : if (!H)
4620 : {
4621 0 : H = ZXM_init_CRT(Hp, degpol(P)-1, p);
4622 0 : q = utoipos(p);
4623 : }
4624 : else
4625 0 : ZXM_incremental_CRT(&H, Hp, &q, p);
4626 0 : Hr = FpXM_ratlift(H, q);
4627 0 : if (DEBUGLEVEL>5) err_printf("ZabM_inv mod %ld (ratlift=%ld)\n", p,!!Hr);
4628 0 : if (Hr) {/* DONE ? */
4629 0 : GEN Hl = Q_remove_denom(Hr, pden);
4630 0 : GEN MH = ZXQM_mul(Hl, M, P);
4631 0 : if (*pden)
4632 0 : { if (RgM_isscalar(MH, *pden)) { H = Hl; break; }}
4633 : else
4634 0 : { if (RgM_isidentity(MH)) { H = Hl; *pden = gen_1; break; } }
4635 : }
4636 :
4637 0 : if (gc_needed(av,2))
4638 : {
4639 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_inv");
4640 0 : gerepileall(av2, 2, &H, &q);
4641 : }
4642 : }
4643 0 : return gc_all(av, 2, &H, pden);
4644 : }
4645 :
4646 : static GEN
4647 1276 : FlkM_ker(GEN M, GEN P, ulong p)
4648 : {
4649 1276 : ulong PI = get_Fl_red(p), pi = SMALL_ULONG(p)? 0: PI;
4650 1276 : GEN R = Flx_roots_pre(P, p, pi);
4651 1276 : long l = lg(R), i, dP = degpol(P), r;
4652 : GEN M1, K, D;
4653 1276 : GEN W = Flv_invVandermonde(R, 1UL, p);
4654 1276 : GEN V = cgetg(l, t_VEC);
4655 1276 : M1 = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,1), dP, p, PI), p, pi);
4656 1276 : K = Flm_ker_sp(M1, p, 2);
4657 1276 : r = lg(gel(K,1)); D = gel(K,2);
4658 1276 : gel(V, 1) = gel(K,1);
4659 2652 : for(i=2; i<l; i++)
4660 : {
4661 1376 : GEN Mi = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,i), dP, p, PI), p, pi);
4662 1376 : GEN K = Flm_ker_sp(Mi, p, 2);
4663 1376 : if (lg(gel(K,1)) != r || !zv_equal(D, gel(K,2))) return NULL;
4664 1376 : gel(V, i) = gel(K,1);
4665 : }
4666 1276 : return mkvec2(FlmV_recover_pre(V, W, p, pi, P[1]), D);
4667 : }
4668 :
4669 : static int
4670 655 : ZabM_ker_check(GEN M, GEN H, ulong p, GEN P, long n)
4671 : {
4672 : GEN pow;
4673 655 : long j, l = lg(H);
4674 : ulong pi, r;
4675 3899 : do p += n; while(!uisprime(p));
4676 655 : pi = get_Fl_red(p);
4677 655 : P = ZX_to_Flx(P, p);
4678 655 : r = Flx_oneroot_pre(P, p, pi);
4679 655 : pow = Fl_powers_pre(r, degpol(P),p, (p & HIGHMASK)? pi: 0);
4680 655 : M = ZXM_to_FlxM(M, p, P[1]); M = FlxM_eval_powers_pre(M, pow, p, pi);
4681 655 : H = ZXM_to_FlxM(H, p, P[1]); H = FlxM_eval_powers_pre(H, pow, p, pi);
4682 2178 : for (j = 1; j < l; j++)
4683 1555 : if (!zv_equal0(Flm_Flc_mul_pre(M, gel(H,j), p, pi))) return 0;
4684 623 : return 1;
4685 : }
4686 :
4687 : GEN
4688 623 : ZabM_ker(GEN M, GEN P, long n)
4689 : {
4690 623 : pari_sp av = avma;
4691 : pari_timer ti;
4692 623 : GEN q, H = NULL, D = NULL;
4693 623 : ulong m = LONG_MAX>>1;
4694 623 : ulong p = 1 + m - (m % n);
4695 :
4696 623 : if (DEBUGLEVEL>5) timer_start(&ti);
4697 : for(;;)
4698 653 : {
4699 : GEN Kp, Hp, Dp, Pp, Mp, Hr;
4700 22341 : do p += n; while(!uisprime(p));
4701 1276 : Pp = ZX_to_Flx(P, p);
4702 1276 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4703 1276 : Kp = FlkM_ker(Mp, Pp, p);
4704 1276 : if (!Kp) continue;
4705 1276 : Hp = gel(Kp,1); Dp = gel(Kp,2);
4706 1276 : if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
4707 1276 : if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
4708 : {
4709 623 : H = ZXM_init_CRT(Hp, degpol(P)-1, p); D = Dp;
4710 623 : q = utoipos(p);
4711 : }
4712 : else
4713 653 : ZXM_incremental_CRT(&H, Hp, &q, p);
4714 1276 : Hr = FpXM_ratlift(H, q);
4715 1276 : if (DEBUGLEVEL>5) timer_printf(&ti,"ZabM_ker mod %ld (ratlift=%ld)", p,!!Hr);
4716 1276 : if (Hr) {/* DONE ? */
4717 655 : GEN Hl = vec_Q_primpart(Hr);
4718 655 : if (ZabM_ker_check(M, Hl, p, P, n)) { H = Hl; break; }
4719 : }
4720 :
4721 653 : if (gc_needed(av,2))
4722 : {
4723 4 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_ker");
4724 4 : gerepileall(av, 3, &H, &D, &q);
4725 : }
4726 : }
4727 623 : return gerepilecopy(av, H);
4728 : }
4729 :
4730 : GEN
4731 2387 : ZabM_indexrank(GEN M, GEN P, long n)
4732 : {
4733 2387 : pari_sp av = avma;
4734 2387 : ulong m = LONG_MAX>>1;
4735 2387 : ulong p = 1+m-(m%n), D = degpol(P);
4736 2387 : long lM = lg(M), lmax = 0, c = 0;
4737 : GEN v;
4738 : for(;;)
4739 735 : {
4740 : GEN R, Pp, Mp, K;
4741 : ulong pi;
4742 : long l;
4743 61415 : do p += n; while (!uisprime(p));
4744 3122 : pi = (p & HIGHMASK)? get_Fl_red(p): 0;
4745 3122 : Pp = ZX_to_Flx(P, p);
4746 3122 : R = Flx_roots_pre(Pp, p, pi);
4747 3122 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4748 3122 : K = FlxM_eval_powers_pre(Mp, Fl_powers_pre(uel(R,1), D,p,pi), p,pi);
4749 3122 : v = Flm_indexrank(K, p);
4750 3122 : l = lg(gel(v,2));
4751 3122 : if (l == lM) break;
4752 980 : if (lmax >= 0 && l > lmax) { lmax = l; c = 0; } else c++;
4753 980 : if (c > 2)
4754 : { /* probably not maximal rank, expensive check */
4755 245 : lM -= lg(ZabM_ker(M, P, n))-1; /* actual rank (+1) */
4756 245 : if (lmax == lM) break;
4757 0 : lmax = -1; /* disable check */
4758 : }
4759 : }
4760 2387 : return gerepileupto(av, v);
4761 : }
4762 :
4763 : #if 0
4764 : GEN
4765 : ZabM_gauss(GEN M, GEN P, long n, GEN *den)
4766 : {
4767 : pari_sp av = avma;
4768 : GEN v, S, W;
4769 : v = ZabM_indexrank(M, P, n);
4770 : S = shallowmatextract(M,gel(v,1),gel(v,2));
4771 : W = ZabM_inv(S, P, n, den);
4772 : return gc_all(av,2,&W,den);
4773 : }
4774 : #endif
4775 :
4776 : GEN
4777 140 : ZabM_pseudoinv(GEN M, GEN P, long n, GEN *pv, GEN *den)
4778 : {
4779 140 : GEN v = ZabM_indexrank(M, P, n);
4780 140 : if (pv) *pv = v;
4781 140 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4782 140 : return ZabM_inv(M, P, n, den);
4783 : }
4784 : GEN
4785 5019 : ZM_pseudoinv(GEN M, GEN *pv, GEN *den)
4786 : {
4787 5019 : GEN v = ZM_indexrank(M);
4788 5019 : if (pv) *pv = v;
4789 5019 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4790 5019 : return ZM_inv(M, den);
4791 : }
4792 :
4793 : /*******************************************************************/
4794 : /* */
4795 : /* Structured Elimination */
4796 : /* */
4797 : /*******************************************************************/
4798 :
4799 : static void
4800 95967 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
4801 : {
4802 95967 : long lc = lg(c), k;
4803 95967 : iscol[i] = 0; (*rcol)--;
4804 891429 : for (k = 1; k < lc; ++k)
4805 : {
4806 795462 : Wrow[c[k]]--;
4807 795462 : if (Wrow[c[k]]==0) (*rrow)--;
4808 : }
4809 95967 : }
4810 :
4811 : static void
4812 7641 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long idx, long *rcol, long *rrow)
4813 : {
4814 : long i, j;
4815 7641 : long nbcol = lg(iscol)-1, last;
4816 : do
4817 : {
4818 9568 : last = 0;
4819 16915472 : for (i = 1; i <= nbcol; ++i)
4820 16905904 : if (iscol[i])
4821 : {
4822 9074581 : GEN c = idx ? gmael(M, i, idx): gel(M,i);
4823 9074581 : long lc = lg(c);
4824 83831551 : for (j = 1; j < lc; ++j)
4825 74775020 : if (Wrow[c[j]] == 1)
4826 : {
4827 18050 : rem_col(c, i, iscol, Wrow, rcol, rrow);
4828 18050 : last=1; break;
4829 : }
4830 : }
4831 9568 : } while (last);
4832 7641 : }
4833 :
4834 : static GEN
4835 7448 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
4836 : {
4837 7448 : long nbcol = lg(iscol)-1;
4838 : long i, j, m, last;
4839 : GEN per;
4840 20552 : for (m = 2, last=0; !last ; m++)
4841 : {
4842 25077619 : for (i = 1; i <= nbcol; ++i)
4843 : {
4844 25064515 : wcol[i] = 0;
4845 25064515 : if (iscol[i])
4846 : {
4847 13863667 : GEN c = gmael(M, i, 1);
4848 13863667 : long lc = lg(c);
4849 123887839 : for (j = 1; j < lc; ++j)
4850 110024172 : if (Wrow[c[j]] == m) { wcol[i]++; last = 1; }
4851 : }
4852 : }
4853 : }
4854 7448 : per = vecsmall_indexsort(wcol);
4855 7448 : *w = wcol[per[nbcol]];
4856 7448 : return per;
4857 : }
4858 :
4859 : /* M is a RgMs with nbrow rows, A a list of row indices.
4860 : Eliminate rows of M with a single entry that do not belong to A,
4861 : and the corresponding columns. Also eliminate columns until #colums=#rows.
4862 : Return pcol and prow:
4863 : pcol is a map from the new columns indices to the old one.
4864 : prow is a map from the old rows indices to the new one (0 if removed).
4865 : */
4866 :
4867 : void
4868 147 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4869 : {
4870 147 : long i, j, k, lA = lg(A);
4871 147 : GEN prow = cgetg(nbrow+1, t_VECSMALL);
4872 147 : GEN pcol = zero_zv(nbcol);
4873 147 : pari_sp av = avma;
4874 147 : long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
4875 147 : GEN iscol = const_vecsmall(nbcol, 1);
4876 147 : GEN Wrow = zero_zv(nbrow);
4877 147 : GEN wcol = cgetg(nbcol+1, t_VECSMALL);
4878 147 : pari_sp av2 = avma;
4879 110397 : for (i = 1; i <= nbcol; ++i)
4880 : {
4881 110250 : GEN F = gmael(M, i, 1);
4882 110250 : long l = lg(F)-1;
4883 924681 : for (j = 1; j <= l; ++j) Wrow[F[j]]++;
4884 : }
4885 147 : for (j = 1; j < lA; ++j)
4886 : {
4887 0 : if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
4888 0 : Wrow[A[j]] = -1;
4889 : }
4890 228354 : for (i = 1; i <= nbrow; ++i)
4891 228207 : if (Wrow[i]) rrow++;
4892 147 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow);
4893 147 : if (rcol < rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
4894 7595 : while (rcol > rrow)
4895 : {
4896 : long w;
4897 7448 : GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
4898 85365 : for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
4899 77917 : rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
4900 7448 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow); set_avma(av2);
4901 : }
4902 110397 : for (j = 1, i = 1; i <= nbcol; ++i)
4903 110250 : if (iscol[i]) pcol[j++] = i;
4904 147 : setlg(pcol,j);
4905 228354 : for (k = 1, i = 1; i <= nbrow; ++i) prow[i] = Wrow[i]? k++: 0;
4906 147 : *p_col = pcol; *p_row = prow; set_avma(av);
4907 : }
4908 :
4909 : void
4910 0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4911 0 : { RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row); }
4912 :
4913 : GEN
4914 46 : F2Ms_colelim(GEN M, long nbrow)
4915 : {
4916 46 : long i,j, nbcol = lg(M)-1, rcol = nbcol, rrow = 0;
4917 46 : GEN pcol = zero_zv(nbcol);
4918 46 : pari_sp av = avma;
4919 46 : GEN iscol = const_vecsmall(nbcol, 1), Wrow = zero_zv(nbrow);
4920 77470 : for (i = 1; i <= nbcol; ++i)
4921 : {
4922 77424 : GEN F = gel(M, i);
4923 77424 : long l = lg(F)-1;
4924 1431938 : for (j = 1; j <= l; ++j) Wrow[F[j]]++;
4925 : }
4926 46 : rem_singleton(M, iscol, Wrow, 0, &rcol, &rrow);
4927 77470 : for (j = 1, i = 1; i <= nbcol; ++i)
4928 77424 : if (iscol[i]) pcol[j++] = i;
4929 46 : fixlg(pcol,j); return gc_const(av, pcol);
4930 : }
4931 :
4932 : /*******************************************************************/
4933 : /* */
4934 : /* EIGENVECTORS */
4935 : /* (independent eigenvectors, sorted by increasing eigenvalue) */
4936 : /* */
4937 : /*******************************************************************/
4938 : /* assume x is square of dimension > 0 */
4939 : static int
4940 53 : RgM_is_symmetric_cx(GEN x, long bit)
4941 : {
4942 53 : pari_sp av = avma;
4943 53 : long i, j, l = lg(x);
4944 239 : for (i = 1; i < l; i++)
4945 708 : for (j = 1; j < i; j++)
4946 : {
4947 522 : GEN a = gcoeff(x,i,j), b = gcoeff(x,j,i), c = gsub(a,b);
4948 522 : if (!gequal0(c) && gexpo(c) - gexpo(a) > -bit) return gc_long(av,0);
4949 : }
4950 21 : return gc_long(av,1);
4951 : }
4952 : static GEN
4953 53 : eigen_err(int exact, GEN x, long flag, long prec)
4954 : {
4955 53 : pari_sp av = avma;
4956 : GEN y;
4957 53 : if (RgM_is_symmetric_cx(x, prec - 10))
4958 : { /* approximately symmetric: recover */
4959 21 : x = jacobi(x, prec); if (flag) return x;
4960 14 : return gerepilecopy(av, gel(x,2));
4961 : }
4962 32 : if (!exact) x = bestappr(x, NULL);
4963 32 : y = mateigen(x, flag, precdbl(prec));
4964 32 : if (exact)
4965 18 : y = gprec_wtrunc(y, prec);
4966 14 : else if (flag)
4967 7 : y = mkvec2(RgV_gtofp(gel(y,1), prec), RgM_gtofp(gel(y,2), prec));
4968 : else
4969 7 : y = RgM_gtofp(y, prec);
4970 32 : return gerepilecopy(av, y);
4971 : }
4972 : GEN
4973 144 : mateigen(GEN x, long flag, long prec)
4974 : {
4975 : GEN y, R, T;
4976 144 : long k, l, ex, n = lg(x);
4977 : int exact;
4978 144 : pari_sp av = avma;
4979 :
4980 144 : if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
4981 144 : if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
4982 144 : if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
4983 144 : if (n == 1)
4984 : {
4985 14 : if (flag) retmkvec2(cgetg(1,t_COL), cgetg(1,t_MAT));
4986 7 : return cgetg(1,t_MAT);
4987 : }
4988 130 : if (n == 2)
4989 : {
4990 14 : if (flag) retmkvec2(mkcolcopy(gcoeff(x,1,1)), matid(1));
4991 7 : return matid(1);
4992 : }
4993 :
4994 116 : ex = 16 - prec;
4995 116 : T = charpoly(x,0);
4996 116 : exact = RgX_is_QX(T);
4997 116 : if (exact)
4998 : {
4999 74 : T = ZX_radical( Q_primpart(T) );
5000 74 : R = nfrootsQ(T); settyp(R, t_COL);
5001 74 : if (lg(R)-1 < degpol(T))
5002 : { /* add missing complex roots */
5003 60 : GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
5004 60 : R = shallowconcat(R, r);
5005 : }
5006 : }
5007 : else
5008 : {
5009 42 : GEN r1, v = vectrunc_init(lg(T));
5010 : long e;
5011 42 : R = cleanroots(T,prec);
5012 42 : r1 = NULL;
5013 266 : for (k = 1; k < lg(R); k++)
5014 : {
5015 224 : GEN r2 = gel(R,k), r = grndtoi(r2, &e);
5016 224 : if (e < ex) r2 = r;
5017 224 : if (r1)
5018 : {
5019 182 : r = gsub(r1,r2);
5020 182 : if (gequal0(r) || gexpo(r) < ex) continue;
5021 : }
5022 182 : vectrunc_append(v, r2);
5023 182 : r1 = r2;
5024 : }
5025 42 : R = v;
5026 : }
5027 : /* R = distinct complex roots of charpoly(x) */
5028 116 : l = lg(R); y = cgetg(l, t_VEC);
5029 452 : for (k = 1; k < l; k++)
5030 : {
5031 389 : GEN F = ker_aux(RgM_Rg_sub_shallow(x, gel(R,k)), x);
5032 389 : long d = lg(F)-1;
5033 389 : if (!d) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5034 336 : gel(y,k) = F;
5035 336 : if (flag) gel(R,k) = const_col(d, gel(R,k));
5036 : }
5037 63 : y = shallowconcat1(y);
5038 63 : if (lg(y) > n) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5039 : /* lg(y) < n if x is not diagonalizable */
5040 63 : if (flag) y = mkvec2(shallowconcat1(R), y);
5041 63 : return gerepilecopy(av,y);
5042 : }
5043 : GEN
5044 0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
5045 :
5046 : /*******************************************************************/
5047 : /* */
5048 : /* DETERMINANT */
5049 : /* */
5050 : /*******************************************************************/
5051 :
5052 : GEN
5053 26593 : det0(GEN a,long flag)
5054 : {
5055 26593 : switch(flag)
5056 : {
5057 26579 : case 0: return det(a);
5058 14 : case 1: return det2(a);
5059 0 : default: pari_err_FLAG("matdet");
5060 : }
5061 : return NULL; /* LCOV_EXCL_LINE */
5062 : }
5063 :
5064 : /* M a 2x2 matrix, returns det(M) */
5065 : static GEN
5066 94256 : RgM_det2(GEN M)
5067 : {
5068 94256 : pari_sp av = avma;
5069 94256 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5070 94256 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5071 94256 : return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
5072 : }
5073 : /* M a 2x2 ZM, returns det(M) */
5074 : static GEN
5075 8659 : ZM_det2(GEN M)
5076 : {
5077 8659 : pari_sp av = avma;
5078 8659 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5079 8659 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5080 8659 : return gerepileuptoint(av, subii(mulii(a,d), mulii(b, c)));
5081 : }
5082 : /* M a 3x3 ZM, return det(M) */
5083 : static GEN
5084 100472 : ZM_det3(GEN M)
5085 : {
5086 100472 : pari_sp av = avma;
5087 100472 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
5088 100472 : GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
5089 100472 : GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
5090 100472 : GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
5091 100472 : if (signe(g))
5092 : {
5093 66202 : t = mulii(subii(mulii(b,f), mulii(c,e)), g);
5094 66202 : D = addii(D, t);
5095 : }
5096 100472 : if (signe(h))
5097 : {
5098 77604 : t = mulii(subii(mulii(c,d), mulii(a,f)), h);
5099 77604 : D = addii(D, t);
5100 : }
5101 100472 : return gerepileuptoint(av, D);
5102 : }
5103 :
5104 : static GEN
5105 58226 : det_simple_gauss(GEN a, GEN data, pivot_fun pivot)
5106 : {
5107 58226 : pari_sp av = avma;
5108 58226 : long i,j,k, s = 1, nbco = lg(a)-1;
5109 58226 : GEN p, x = gen_1;
5110 :
5111 58226 : a = RgM_shallowcopy(a);
5112 342023 : for (i=1; i<nbco; i++)
5113 : {
5114 283805 : k = pivot(a, data, i, NULL);
5115 283808 : if (k > nbco) return gerepilecopy(av, gcoeff(a,i,i));
5116 283801 : if (k != i)
5117 : { /* exchange the lines s.t. k = i */
5118 1159671 : for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
5119 119171 : s = -s;
5120 : }
5121 283801 : p = gcoeff(a,i,i);
5122 :
5123 283801 : x = gmul(x,p);
5124 1786859 : for (k=i+1; k<=nbco; k++)
5125 : {
5126 1503063 : GEN m = gcoeff(a,i,k);
5127 1503063 : if (gequal0(m)) continue;
5128 :
5129 1067758 : m = gdiv(m,p);
5130 9952421 : for (j=i+1; j<=nbco; j++)
5131 8884665 : gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
5132 : }
5133 283796 : if (gc_needed(av,2))
5134 : {
5135 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5136 0 : gerepileall(av,2, &a,&x);
5137 : }
5138 : }
5139 58218 : if (s < 0) x = gneg_i(x);
5140 58218 : return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
5141 : }
5142 :
5143 : GEN
5144 134468 : det2(GEN a)
5145 : {
5146 : GEN data;
5147 : pivot_fun pivot;
5148 134468 : long n = lg(a)-1;
5149 134468 : if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
5150 134468 : if (!n) return gen_1;
5151 134468 : if (n != nbrows(a)) pari_err_DIM("det2");
5152 134468 : if (n == 1) return gcopy(gcoeff(a,1,1));
5153 85795 : if (n == 2) return RgM_det2(a);
5154 26706 : pivot = get_pivot_fun(a, a, &data);
5155 26706 : return det_simple_gauss(a, data, pivot);
5156 : }
5157 :
5158 : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
5159 : * Gauss-Bareiss. */
5160 : static GEN
5161 462 : det_bareiss(GEN a)
5162 : {
5163 462 : pari_sp av = avma;
5164 462 : long nbco = lg(a)-1,i,j,k,s = 1;
5165 : GEN p, pprec;
5166 :
5167 462 : a = RgM_shallowcopy(a);
5168 1337 : for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
5169 : {
5170 882 : int diveuc = (gequal1(pprec)==0);
5171 : GEN ci;
5172 :
5173 882 : p = gcoeff(a,i,i);
5174 882 : if (gequal0(p))
5175 : {
5176 14 : k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
5177 7 : if (k>nbco) return gerepilecopy(av, p);
5178 0 : swap(gel(a,k), gel(a,i)); s = -s;
5179 0 : p = gcoeff(a,i,i);
5180 : }
5181 875 : ci = gel(a,i);
5182 2373 : for (k=i+1; k<=nbco; k++)
5183 : {
5184 1498 : GEN ck = gel(a,k), m = gel(ck,i);
5185 1498 : if (gequal0(m))
5186 : {
5187 7 : if (gequal1(p))
5188 : {
5189 0 : if (diveuc)
5190 0 : gel(a,k) = gdiv(gel(a,k), pprec);
5191 : }
5192 : else
5193 42 : for (j=i+1; j<=nbco; j++)
5194 : {
5195 35 : GEN p1 = gmul(p, gel(ck,j));
5196 35 : if (diveuc) p1 = gdiv(p1,pprec);
5197 35 : gel(ck,j) = p1;
5198 : }
5199 : }
5200 : else
5201 4662 : for (j=i+1; j<=nbco; j++)
5202 : {
5203 3171 : pari_sp av2 = avma;
5204 3171 : GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
5205 3171 : if (diveuc) p1 = gdiv(p1,pprec);
5206 3171 : gel(ck,j) = gerepileupto(av2, p1);
5207 : }
5208 1498 : if (gc_needed(av,2))
5209 : {
5210 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5211 0 : gerepileall(av,2, &a,&pprec);
5212 0 : ci = gel(a,i);
5213 0 : p = gcoeff(a,i,i);
5214 : }
5215 : }
5216 : }
5217 455 : p = gcoeff(a,nbco,nbco);
5218 455 : p = (s < 0)? gneg(p): gcopy(p);
5219 455 : return gerepileupto(av, p);
5220 : }
5221 :
5222 : /* count nonzero entries in col j, at most 'max' of them.
5223 : * Return their indices */
5224 : static GEN
5225 1470 : col_count_non_zero(GEN a, long j, long max)
5226 : {
5227 1470 : GEN v = cgetg(max+1, t_VECSMALL);
5228 1470 : GEN c = gel(a,j);
5229 1470 : long i, l = lg(a), k = 1;
5230 5614 : for (i = 1; i < l; i++)
5231 5376 : if (!gequal0(gel(c,i)))
5232 : {
5233 5110 : if (k > max) return NULL; /* fail */
5234 3878 : v[k++] = i;
5235 : }
5236 238 : setlg(v, k); return v;
5237 : }
5238 : /* count nonzero entries in row i, at most 'max' of them.
5239 : * Return their indices */
5240 : static GEN
5241 1456 : row_count_non_zero(GEN a, long i, long max)
5242 : {
5243 1456 : GEN v = cgetg(max+1, t_VECSMALL);
5244 1456 : long j, l = lg(a), k = 1;
5245 5558 : for (j = 1; j < l; j++)
5246 5334 : if (!gequal0(gcoeff(a,i,j)))
5247 : {
5248 5096 : if (k > max) return NULL; /* fail */
5249 3864 : v[k++] = j;
5250 : }
5251 224 : setlg(v, k); return v;
5252 : }
5253 :
5254 : static GEN det_develop(GEN a, long max, double bound);
5255 : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
5256 : static GEN
5257 406 : coeff_det(GEN a, long i, long j, long max, double bound)
5258 : {
5259 406 : GEN c = gcoeff(a, i, j);
5260 406 : c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
5261 406 : if (odd(i+j)) c = gneg(c);
5262 406 : return c;
5263 : }
5264 : /* a square t_MAT, 'bound' a rough upper bound for the number of
5265 : * multiplications we are willing to pay while developing rows/columns before
5266 : * switching to Gaussian elimination */
5267 : static GEN
5268 658 : det_develop(GEN M, long max, double bound)
5269 : {
5270 658 : pari_sp av = avma;
5271 658 : long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
5272 658 : GEN best = NULL;
5273 :
5274 658 : if (bound < 1.) return det_bareiss(M); /* too costly now */
5275 :
5276 434 : switch(n)
5277 : {
5278 0 : case 0: return gen_1;
5279 0 : case 1: return gcopy(gcoeff(M,1,1));
5280 14 : case 2: return RgM_det2(M);
5281 : }
5282 420 : if (max > ((n+2)>>1)) max = (n+2)>>1;
5283 1876 : for (j = 1; j <= n; j++)
5284 : {
5285 1470 : pari_sp av2 = avma;
5286 1470 : GEN v = col_count_non_zero(M, j, max);
5287 : long lv;
5288 1470 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5289 182 : if (lv == 1) { set_avma(av); return gen_0; }
5290 182 : if (lv == 2) {
5291 14 : set_avma(av);
5292 14 : return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
5293 : }
5294 168 : best = v; lbest = lv; best_col = j;
5295 : }
5296 1862 : for (i = 1; i <= n; i++)
5297 : {
5298 1456 : pari_sp av2 = avma;
5299 1456 : GEN v = row_count_non_zero(M, i, max);
5300 : long lv;
5301 1456 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5302 0 : if (lv == 1) { set_avma(av); return gen_0; }
5303 0 : if (lv == 2) {
5304 0 : set_avma(av);
5305 0 : return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
5306 : }
5307 0 : best = v; lbest = lv; best_row = i;
5308 : }
5309 406 : if (best_row)
5310 : {
5311 0 : double d = lbest-1;
5312 0 : GEN s = NULL;
5313 : long k;
5314 0 : bound /= d*d*d;
5315 0 : for (k = 1; k < lbest; k++)
5316 : {
5317 0 : GEN c = coeff_det(M, best_row, best[k], max, bound);
5318 0 : s = s? gadd(s, c): c;
5319 : }
5320 0 : return gerepileupto(av, s);
5321 : }
5322 406 : if (best_col)
5323 : {
5324 168 : double d = lbest-1;
5325 168 : GEN s = NULL;
5326 : long k;
5327 168 : bound /= d*d*d;
5328 560 : for (k = 1; k < lbest; k++)
5329 : {
5330 392 : GEN c = coeff_det(M, best[k], best_col, max, bound);
5331 392 : s = s? gadd(s, c): c;
5332 : }
5333 168 : return gerepileupto(av, s);
5334 : }
5335 238 : return det_bareiss(M);
5336 : }
5337 :
5338 : /* area of parallelogram bounded by (v1,v2) */
5339 : static GEN
5340 64282 : parallelogramarea(GEN v1, GEN v2)
5341 64282 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
5342 :
5343 : /* Square of Hadamard bound for det(a), a square matrix.
5344 : * Slight improvement: instead of using the column norms, use the area of
5345 : * the parallelogram formed by pairs of consecutive vectors */
5346 : GEN
5347 19996 : RgM_Hadamard(GEN a)
5348 : {
5349 19996 : pari_sp av = avma;
5350 19996 : long n = lg(a)-1, i;
5351 : GEN B;
5352 19996 : if (n == 0) return gen_1;
5353 19996 : if (n == 1) return gsqr(gcoeff(a,1,1));
5354 19996 : a = RgM_gtofp(a, LOWDEFAULTPREC);
5355 19996 : B = gen_1;
5356 84278 : for (i = 1; i <= n/2; i++)
5357 64282 : B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
5358 19996 : if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
5359 19996 : return gerepileuptoint(av, ceil_safe(B));
5360 : }
5361 :
5362 : /* If B=NULL, assume B=A' */
5363 : static GEN
5364 20875 : ZM_det_slice(GEN A, GEN P, GEN *mod)
5365 : {
5366 20875 : pari_sp av = avma;
5367 20875 : long i, n = lg(P)-1;
5368 : GEN H, T;
5369 20875 : if (n == 1)
5370 : {
5371 0 : ulong Hp, p = uel(P,1);
5372 0 : GEN a = ZM_to_Flm(A, p);
5373 0 : Hp = Flm_det_sp(a, p);
5374 0 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
5375 : }
5376 20875 : T = ZV_producttree(P);
5377 20875 : A = ZM_nv_mod_tree(A, P, T);
5378 20875 : H = cgetg(n+1, t_VECSMALL);
5379 87556 : for(i=1; i <= n; i++)
5380 : {
5381 66681 : ulong p = P[i];
5382 66681 : GEN a = gel(A,i);
5383 66681 : H[i] = Flm_det_sp(a, p);
5384 : }
5385 20875 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
5386 20875 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
5387 : }
5388 :
5389 : GEN
5390 20875 : ZM_det_worker(GEN P, GEN A)
5391 : {
5392 20875 : GEN V = cgetg(3, t_VEC);
5393 20875 : gel(V,1) = ZM_det_slice(A, P, &gel(V,2));
5394 20875 : return V;
5395 : }
5396 :
5397 : GEN
5398 130695 : ZM_det(GEN M)
5399 : {
5400 : pari_sp av, av2;
5401 130695 : long n = lg(M)-1;
5402 : ulong p, Dp;
5403 : forprime_t S;
5404 : pari_timer ti;
5405 : GEN H, mod, h, q, worker;
5406 : #ifdef LONG_IS_64BIT
5407 112032 : const ulong PMAX = 18446744073709551557UL;
5408 : #else
5409 18663 : const ulong PMAX = 4294967291UL;
5410 : #endif
5411 :
5412 130695 : switch(n)
5413 : {
5414 7 : case 0: return gen_1;
5415 1561 : case 1: return icopy(gcoeff(M,1,1));
5416 8659 : case 2: return ZM_det2(M);
5417 100472 : case 3: return ZM_det3(M);
5418 : }
5419 19996 : if (DEBUGLEVEL>=4) timer_start(&ti);
5420 19996 : av = avma; h = RgM_Hadamard(M); /* |D| <= sqrt(h) */
5421 19996 : if (!signe(h)) { set_avma(av); return gen_0; }
5422 19996 : h = sqrti(h);
5423 19996 : if (lgefint(h) == 3 && (ulong)h[2] <= (PMAX >> 1))
5424 : { /* h < p/2 => direct result */
5425 7209 : p = PMAX;
5426 7209 : Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5427 7209 : set_avma(av);
5428 7209 : if (!Dp) return gen_0;
5429 7209 : return (Dp <= (p>>1))? utoipos(Dp): utoineg(p - Dp);
5430 : }
5431 12787 : q = gen_1; Dp = 1;
5432 12787 : init_modular_big(&S);
5433 12787 : p = 0; /* -Wall */
5434 12787 : while (cmpii(q, h) <= 0 && (p = u_forprime_next(&S)))
5435 : {
5436 12787 : av2 = avma; Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5437 12787 : set_avma(av2);
5438 12787 : if (Dp) break;
5439 0 : q = muliu(q, p);
5440 : }
5441 12787 : if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
5442 12787 : if (!Dp) { set_avma(av); return gen_0; }
5443 12787 : worker = snm_closure(is_entry("_ZM_det_worker"), mkvec(M));
5444 12787 : H = gen_crt("ZM_det", worker, &S, NULL, expi(h)+1, 0, &mod,
5445 : ZV_chinese, NULL);
5446 : /* H = det(M) modulo mod, (mod,D) = 1; |det(M) / D| <= h */
5447 12787 : H = Fp_center(H, mod, shifti(mod,-1));
5448 12787 : return gerepileuptoint(av, H);
5449 : }
5450 :
5451 : static GEN
5452 1519 : RgM_det_FpM(GEN a, GEN p)
5453 : {
5454 1519 : pari_sp av = avma;
5455 : ulong pp, d;
5456 1519 : a = RgM_Fp_init(a,p,&pp);
5457 1519 : switch(pp)
5458 : {
5459 70 : case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
5460 14 : case 2: d = F2m_det_sp(a); break;
5461 1435 : default:d = Flm_det_sp(a, pp); break;
5462 : }
5463 1449 : set_avma(av); return mkintmodu(d, pp);
5464 : }
5465 :
5466 : static GEN
5467 42 : RgM_det_FqM(GEN x, GEN pol, GEN p)
5468 : {
5469 42 : pari_sp av = avma;
5470 42 : GEN b, T = RgX_to_FpX(pol, p);
5471 42 : if (signe(T) == 0) pari_err_OP("%",x,pol);
5472 42 : b = FqM_det(RgM_to_FqM(x, T, p), T, p);
5473 42 : if (!b) return gc_NULL(av);
5474 42 : return gerepilecopy(av, mkpolmod(FpX_to_mod(b, p), FpX_to_mod(T, p)));
5475 : }
5476 :
5477 : #define code(t1,t2) ((t1 << 6) | t2)
5478 : static GEN
5479 33900 : RgM_det_fast(GEN x)
5480 : {
5481 : GEN p, pol;
5482 : long pa;
5483 33900 : long t = RgM_type(x, &p,&pol,&pa);
5484 33900 : switch(t)
5485 : {
5486 301 : case t_INT: return ZM_det(x);
5487 203 : case t_FRAC: return QM_det(x);
5488 63 : case t_FFELT: return FFM_det(x, pol);
5489 1519 : case t_INTMOD: return RgM_det_FpM(x, p);
5490 42 : case code(t_POLMOD, t_INTMOD):
5491 42 : return RgM_det_FqM(x, pol, p);
5492 31772 : default: return NULL;
5493 : }
5494 : }
5495 : #undef code
5496 :
5497 : static long
5498 252 : det_init_max(long n)
5499 : {
5500 252 : if (n > 100) return 0;
5501 252 : if (n > 50) return 1;
5502 252 : if (n > 30) return 4;
5503 252 : return 7;
5504 : }
5505 :
5506 : GEN
5507 246213 : det(GEN a)
5508 : {
5509 246213 : long n = lg(a)-1;
5510 : double B;
5511 : GEN data, b;
5512 : pivot_fun pivot;
5513 :
5514 246213 : if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
5515 246213 : if (!n) return gen_1;
5516 246171 : if (n != nbrows(a)) pari_err_DIM("det");
5517 246164 : if (n == 1) return gcopy(gcoeff(a,1,1));
5518 69053 : if (n == 2) return RgM_det2(a);
5519 33900 : b = RgM_det_fast(a);
5520 33900 : if (b) return b;
5521 31772 : pivot = get_pivot_fun(a, a, &data);
5522 31772 : if (pivot != gauss_get_pivot_NZ) return det_simple_gauss(a, data, pivot);
5523 252 : B = (double)n;
5524 252 : return det_develop(a, det_init_max(n), B*B*B);
5525 : }
5526 :
5527 : GEN
5528 203 : QM_det(GEN M)
5529 : {
5530 203 : pari_sp av = avma;
5531 203 : GEN cM, pM = Q_primitive_part(M, &cM);
5532 203 : GEN b = ZM_det(pM);
5533 203 : if (cM) b = gmul(b, gpowgs(cM, lg(M)-1));
5534 203 : return gerepileupto(av, b);
5535 : }
|
