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