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