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 1532 : gen_ker(GEN x, long deplin, void *E, const struct bb_field *ff)
124 : {
125 1532 : pari_sp av0 = avma, av, tetpil;
126 : GEN y, c, d;
127 : long i, j, k, r, t, n, m;
128 :
129 1532 : n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
130 1532 : m=nbrows(x); r=0;
131 1532 : x = RgM_shallowcopy(x);
132 1532 : c = zero_zv(m);
133 1532 : d=new_chunk(n+1);
134 1532 : av=avma;
135 5444 : for (k=1; k<=n; k++)
136 : {
137 11113 : for (j=1; j<=m; j++)
138 9394 : if (!c[j])
139 : {
140 6435 : gcoeff(x,j,k) = ff->red(E, gcoeff(x,j,k));
141 6435 : if (!ff->equal0(gcoeff(x,j,k))) break;
142 : }
143 3940 : if (j>m)
144 : {
145 1719 : if (deplin)
146 : {
147 28 : GEN c = cgetg(n+1, t_COL), g0 = ff->s(E,0), g1=ff->s(E,1);
148 63 : for (i=1; i<k; i++) gel(c,i) = ff->red(E, gcoeff(x,d[i],k));
149 28 : gel(c,k) = g1; for (i=k+1; i<=n; i++) gel(c,i) = g0;
150 28 : return gerepileupto(av0, c);
151 : }
152 1691 : r++; d[k]=0;
153 4199 : for(j=1; j<k; j++)
154 2508 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
155 : }
156 : else
157 : {
158 2221 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
159 2221 : c[j] = k; d[k] = j;
160 2221 : gcoeff(x,j,k) = ff->s(E,-1);
161 5188 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
162 10735 : for (t=1; t<=m; t++)
163 : {
164 8514 : if (t==j) continue;
165 :
166 6293 : piv = ff->red(E,gcoeff(x,t,k));
167 6293 : if (ff->equal0(piv)) continue;
168 :
169 1631 : gcoeff(x,t,k) = ff->s(E,0);
170 4224 : for (i=k+1; i<=n; i++)
171 2593 : gcoeff(x,t,i) = ff->red(E, ff->add(E, gcoeff(x,t,i),
172 2593 : ff->mul(E,piv,gcoeff(x,j,i))));
173 1631 : if (gc_needed(av,1))
174 0 : gen_gerepile_gauss_ker(x,k,t,av,E,ff->red);
175 : }
176 : }
177 : }
178 1504 : if (deplin) return gc_NULL(av0);
179 :
180 1476 : tetpil=avma; y=cgetg(r+1,t_MAT);
181 3167 : for (j=k=1; j<=r; j++,k++)
182 : {
183 1691 : GEN C = cgetg(n+1,t_COL);
184 1691 : GEN g0 = ff->s(E,0), g1 = ff->s(E,1);
185 3356 : gel(y,j) = C; while (d[k]) k++;
186 4199 : for (i=1; i<k; i++)
187 2508 : if (d[i])
188 : {
189 2101 : GEN p1=gcoeff(x,d[i],k);
190 2101 : gel(C,i) = ff->red(E,p1); gunclone(p1);
191 : }
192 : else
193 407 : gel(C,i) = g0;
194 2255 : gel(C,k) = g1; for (i=k+1; i<=n; i++) gel(C,i) = g0;
195 : }
196 1476 : return gerepile(av0,tetpil,y);
197 : }
198 :
199 : GEN
200 1548 : gen_Gauss_pivot(GEN x, long *rr, void *E, const struct bb_field *ff)
201 : {
202 : pari_sp av;
203 : GEN c, d;
204 1548 : long i, j, k, r, t, m, n = lg(x)-1;
205 :
206 1548 : if (!n) { *rr = 0; return NULL; }
207 :
208 1548 : m=nbrows(x); r=0;
209 1548 : d = cgetg(n+1, t_VECSMALL);
210 1548 : x = RgM_shallowcopy(x);
211 1548 : c = zero_zv(m);
212 1548 : av=avma;
213 5721 : for (k=1; k<=n; k++)
214 : {
215 10909 : for (j=1; j<=m; j++)
216 10621 : if (!c[j])
217 : {
218 7156 : gcoeff(x,j,k) = ff->red(E,gcoeff(x,j,k));
219 7156 : if (!ff->equal0(gcoeff(x,j,k))) break;
220 : }
221 4173 : if (j>m) { r++; d[k]=0; }
222 : else
223 : {
224 3885 : GEN piv = ff->neg(E,ff->inv(E,gcoeff(x,j,k)));
225 3885 : GEN g0 = ff->s(E,0);
226 3885 : c[j] = k; d[k] = j;
227 8092 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = ff->red(E,ff->mul(E,piv,gcoeff(x,j,i)));
228 23961 : for (t=1; t<=m; t++)
229 : {
230 20076 : if (c[t]) continue; /* already a pivot on that line */
231 :
232 12292 : piv = ff->red(E,gcoeff(x,t,k));
233 12292 : if (ff->equal0(piv)) continue;
234 4781 : gcoeff(x,t,k) = g0;
235 8694 : for (i=k+1; i<=n; i++)
236 3913 : gcoeff(x,t,i) = ff->red(E, ff->add(E,gcoeff(x,t,i), ff->mul(E,piv,gcoeff(x,j,i))));
237 4781 : if (gc_needed(av,1))
238 0 : gerepile_gauss(x,k,t,av,j,c);
239 : }
240 11977 : for (i=k; i<=n; i++) gcoeff(x,j,i) = g0; /* dummy */
241 : }
242 : }
243 1548 : *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 143373 : _gen_addmul(GEN b, long k, long i, GEN m, void *E, const struct bb_field *ff)
292 : {
293 143373 : gel(b,i) = ff->red(E,gel(b,i));
294 143373 : gel(b,k) = ff->add(E,gel(b,k), ff->mul(E,m, gel(b,i)));
295 143373 : }
296 :
297 : static GEN
298 54606 : _gen_get_col(GEN a, GEN b, long li, void *E, const struct bb_field *ff)
299 : {
300 54606 : GEN u = cgetg(li+1,t_COL);
301 54606 : pari_sp av = avma;
302 : long i, j;
303 :
304 54606 : gel(u,li) = gerepileupto(av, ff->red(E,ff->mul(E,gel(b,li), gcoeff(a,li,li))));
305 284283 : for (i=li-1; i>0; i--)
306 : {
307 229677 : pari_sp av = avma;
308 229677 : GEN m = gel(b,i);
309 919952 : 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 229677 : m = ff->red(E, m);
311 229677 : gel(u,i) = gerepileupto(av, ff->red(E,ff->mul(E,m, gcoeff(a,i,i))));
312 : }
313 54606 : return u;
314 : }
315 :
316 : GEN
317 12119 : gen_Gauss(GEN a, GEN b, void *E, const struct bb_field *ff)
318 : {
319 : long i, j, k, li, bco, aco;
320 12119 : GEN u, g0 = ff->s(E,0);
321 12119 : pari_sp av = avma;
322 12119 : a = RgM_shallowcopy(a);
323 12119 : b = RgM_shallowcopy(b);
324 12119 : aco = lg(a)-1; bco = lg(b)-1; li = nbrows(a);
325 53465 : for (i=1; i<=aco; i++)
326 : {
327 : GEN invpiv;
328 64413 : for (k = i; k <= li; k++)
329 : {
330 64357 : GEN piv = ff->red(E,gcoeff(a,k,i));
331 64357 : if (!ff->equal0(piv)) { gcoeff(a,k,i) = ff->inv(E,piv); break; }
332 10948 : gcoeff(a,k,i) = g0;
333 : }
334 : /* found a pivot on line k */
335 53465 : if (k > li) return NULL;
336 53409 : if (k != i)
337 : { /* swap lines so that k = i */
338 44149 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
339 59234 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
340 : }
341 53409 : if (i == aco) break;
342 :
343 41346 : invpiv = gcoeff(a,i,i); /* 1/piv mod p */
344 156937 : for (k=i+1; k<=li; k++)
345 : {
346 115591 : GEN m = ff->red(E,gcoeff(a,k,i)); gcoeff(a,k,i) = g0;
347 115591 : if (ff->equal0(m)) continue;
348 :
349 17591 : m = ff->red(E,ff->neg(E,ff->mul(E,m, invpiv)));
350 70596 : for (j=i+1; j<=aco; j++) _gen_addmul(gel(a,j),k,i,m,E,ff);
351 107959 : for (j=1 ; j<=bco; j++) _gen_addmul(gel(b,j),k,i,m,E,ff);
352 : }
353 41346 : 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 12063 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
361 12063 : u = cgetg(bco+1,t_MAT);
362 66669 : for (j=1; j<=bco; j++) gel(u,j) = _gen_get_col(a, gel(b,j), aco, E, ff);
363 12063 : return u;
364 : }
365 :
366 : /* compatible t_MAT * t_COL, lgA = lg(A) = lg(B) > 1, l = lgcols(A) */
367 : static GEN
368 575767 : gen_matcolmul_i(GEN A, GEN B, ulong lgA, ulong l,
369 : void *E, const struct bb_field *ff)
370 : {
371 575767 : GEN C = cgetg(l, t_COL);
372 : ulong i;
373 3852221 : for (i = 1; i < l; i++) {
374 3276454 : pari_sp av = avma;
375 3276454 : GEN e = ff->mul(E, gcoeff(A, i, 1), gel(B, 1));
376 : ulong k;
377 14080464 : for(k = 2; k < lgA; k++)
378 10804010 : e = ff->add(E, e, ff->mul(E, gcoeff(A, i, k), gel(B, k)));
379 3276454 : gel(C, i) = gerepileupto(av, ff->red(E, e));
380 : }
381 575767 : return C;
382 : }
383 :
384 : GEN
385 179676 : gen_matcolmul(GEN A, GEN B, void *E, const struct bb_field *ff)
386 : {
387 179676 : ulong lgA = lg(A);
388 179676 : if (lgA != (ulong)lg(B))
389 0 : pari_err_OP("operation 'gen_matcolmul'", A, B);
390 179676 : if (lgA == 1)
391 0 : return cgetg(1, t_COL);
392 179676 : return gen_matcolmul_i(A, B, lgA, lgcols(A), E, ff);
393 : }
394 :
395 : static GEN
396 76264 : 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 76264 : GEN C = cgetg(lb, t_MAT);
401 472355 : for(j = 1; j < lb; j++)
402 396091 : gel(C, j) = gen_matcolmul_i(A, gel(B, j), la, l, E, ff);
403 76264 : 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 76264 : gen_matmul_i(GEN A, GEN B, long l, long la, long lb,
587 : void *E, const struct bb_field *ff)
588 : {
589 76264 : if (l <= gen_matmul_sw_bound
590 14 : || la <= gen_matmul_sw_bound
591 0 : || lb <= gen_matmul_sw_bound)
592 76264 : 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 76264 : gen_matmul(GEN A, GEN B, void *E, const struct bb_field *ff)
599 : {
600 76264 : ulong lgA, lgB = lg(B);
601 76264 : if (lgB == 1)
602 0 : return cgetg(1, t_MAT);
603 76264 : lgA = lg(A);
604 76264 : if (lgA != (ulong)lgcols(B))
605 0 : pari_err_OP("operation 'gen_matmul'", A, B);
606 76264 : if (lgA == 1)
607 0 : return zeromat(0, lgB - 1);
608 76264 : return gen_matmul_i(A, B, lgcols(A), lgA, lgB, E, ff);
609 : }
610 :
611 : static GEN
612 18217 : gen_colneg(GEN A, void *E, const struct bb_field *ff)
613 : {
614 : long i, l;
615 18217 : GEN B = cgetg_copy(A, &l);
616 72147 : for (i = 1; i < l; i++)
617 53930 : gel(B, i) = ff->neg(E, gel(A, i));
618 18217 : return B;
619 : }
620 :
621 : static GEN
622 3945 : gen_matneg(GEN A, void *E, const struct bb_field *ff)
623 : {
624 : long i, l;
625 3945 : GEN B = cgetg_copy(A, &l);
626 22106 : for (i = 1; i < l; i++)
627 18161 : gel(B, i) = gen_colneg(gel(A, i), E, ff);
628 3945 : return B;
629 : }
630 :
631 : static GEN
632 253851 : gen_colscalmul(GEN A, GEN b, void *E, const struct bb_field *ff)
633 : {
634 : long i, l;
635 253851 : GEN B = cgetg_copy(A, &l);
636 612771 : for (i = 1; i < l; i++)
637 358920 : gel(B, i) = ff->red(E, ff->mul(E, gel(A, i), b));
638 253851 : return B;
639 : }
640 :
641 : static GEN
642 48849 : gen_matscalmul(GEN A, GEN b, void *E, const struct bb_field *ff)
643 : {
644 : long i, l;
645 48849 : GEN B = cgetg_copy(A, &l);
646 302700 : for (i = 1; i < l; i++)
647 253851 : gel(B, i) = gen_colscalmul(gel(A, i), b, E, ff);
648 48849 : return B;
649 : }
650 :
651 : static GEN
652 481623 : gen_colsub(GEN A, GEN C, void *E, const struct bb_field *ff)
653 : {
654 : long i, l;
655 481623 : GEN B = cgetg_copy(A, &l);
656 1737253 : for (i = 1; i < l; i++)
657 1255630 : gel(B, i) = ff->add(E, gel(A, i), ff->neg(E, gel(C, i)));
658 481623 : return B;
659 : }
660 :
661 : static GEN
662 70746 : gen_matsub(GEN A, GEN C, void *E, const struct bb_field *ff)
663 : {
664 : long i, l;
665 70746 : GEN B = cgetg_copy(A, &l);
666 552369 : for (i = 1; i < l; i++)
667 481623 : gel(B, i) = gen_colsub(gel(A, i), gel(C, i), E, ff);
668 70746 : return B;
669 : }
670 :
671 : static GEN
672 42745 : gen_zerocol(long n, void* data, const struct bb_field *R)
673 : {
674 42745 : GEN C = cgetg(n+1,t_COL), zero = R->s(data, 0);
675 : long i;
676 253906 : for (i=1; i<=n; i++) gel(C,i) = zero;
677 42745 : return C;
678 : }
679 :
680 : static GEN
681 13599 : gen_zeromat(long m, long n, void* data, const struct bb_field *R)
682 : {
683 13599 : GEN M = cgetg(n+1,t_MAT);
684 : long i;
685 56344 : for (i=1; i<=n; i++) gel(M,i) = gen_zerocol(m, data, R);
686 13599 : return M;
687 : }
688 :
689 : static GEN
690 140 : gen_colei(long n, long i, void *E, const struct bb_field *S)
691 : {
692 140 : GEN y = cgetg(n+1,t_COL), _0, _1;
693 : long j;
694 140 : if (n < 0) pari_err_DOMAIN("gen_colei", "dimension","<",gen_0,stoi(n));
695 140 : _0 = S->s(E,0);
696 140 : _1 = S->s(E,1);
697 2268 : for (j=1; j<=n; j++)
698 2128 : gel(y, j) = i==j ? _1: _0;
699 140 : 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 317153 : image_from_pivot(GEN x, GEN d, long r)
785 : {
786 : GEN y;
787 : long j, k;
788 :
789 317153 : if (!d) return gcopy(x);
790 : /* d left on stack for efficiency */
791 314493 : r = lg(x)-1 - r; /* = dim Im(x) */
792 314493 : y = cgetg(r+1,t_MAT);
793 1908625 : for (j=k=1; j<=r; k++)
794 1594134 : if (d[k]) gel(y,j++) = gcopy(gel(x,k));
795 314491 : return y;
796 : }
797 :
798 : /* r = dim Ker x, n = nbrows(x) */
799 : static GEN
800 244030 : get_suppl(GEN x, GEN d, long n, long r, GEN(*ei)(long,long))
801 : {
802 : pari_sp av;
803 : GEN y, c;
804 244030 : long j, k, rx = lg(x)-1; /* != 0 due to init_suppl() */
805 :
806 244030 : if (rx == n && r == 0) return gcopy(x);
807 177701 : y = cgetg(n+1, t_MAT);
808 177704 : 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 755188 : for (k = j = 1; j<=rx; j++)
812 577484 : if (d[j]) { c[ d[j] ] = 1; gel(y,k++) = gel(x,j); }
813 1063525 : for (j=1; j<=n; j++)
814 885821 : if (!c[j]) gel(y,k++) = (GEN)j; /* HACK */
815 177704 : set_avma(av);
816 :
817 177703 : rx -= r;
818 755117 : for (j=1; j<=rx; j++) gel(y,j) = gcopy(gel(y,j));
819 486111 : for ( ; j<=n; j++) gel(y,j) = ei(n, y[j]);
820 177704 : return y;
821 : }
822 :
823 : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
824 : static GEN
825 191783 : indexrank0(long n, long r, GEN d)
826 : {
827 191783 : GEN p1, p2, res = cgetg(3,t_VEC);
828 : long i, j;
829 :
830 191784 : r = n - r; /* now r = dim Im(x) */
831 191784 : p1 = cgetg(r+1,t_VECSMALL); gel(res,1) = p1;
832 191784 : p2 = cgetg(r+1,t_VECSMALL); gel(res,2) = p2;
833 191784 : if (d)
834 : {
835 1087968 : for (i=0,j=1; j<=n; j++)
836 897500 : if (d[j]) { i++; p1[i] = d[j]; p2[i] = j; }
837 190468 : vecsmall_sort(p1);
838 : }
839 191784 : 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 16919 : indexcompl(GEN v, long n)
867 : {
868 16919 : long i, j, k, m = lg(v) - 1;
869 16919 : GEN w = cgetg(n - m + 1, t_VECSMALL);
870 154243 : for (i = j = k = 1; i <= n; i++)
871 137324 : if (j <= m && v[j] == i) j++; else w[k++] = i;
872 16919 : return w;
873 : }
874 :
875 : static GEN
876 3727 : gen_solve_upper_1(GEN U, GEN B, void *E, const struct bb_field *ff)
877 3727 : { return gen_matscalmul(B, ff->inv(E, gcoeff(U, 1, 1)), E, ff); }
878 :
879 : static GEN
880 1948 : gen_rsolve_upper_2(GEN U, GEN B, void *E, const struct bb_field *ff)
881 : {
882 1948 : GEN a = gcoeff(U, 1, 1), b = gcoeff(U, 1, 2), d = gcoeff(U, 2, 2);
883 1948 : GEN D = ff->red(E, ff->mul(E, a, d)), Dinv = ff->inv(E, D);
884 1948 : GEN ainv = ff->red(E, ff->mul(E, d, Dinv));
885 1948 : GEN dinv = ff->red(E, ff->mul(E, a, Dinv));
886 1948 : GEN B1 = rowslice(B, 1, 1);
887 1948 : GEN B2 = rowslice(B, 2, 2);
888 1948 : GEN X2 = gen_matscalmul(B2, dinv, E, ff);
889 1948 : GEN X1 = gen_matscalmul(gen_matsub(B1, gen_matscalmul(X2, b, E, ff), E, ff),
890 : ainv, E, ff);
891 1948 : return vconcat(X1, X2);
892 : }
893 :
894 : /* solve U*X = B, U upper triangular and invertible */
895 : static GEN
896 5132 : gen_rsolve_upper(GEN U, GEN B, void *E, const struct bb_field *ff,
897 : GEN (*mul)(void *E, GEN a, GEN))
898 : {
899 5132 : long n = lg(U) - 1, n1;
900 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
901 5132 : pari_sp av = avma;
902 :
903 5132 : if (n == 0) return B;
904 5132 : if (n == 1) return gen_solve_upper_1(U, B, E, ff);
905 4280 : if (n == 2) return gen_rsolve_upper_2(U, B, E, ff);
906 2332 : n1 = (n + 1)/2;
907 2332 : U2 = vecslice(U, n1 + 1, n);
908 2332 : U11 = matslice(U, 1,n1, 1,n1);
909 2332 : U12 = rowslice(U2, 1, n1);
910 2332 : U22 = rowslice(U2, n1 + 1, n);
911 2332 : B1 = rowslice(B, 1, n1);
912 2332 : B2 = rowslice(B, n1 + 1, n);
913 2332 : X2 = gen_rsolve_upper(U22, B2, E, ff, mul);
914 2332 : B1 = gen_matsub(B1, mul(E, U12, X2), E, ff);
915 2332 : if (gc_needed(av, 1)) gerepileall(av, 3, &B1, &U11, &X2);
916 2332 : X1 = gen_rsolve_upper(U11, B1, E, ff, mul);
917 2332 : X = vconcat(X1, X2);
918 2332 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
919 2332 : return X;
920 : }
921 :
922 : static GEN
923 5615 : gen_lsolve_upper_2(GEN U, GEN B, void *E, const struct bb_field *ff)
924 : {
925 5615 : GEN a = gcoeff(U, 1, 1), b = gcoeff(U, 1, 2), d = gcoeff(U, 2, 2);
926 5615 : GEN D = ff->red(E, ff->mul(E, a, d)), Dinv = ff->inv(E, D);
927 5615 : GEN ainv = ff->red(E, ff->mul(E, d, Dinv)), dinv = ff->red(E, ff->mul(E, a, Dinv));
928 5615 : GEN B1 = vecslice(B, 1, 1);
929 5615 : GEN B2 = vecslice(B, 2, 2);
930 5615 : GEN X1 = gen_matscalmul(B1, ainv, E, ff);
931 5615 : GEN X2 = gen_matscalmul(gen_matsub(B2, gen_matscalmul(X1, b, E, ff), E, ff), dinv, E, ff);
932 5615 : return shallowconcat(X1, X2);
933 : }
934 :
935 : /* solve X*U = B, U upper triangular and invertible */
936 : static GEN
937 12988 : gen_lsolve_upper(GEN U, GEN B, void *E, const struct bb_field *ff,
938 : GEN (*mul)(void *E, GEN a, GEN))
939 : {
940 12988 : long n = lg(U) - 1, n1;
941 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
942 12988 : pari_sp av = avma;
943 :
944 12988 : if (n == 0) return B;
945 12988 : if (n == 1) return gen_solve_upper_1(U, B, E, ff);
946 10113 : if (n == 2) return gen_lsolve_upper_2(U, B, E, ff);
947 4498 : n1 = (n + 1)/2;
948 4498 : U2 = vecslice(U, n1 + 1, n);
949 4498 : U11 = matslice(U, 1,n1, 1,n1);
950 4498 : U12 = rowslice(U2, 1, n1);
951 4498 : U22 = rowslice(U2, n1 + 1, n);
952 4498 : B1 = vecslice(B, 1, n1);
953 4498 : B2 = vecslice(B, n1 + 1, n);
954 4498 : X1 = gen_lsolve_upper(U11, B1, E, ff, mul);
955 4498 : B2 = gen_matsub(B2, mul(E, X1, U12), E, ff);
956 4498 : if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
957 4498 : X2 = gen_lsolve_upper(U22, B2, E, ff, mul);
958 4498 : X = shallowconcat(X1, X2);
959 4498 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
960 4498 : return X;
961 : }
962 :
963 : static GEN
964 15175 : gen_rsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
965 : {
966 15175 : GEN X1 = rowslice(A, 1, 1);
967 15175 : GEN X2 = gen_matsub(rowslice(A, 2, 2), gen_matscalmul(X1, gcoeff(L, 2, 1), E, ff), E, ff);
968 15175 : 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 35527 : 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 35527 : long m = lg(L) - 1, m1, n;
978 : GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
979 35527 : pari_sp av = avma;
980 :
981 35527 : if (m == 0) return zeromat(0, lg(A) - 1);
982 35527 : if (m == 1) return rowslice(A, 1, 1);
983 27733 : if (m == 2) return gen_rsolve_lower_unit_2(L, A, E, ff);
984 12558 : m1 = (m + 1)/2;
985 12558 : n = nbrows(L);
986 12558 : L1 = vecslice(L, 1, m1);
987 12558 : L11 = rowslice(L1, 1, m1);
988 12558 : L21 = rowslice(L1, m1 + 1, n);
989 12558 : A1 = rowslice(A, 1, m1);
990 12558 : X1 = gen_rsolve_lower_unit(L11, A1, E, ff, mul);
991 12558 : A2 = rowslice(A, m1 + 1, n);
992 12558 : A2 = gen_matsub(A2, mul(E, L21, X1), E, ff);
993 12558 : if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
994 12558 : L22 = matslice(L, m1+1,n, m1+1,m);
995 12558 : X2 = gen_rsolve_lower_unit(L22, A2, E, ff, mul);
996 12558 : X = vconcat(X1, X2);
997 12558 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
998 12558 : return X;
999 : }
1000 :
1001 : static GEN
1002 7258 : gen_lsolve_lower_unit_2(GEN L, GEN A, void *E, const struct bb_field *ff)
1003 : {
1004 7258 : GEN X2 = vecslice(A, 2, 2);
1005 7258 : GEN X1 = gen_matsub(vecslice(A, 1, 1),
1006 7258 : gen_matscalmul(X2, gcoeff(L, 2, 1), E, ff), E, ff);
1007 7258 : 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 18862 : 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 18862 : long m = lg(L) - 1, m1;
1017 : GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
1018 18862 : pari_sp av = avma;
1019 :
1020 18862 : if (m <= 1) return A;
1021 14727 : if (m == 2) return gen_lsolve_lower_unit_2(L, A, E, ff);
1022 7469 : m1 = (m + 1)/2;
1023 7469 : L2 = vecslice(L, m1 + 1, m);
1024 7469 : L22 = rowslice(L2, m1 + 1, m);
1025 7469 : A2 = vecslice(A, m1 + 1, m);
1026 7469 : X2 = gen_lsolve_lower_unit(L22, A2, E, ff, mul);
1027 7469 : if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
1028 7469 : L1 = vecslice(L, 1, m1);
1029 7469 : L21 = rowslice(L1, m1 + 1, m);
1030 7469 : A1 = vecslice(A, 1, m1);
1031 7469 : A1 = gen_matsub(A1, mul(E, X2, L21), E, ff);
1032 7469 : L11 = rowslice(L1, 1, m1);
1033 7469 : if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
1034 7469 : X1 = gen_lsolve_lower_unit(L11, A1, E, ff, mul);
1035 7469 : X = shallowconcat(X1, X2);
1036 7469 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
1037 7469 : return X;
1038 : }
1039 :
1040 : /* destroy A */
1041 : static long
1042 20886 : gen_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, void *E, const struct bb_field *ff)
1043 : {
1044 20886 : long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
1045 : pari_sp av;
1046 : GEN u, v;
1047 :
1048 20886 : if (P) *P = identity_perm(n);
1049 20886 : *R = cgetg(m + 1, t_VECSMALL);
1050 20886 : av = avma;
1051 54390 : for (j = 1, pr = 0; j <= n; j++)
1052 : {
1053 123897 : for (pr++, pc = 0; pr <= m; pr++)
1054 : {
1055 568981 : for (k = j; k <= n; k++)
1056 : {
1057 461422 : v = ff->red(E, gcoeff(A, pr, k));
1058 461422 : gcoeff(A, pr, k) = v;
1059 461422 : if (!pc && !ff->equal0(v)) pc = k;
1060 : }
1061 107559 : if (pc) break;
1062 : }
1063 49842 : if (!pc) break;
1064 33504 : (*R)[j] = pr;
1065 33504 : if (pc != j)
1066 : {
1067 4400 : swap(gel(A, j), gel(A, pc));
1068 4400 : if (P) lswap((*P)[j], (*P)[pc]);
1069 : }
1070 33504 : u = ff->inv(E, gcoeff(A, pr, j));
1071 162603 : for (i = pr + 1; i <= m; i++)
1072 : {
1073 129099 : v = ff->red(E, ff->mul(E, gcoeff(A, i, j), u));
1074 129099 : gcoeff(A, i, j) = v;
1075 129099 : v = ff->neg(E, v);
1076 407065 : for (k = j + 1; k <= n; k++)
1077 277966 : gcoeff(A, i, k) = ff->add(E, gcoeff(A, i, k),
1078 277966 : ff->red(E, ff->mul(E, gcoeff(A, pr, k), v)));
1079 : }
1080 33504 : if (gc_needed(av, 2)) A = gerepilecopy(av, A);
1081 : }
1082 20886 : setlg(*R, j);
1083 20886 : *C = vecslice(A, 1, j - 1);
1084 20886 : if (U) *U = rowpermute(A, *R);
1085 20886 : return j - 1;
1086 : }
1087 :
1088 : static const long gen_CUP_LIMIT = 5;
1089 :
1090 : static long
1091 10346 : 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 10346 : 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 10346 : pari_sp av = avma;
1098 :
1099 10346 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1100 : /* destroy A; not called at the outermost recursion level */
1101 5940 : return gen_CUP_basecase(A, R, C, U, P, E, ff);
1102 4406 : m1 = (minss(m, n) + 1)/2;
1103 4406 : A1 = rowslice(A, 1, m1);
1104 4406 : A2 = rowslice(A, m1 + 1, m);
1105 4406 : r1 = gen_CUP(A1, &R1, &C1, &U1, &P1, E, ff, mul);
1106 4406 : if (r1 == 0)
1107 : {
1108 414 : r2 = gen_CUP(A2, &R2, &C2, &U2, &P2, E, ff, mul);
1109 414 : *R = cgetg(r2 + 1, t_VECSMALL);
1110 717 : for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
1111 414 : *C = vconcat(gen_zeromat(m1, r2, E, ff), C2);
1112 414 : *U = U2;
1113 414 : *P = P2;
1114 414 : r = r2;
1115 : }
1116 : else
1117 : {
1118 3992 : U11 = vecslice(U1, 1, r1);
1119 3992 : U12 = vecslice(U1, r1 + 1, n);
1120 3992 : T21 = vecslicepermute(A2, P1, 1, r1);
1121 3992 : T22 = vecslicepermute(A2, P1, r1 + 1, n);
1122 3992 : C21 = gen_lsolve_upper(U11, T21, E, ff, mul);
1123 3992 : if (gc_needed(av, 1))
1124 0 : gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
1125 3992 : B2 = gen_matsub(T22, mul(E, C21, U12), E, ff);
1126 3992 : r2 = gen_CUP(B2, &R2, &C2, &U2, &P2, E, ff, mul);
1127 3992 : r = r1 + r2;
1128 3992 : *R = cgetg(r + 1, t_VECSMALL);
1129 18097 : for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
1130 19543 : for ( ; i <= r; i++) (*R)[i] = R2[i - r1] + m1;
1131 3992 : *C = shallowconcat(vconcat(C1, C21),
1132 : vconcat(gen_zeromat(m1, r2, E, ff), C2));
1133 3992 : *U = shallowconcat(vconcat(U11, gen_zeromat(r2, r1, E, ff)),
1134 : vconcat(vecpermute(U12, P2), U2));
1135 :
1136 3992 : *P = cgetg(n + 1, t_VECSMALL);
1137 18097 : for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
1138 43341 : for ( ; i <= n; i++) (*P)[i] = P1[P2[i - r1] + r1];
1139 : }
1140 4406 : if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
1141 4406 : return r;
1142 : }
1143 :
1144 : /* column echelon form */
1145 : static long
1146 26063 : gen_echelon(GEN A, GEN *R, GEN *C, void *E, const struct bb_field *ff,
1147 : GEN (*mul)(void*, GEN, GEN))
1148 : {
1149 26063 : 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 26063 : pari_sp av = avma;
1153 :
1154 26063 : if (m < gen_CUP_LIMIT || n < gen_CUP_LIMIT)
1155 14946 : return gen_CUP_basecase(shallowcopy(A), R, C, NULL, NULL, E, ff);
1156 :
1157 11117 : n1 = (n + 1)/2;
1158 11117 : A1 = vecslice(A, 1, n1);
1159 11117 : A2 = vecslice(A, n1 + 1, n);
1160 11117 : r1 = gen_echelon(A1, &R1, &C1, E, ff, mul);
1161 11117 : if (!r1) return gen_echelon(A2, R, C, E, ff, mul);
1162 10052 : if (r1 == m) { *R = R1; *C = C1; return r1; }
1163 9901 : R1c = indexcompl(R1, m);
1164 9901 : C11 = rowpermute(C1, R1);
1165 9901 : C21 = rowpermute(C1, R1c);
1166 9901 : A12 = rowpermute(A2, R1);
1167 9901 : A22 = rowpermute(A2, R1c);
1168 9901 : M12 = gen_rsolve_lower_unit(C11, A12, E, ff, mul);
1169 9901 : B2 = gen_matsub(A22, mul(E, C21, M12), E, ff);
1170 9901 : r2 = gen_echelon(B2, &R2, &C2, E, ff, mul);
1171 9901 : if (!r2) { *R = R1; *C = C1; r = r1; }
1172 : else
1173 : {
1174 5180 : R2 = perm_mul(R1c, R2);
1175 5180 : C2 = rowpermute(vconcat(gen_zeromat(r1, r2, E, ff), C2),
1176 : perm_inv(vecsmall_concat(R1, R1c)));
1177 5180 : r = r1 + r2;
1178 5180 : *R = cgetg(r + 1, t_VECSMALL);
1179 5180 : *C = cgetg(r + 1, t_MAT);
1180 35586 : for (j = j1 = j2 = 1; j <= r; j++)
1181 30406 : if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
1182 : {
1183 17872 : gel(*C, j) = gel(C1, j1);
1184 17872 : (*R)[j] = R1[j1++];
1185 : }
1186 : else
1187 : {
1188 12534 : gel(*C, j) = gel(C2, j2);
1189 12534 : (*R)[j] = R2[j2++];
1190 : }
1191 : }
1192 9901 : if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
1193 9901 : return r;
1194 : }
1195 :
1196 : static GEN
1197 828 : 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 828 : long i, n = lg(x) - 1, r;
1202 828 : GEN R, C, U, P, d = zero_zv(n);
1203 828 : av = avma;
1204 828 : r = gen_CUP(x, &R, &C, &U, &P, E, ff, mul);
1205 6383 : for(i = 1; i <= r; i++)
1206 5555 : d[P[i]] = R[i];
1207 828 : set_avma(av);
1208 828 : *rr = n - r;
1209 828 : 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 3868 : gen_ker_echelon(GEN x, void *E, const struct bb_field *ff,
1268 : GEN (*mul)(void*, GEN, GEN))
1269 : {
1270 3868 : pari_sp av = avma;
1271 : GEN R, Rc, C, C1, C2, S, K;
1272 3868 : long n = lg(x) - 1, r;
1273 3868 : r = gen_echelon(shallowtrans(x), &R, &C, E, ff, mul);
1274 3868 : Rc = indexcompl(R, n);
1275 3868 : C1 = rowpermute(C, R);
1276 3868 : C2 = rowpermute(C, Rc);
1277 3868 : S = gen_lsolve_lower_unit(C1, C2, E, ff, mul);
1278 3868 : K = vecpermute(shallowconcat(gen_matneg(S, E, ff), gen_matid(n - r, E, ff)),
1279 : perm_inv(vecsmall_concat(R, Rc)));
1280 3868 : K = shallowtrans(K);
1281 3868 : return gerepilecopy(av, K);
1282 : }
1283 :
1284 : static GEN
1285 84 : gen_deplin_echelon(GEN x, void *E, const struct bb_field *ff,
1286 : GEN (*mul)(void*, GEN, GEN))
1287 : {
1288 84 : pari_sp av = avma;
1289 : GEN R, Rc, C, C1, C2, s, v;
1290 84 : long i, n = lg(x) - 1, r;
1291 84 : r = gen_echelon(shallowtrans(x), &R, &C, E, ff, mul);
1292 84 : if (r == n) return gc_NULL(av);
1293 56 : Rc = indexcompl(R, n);
1294 56 : i = Rc[1];
1295 56 : C1 = rowpermute(C, R);
1296 56 : C2 = rowslice(C, i, i);
1297 56 : s = row(gen_lsolve_lower_unit(C1, C2, E, ff, mul), 1);
1298 56 : settyp(s, t_COL);
1299 56 : v = vecpermute(shallowconcat(gen_colneg(s, E, ff), gen_colei(n - r, 1, E, ff)),
1300 : perm_inv(vecsmall_concat(R, Rc)));
1301 56 : return gerepilecopy(av, v);
1302 : }
1303 :
1304 : static GEN
1305 503 : 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 503 : long n = lg(a) - 1, r;
1310 503 : if (nbrows(a) < n || (r = gen_CUP(a, &R, &C, &U, &P, E, ff, mul)) < n)
1311 56 : return NULL;
1312 447 : Y = gen_rsolve_lower_unit(rowpermute(C, R), rowpermute(b, R), E, ff, mul);
1313 447 : return rowpermute(gen_rsolve_upper(U, Y, E, ff, mul), perm_inv(P));
1314 : }
1315 :
1316 : static GEN
1317 3830 : gen_gauss(GEN a, GEN b, void *E, const struct bb_field *ff,
1318 : GEN (*mul)(void*, GEN, GEN))
1319 : {
1320 3830 : if (lg(a) - 1 >= gen_CUP_LIMIT)
1321 503 : return gen_gauss_CUP(a, b, E, ff, mul);
1322 3327 : return gen_Gauss(a, b, E, ff);
1323 : }
1324 :
1325 : static GEN
1326 5407 : gen_ker_i(GEN x, long deplin, void *E, const struct bb_field *ff,
1327 : GEN (*mul)(void*, GEN, GEN)) {
1328 5407 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1329 3952 : return deplin? gen_deplin_echelon(x, E, ff, mul): gen_ker_echelon(x, E, ff, mul);
1330 1455 : 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 70 : gen_matcolinvimage_i(GEN A, GEN y, void *E, const struct bb_field *ff,
1348 : GEN (*mul)(void*, GEN, GEN))
1349 : {
1350 70 : pari_sp av = avma;
1351 70 : long i, l = lg(A);
1352 : GEN M, x, t;
1353 :
1354 70 : M = gen_ker_i(shallowconcat(A, y), 0, E, ff, mul);
1355 70 : i = lg(M) - 1;
1356 70 : if (!i) return gc_NULL(av);
1357 :
1358 70 : x = gel(M, i);
1359 70 : t = gel(x, l);
1360 70 : if (ff->equal0(t)) return gc_NULL(av);
1361 :
1362 49 : t = ff->neg(E, ff->inv(E, t));
1363 49 : setlg(x, l);
1364 175 : for (i = 1; i < l; i++)
1365 126 : gel(x, i) = ff->red(E, ff->mul(E, t, gel(x, i)));
1366 49 : 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 2369 : gen_pivots(GEN x, long *rr, void *E, const struct bb_field *ff,
1381 : GEN (*mul)(void*, GEN, GEN))
1382 : {
1383 2369 : if (lg(x) - 1 >= gen_CUP_LIMIT && nbrows(x) >= gen_CUP_LIMIT)
1384 828 : return gen_pivots_CUP(x, rr, E, ff, mul);
1385 1541 : 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 49 : FqM_FqC_mul(GEN A, GEN B, GEN T, GEN p) {
1441 : void *E;
1442 49 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1443 49 : 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 150196 : FlxqM_mul(GEN A, GEN B, GEN T, ulong p) {
1455 : void *E;
1456 : const struct bb_field *ff;
1457 150196 : long n = lg(A) - 1;
1458 :
1459 150196 : if (n == 0)
1460 0 : return cgetg(1, t_MAT);
1461 150196 : if (n > 1)
1462 87512 : return FlxqM_mul_Kronecker(A, B, T, p);
1463 62684 : ff = get_Flxq_field(&E, T, p);
1464 62684 : return gen_matmul(A, B, E, ff);
1465 : }
1466 :
1467 : GEN
1468 66423 : FqM_mul(GEN A, GEN B, GEN T, GEN p) {
1469 : void *E;
1470 66423 : long n = lg(A) - 1;
1471 : const struct bb_field *ff;
1472 66423 : if (n == 0)
1473 0 : return cgetg(1, t_MAT);
1474 66423 : if (n > 1)
1475 62804 : 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 16422 : _FlxqM_mul(void *E, GEN A, GEN B)
1498 : {
1499 16422 : struct _Flxq *D = (struct _Flxq*)E;
1500 16422 : return FlxqM_mul(A, B, D->T, D->p);
1501 : }
1502 :
1503 : static GEN
1504 18042 : _FpM_mul(void *E, GEN A, GEN B)
1505 18042 : { 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 2218705 : FpM_init(GEN a, GEN p, ulong *pp)
1521 : {
1522 2218705 : if (lgefint(p) == 3)
1523 : {
1524 2211872 : *pp = uel(p,2);
1525 2211872 : return (*pp==2)? ZM_to_F2m(a): ZM_to_Flm(a, *pp);
1526 : }
1527 6833 : *pp = 0; return a;
1528 : }
1529 : GEN
1530 2359 : RgM_Fp_init(GEN a, GEN p, ulong *pp)
1531 : {
1532 2359 : if (lgefint(p) == 3)
1533 : {
1534 2009 : *pp = uel(p,2);
1535 2009 : return (*pp==2)? RgM_to_F2m(a): RgM_to_Flm(a, *pp);
1536 : }
1537 350 : *pp = 0; return RgM_to_FpM(a,p);
1538 : }
1539 :
1540 : static GEN
1541 315 : FpM_det_gen(GEN a, GEN p)
1542 : {
1543 : void *E;
1544 315 : const struct bb_field *S = get_Fp_field(&E,p);
1545 315 : return gen_det_i(a, E, S, _FpM_mul);
1546 : }
1547 : GEN
1548 3948 : FpM_det(GEN a, GEN p)
1549 : {
1550 3948 : pari_sp av = avma;
1551 : ulong pp, d;
1552 3948 : a = FpM_init(a, p, &pp);
1553 3948 : switch(pp)
1554 : {
1555 315 : case 0: return FpM_det_gen(a, p);
1556 1617 : case 2: d = F2m_det_sp(a); break;
1557 2016 : default:d = Flm_det_sp(a,pp); break;
1558 : }
1559 3633 : set_avma(av); return utoi(d);
1560 : }
1561 :
1562 : GEN
1563 7 : F2xqM_det(GEN a, GEN T)
1564 : {
1565 : void *E;
1566 7 : const struct bb_field *S = get_F2xq_field(&E, T);
1567 7 : return gen_det_i(a, E, S, _F2xqM_mul);
1568 : }
1569 :
1570 : GEN
1571 28 : FlxqM_det(GEN a, GEN T, ulong p) {
1572 : void *E;
1573 28 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1574 28 : return gen_det_i(a, E, S, _FlxqM_mul);
1575 : }
1576 :
1577 : GEN
1578 70 : FqM_det(GEN x, GEN T, GEN p)
1579 : {
1580 : void *E;
1581 70 : const struct bb_field *S = get_Fq_field(&E,T,p);
1582 70 : return gen_det_i(x, E, S, _FqM_mul);
1583 : }
1584 :
1585 : static GEN
1586 815 : FpM_gauss_pivot_gen(GEN x, GEN p, long *rr)
1587 : {
1588 : void *E;
1589 815 : const struct bb_field *S = get_Fp_field(&E,p);
1590 815 : return gen_pivots(x, rr, E, S, _FpM_mul);
1591 : }
1592 :
1593 : static GEN
1594 562461 : FpM_gauss_pivot(GEN x, GEN p, long *rr)
1595 : {
1596 : ulong pp;
1597 562461 : if (lg(x)==1) { *rr = 0; return NULL; }
1598 560298 : x = FpM_init(x, p, &pp);
1599 560296 : switch(pp)
1600 : {
1601 815 : case 0: return FpM_gauss_pivot_gen(x, p, rr);
1602 334367 : case 2: return F2m_gauss_pivot(x, rr);
1603 225114 : default:return Flm_pivots(x, pp, rr, 1);
1604 : }
1605 : }
1606 :
1607 : static GEN
1608 21 : F2xqM_gauss_pivot(GEN x, GEN T, long *rr)
1609 : {
1610 : void *E;
1611 21 : const struct bb_field *S = get_F2xq_field(&E,T);
1612 21 : return gen_pivots(x, rr, E, S, _F2xqM_mul);
1613 : }
1614 :
1615 : static GEN
1616 1407 : FlxqM_gauss_pivot(GEN x, GEN T, ulong p, long *rr) {
1617 : void *E;
1618 1407 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1619 1407 : return gen_pivots(x, rr, E, S, _FlxqM_mul);
1620 : }
1621 :
1622 : static GEN
1623 105 : FqM_gauss_pivot_gen(GEN x, GEN T, GEN p, long *rr)
1624 : {
1625 : void *E;
1626 105 : const struct bb_field *S = get_Fq_field(&E,T,p);
1627 105 : return gen_pivots(x, rr, E, S, _FqM_mul);
1628 : }
1629 : static GEN
1630 1484 : FqM_gauss_pivot(GEN x, GEN T, GEN p, long *rr)
1631 : {
1632 1484 : if (lg(x)==1) { *rr = 0; return NULL; }
1633 1484 : if (!T) return FpM_gauss_pivot(x, p, rr);
1634 1484 : if (lgefint(p) == 3)
1635 : {
1636 1379 : pari_sp av = avma;
1637 1379 : ulong pp = uel(p,2);
1638 1379 : GEN Tp = ZXT_to_FlxT(T, pp);
1639 1379 : GEN d = FlxqM_gauss_pivot(ZXM_to_FlxM(x, pp, get_Flx_var(Tp)), Tp, pp, rr);
1640 1379 : return d ? gerepileuptoleaf(av, d): d;
1641 : }
1642 105 : return FqM_gauss_pivot_gen(x, T, p, rr);
1643 : }
1644 :
1645 : GEN
1646 289642 : FpM_image(GEN x, GEN p)
1647 : {
1648 : long r;
1649 289642 : GEN d = FpM_gauss_pivot(x,p,&r); /* d left on stack for efficiency */
1650 289640 : return image_from_pivot(x,d,r);
1651 : }
1652 :
1653 : GEN
1654 27415 : Flm_image(GEN x, ulong p)
1655 : {
1656 : long r;
1657 27415 : GEN d = Flm_pivots(x, p, &r, 0); /* d left on stack for efficiency */
1658 27415 : return image_from_pivot(x,d,r);
1659 : }
1660 :
1661 : GEN
1662 7 : F2m_image(GEN x)
1663 : {
1664 : long r;
1665 7 : GEN d = F2m_gauss_pivot(F2m_copy(x),&r); /* d left on stack for efficiency */
1666 7 : return image_from_pivot(x,d,r);
1667 : }
1668 :
1669 : GEN
1670 7 : F2xqM_image(GEN x, GEN T)
1671 : {
1672 : long r;
1673 7 : GEN d = F2xqM_gauss_pivot(x,T,&r); /* d left on stack for efficiency */
1674 7 : return image_from_pivot(x,d,r);
1675 : }
1676 :
1677 : GEN
1678 21 : FlxqM_image(GEN x, GEN T, ulong p)
1679 : {
1680 : long r;
1681 21 : GEN d = FlxqM_gauss_pivot(x, T, p, &r); /* d left on stack for efficiency */
1682 21 : return image_from_pivot(x,d,r);
1683 : }
1684 :
1685 : GEN
1686 49 : FqM_image(GEN x, GEN T, GEN p)
1687 : {
1688 : long r;
1689 49 : GEN d = FqM_gauss_pivot(x,T,p,&r); /* d left on stack for efficiency */
1690 49 : return image_from_pivot(x,d,r);
1691 : }
1692 :
1693 : long
1694 56 : FpM_rank(GEN x, GEN p)
1695 : {
1696 56 : pari_sp av = avma;
1697 : long r;
1698 56 : (void)FpM_gauss_pivot(x,p,&r);
1699 56 : return gc_long(av, lg(x)-1 - r);
1700 : }
1701 :
1702 : long
1703 7 : F2xqM_rank(GEN x, GEN T)
1704 : {
1705 7 : pari_sp av = avma;
1706 : long r;
1707 7 : (void)F2xqM_gauss_pivot(x,T,&r);
1708 7 : return gc_long(av, lg(x)-1 - r);
1709 : }
1710 :
1711 : long
1712 35 : FlxqM_rank(GEN x, GEN T, ulong p)
1713 : {
1714 : void *E;
1715 35 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1716 35 : return gen_matrank(x, E, S, _FlxqM_mul);
1717 : }
1718 :
1719 : long
1720 70 : FqM_rank(GEN x, GEN T, GEN p)
1721 : {
1722 70 : pari_sp av = avma;
1723 : long r;
1724 70 : (void)FqM_gauss_pivot(x,T,p,&r);
1725 70 : return gc_long(av, lg(x)-1 - r);
1726 : }
1727 :
1728 : static GEN
1729 35 : FpM_invimage_gen(GEN A, GEN B, GEN p)
1730 : {
1731 : void *E;
1732 35 : const struct bb_field *ff = get_Fp_field(&E, p);
1733 35 : return gen_invimage(A, B, E, ff, _FpM_mul);
1734 : }
1735 :
1736 : GEN
1737 0 : FpM_invimage(GEN A, GEN B, GEN p)
1738 : {
1739 0 : pari_sp av = avma;
1740 : ulong pp;
1741 : GEN y;
1742 :
1743 0 : A = FpM_init(A, p, &pp);
1744 0 : switch(pp)
1745 : {
1746 0 : case 0: return FpM_invimage_gen(A, B, p);
1747 0 : case 2:
1748 0 : y = F2m_invimage(A, ZM_to_F2m(B));
1749 0 : if (!y) return gc_NULL(av);
1750 0 : y = F2m_to_ZM(y);
1751 0 : return gerepileupto(av, y);
1752 0 : default:
1753 0 : y = Flm_invimage_i(A, ZM_to_Flm(B, pp), pp);
1754 0 : if (!y) return gc_NULL(av);
1755 0 : y = Flm_to_ZM(y);
1756 0 : return gerepileupto(av, y);
1757 : }
1758 : }
1759 :
1760 : GEN
1761 21 : F2xqM_invimage(GEN A, GEN B, GEN T) {
1762 : void *E;
1763 21 : const struct bb_field *ff = get_F2xq_field(&E, T);
1764 21 : return gen_invimage(A, B, E, ff, _F2xqM_mul);
1765 : }
1766 :
1767 : GEN
1768 42 : FlxqM_invimage(GEN A, GEN B, GEN T, ulong p) {
1769 : void *E;
1770 42 : const struct bb_field *ff = get_Flxq_field(&E, T, p);
1771 42 : return gen_invimage(A, B, E, ff, _FlxqM_mul);
1772 : }
1773 :
1774 : GEN
1775 42 : FqM_invimage(GEN A, GEN B, GEN T, GEN p) {
1776 : void *E;
1777 42 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1778 42 : return gen_invimage(A, B, E, ff, _FqM_mul);
1779 : }
1780 :
1781 : static GEN
1782 7 : FpM_FpC_invimage_gen(GEN A, GEN y, GEN p)
1783 : {
1784 : void *E;
1785 7 : const struct bb_field *ff = get_Fp_field(&E, p);
1786 7 : return gen_matcolinvimage_i(A, y, E, ff, _FpM_mul);
1787 : }
1788 :
1789 : GEN
1790 283758 : FpM_FpC_invimage(GEN A, GEN x, GEN p)
1791 : {
1792 283758 : pari_sp av = avma;
1793 : ulong pp;
1794 : GEN y;
1795 :
1796 283758 : A = FpM_init(A, p, &pp);
1797 283766 : switch(pp)
1798 : {
1799 7 : case 0: return FpM_FpC_invimage_gen(A, x, p);
1800 187760 : case 2:
1801 187760 : y = F2m_F2c_invimage(A, ZV_to_F2v(x));
1802 187761 : if (!y) return y;
1803 187761 : y = F2c_to_ZC(y);
1804 187759 : return gerepileupto(av, y);
1805 95999 : default:
1806 95999 : y = Flm_Flc_invimage(A, ZV_to_Flv(x, pp), pp);
1807 95999 : if (!y) return y;
1808 95999 : y = Flc_to_ZC(y);
1809 95998 : return gerepileupto(av, y);
1810 : }
1811 : }
1812 :
1813 : GEN
1814 21 : F2xqM_F2xqC_invimage(GEN A, GEN B, GEN T) {
1815 : void *E;
1816 21 : const struct bb_field *ff = get_F2xq_field(&E, T);
1817 21 : return gen_matcolinvimage_i(A, B, E, ff, _F2xqM_mul);
1818 : }
1819 :
1820 : GEN
1821 21 : FlxqM_FlxqC_invimage(GEN A, GEN B, GEN T, ulong p) {
1822 : void *E;
1823 21 : const struct bb_field *ff = get_Flxq_field(&E, T, p);
1824 21 : return gen_matcolinvimage_i(A, B, E, ff, _FlxqM_mul);
1825 : }
1826 :
1827 : GEN
1828 21 : FqM_FqC_invimage(GEN A, GEN B, GEN T, GEN p) {
1829 : void *E;
1830 21 : const struct bb_field *ff = get_Fq_field(&E, T, p);
1831 21 : return gen_matcolinvimage_i(A, B, E, ff, _FqM_mul);
1832 : }
1833 :
1834 : static GEN
1835 2250 : FpM_ker_gen(GEN x, GEN p, long deplin)
1836 : {
1837 : void *E;
1838 2250 : const struct bb_field *S = get_Fp_field(&E,p);
1839 2250 : return gen_ker_i(x, deplin, E, S, _FpM_mul);
1840 : }
1841 : static GEN
1842 1059256 : FpM_ker_i(GEN x, GEN p, long deplin)
1843 : {
1844 1059256 : pari_sp av = avma;
1845 : ulong pp;
1846 : GEN y;
1847 :
1848 1059256 : if (lg(x)==1) return cgetg(1,t_MAT);
1849 1059256 : x = FpM_init(x, p, &pp);
1850 1059262 : switch(pp)
1851 : {
1852 2180 : case 0: return FpM_ker_gen(x,p,deplin);
1853 669036 : case 2:
1854 669036 : y = F2m_ker_sp(x, deplin);
1855 669049 : if (!y) return gc_NULL(av);
1856 669055 : y = deplin? F2c_to_ZC(y): F2m_to_ZM(y);
1857 669054 : return gerepileupto(av, y);
1858 388046 : default:
1859 388046 : y = Flm_ker_sp(x, pp, deplin);
1860 388045 : if (!y) return gc_NULL(av);
1861 388045 : y = deplin? Flc_to_ZC(y): Flm_to_ZM(y);
1862 388046 : return gerepileupto(av, y);
1863 : }
1864 : }
1865 :
1866 : GEN
1867 624233 : FpM_ker(GEN x, GEN p) { return FpM_ker_i(x,p,0); }
1868 :
1869 : static GEN
1870 35 : F2xqM_ker_i(GEN x, GEN T, long deplin)
1871 : {
1872 : const struct bb_field *ff;
1873 : void *E;
1874 :
1875 35 : if (lg(x)==1) return cgetg(1,t_MAT);
1876 35 : ff = get_F2xq_field(&E,T);
1877 35 : return gen_ker_i(x,deplin, E, ff, _F2xqM_mul);
1878 : }
1879 :
1880 : GEN
1881 21 : F2xqM_ker(GEN x, GEN T)
1882 : {
1883 21 : return F2xqM_ker_i(x, T, 0);
1884 : }
1885 :
1886 : static GEN
1887 2926 : FlxqM_ker_i(GEN x, GEN T, ulong p, long deplin) {
1888 : void *E;
1889 2926 : const struct bb_field *S = get_Flxq_field(&E, T, p);
1890 2926 : return gen_ker_i(x, deplin, E, S, _FlxqM_mul);
1891 : }
1892 :
1893 : GEN
1894 2891 : FlxqM_ker(GEN x, GEN T, ulong p)
1895 : {
1896 2891 : return FlxqM_ker_i(x, T, p, 0);
1897 : }
1898 :
1899 : static GEN
1900 126 : FqM_ker_gen(GEN x, GEN T, GEN p, long deplin)
1901 : {
1902 : void *E;
1903 126 : const struct bb_field *S = get_Fq_field(&E,T,p);
1904 126 : return gen_ker_i(x,deplin,E,S,_FqM_mul);
1905 : }
1906 : static GEN
1907 8981 : FqM_ker_i(GEN x, GEN T, GEN p, long deplin)
1908 : {
1909 8981 : if (!T) return FpM_ker_i(x,p,deplin);
1910 2989 : if (lg(x)==1) return cgetg(1,t_MAT);
1911 :
1912 2989 : if (lgefint(p)==3)
1913 : {
1914 2863 : pari_sp ltop=avma;
1915 2863 : ulong l= p[2];
1916 2863 : GEN Tl = ZXT_to_FlxT(T,l);
1917 2863 : GEN Ml = ZXM_to_FlxM(x, l, get_Flx_var(Tl));
1918 2863 : GEN p1 = FlxM_to_ZXM(FlxqM_ker(Ml,Tl,l));
1919 2863 : return gerepileupto(ltop,p1);
1920 : }
1921 126 : return FqM_ker_gen(x, T, p, deplin);
1922 : }
1923 :
1924 : GEN
1925 8904 : FqM_ker(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,0); }
1926 :
1927 : GEN
1928 429045 : FpM_deplin(GEN x, GEN p) { return FpM_ker_i(x,p,1); }
1929 :
1930 : GEN
1931 14 : F2xqM_deplin(GEN x, GEN T)
1932 : {
1933 14 : return F2xqM_ker_i(x, T, 1);
1934 : }
1935 :
1936 : GEN
1937 35 : FlxqM_deplin(GEN x, GEN T, ulong p)
1938 : {
1939 35 : return FlxqM_ker_i(x, T, p, 1);
1940 : }
1941 :
1942 : GEN
1943 77 : FqM_deplin(GEN x, GEN T, GEN p) { return FqM_ker_i(x,T,p,1); }
1944 :
1945 : static GEN
1946 3522 : FpM_gauss_gen(GEN a, GEN b, GEN p)
1947 : {
1948 : void *E;
1949 3522 : const struct bb_field *S = get_Fp_field(&E,p);
1950 3522 : return gen_gauss(a,b, E, S, _FpM_mul);
1951 : }
1952 : /* a an FpM, lg(a)>1; b an FpM or NULL (replace by identity) */
1953 : static GEN
1954 311549 : FpM_gauss_i(GEN a, GEN b, GEN p, ulong *pp)
1955 : {
1956 311549 : long n = nbrows(a);
1957 311547 : a = FpM_init(a,p,pp);
1958 311544 : switch(*pp)
1959 : {
1960 3522 : case 0:
1961 3522 : if (!b) b = matid(n);
1962 3522 : return FpM_gauss_gen(a,b,p);
1963 219039 : case 2:
1964 219039 : if (b) b = ZM_to_F2m(b); else b = matid_F2m(n);
1965 219036 : return F2m_gauss_sp(a,b);
1966 88983 : default:
1967 88983 : if (b) b = ZM_to_Flm(b, *pp); else b = matid_Flm(n);
1968 88983 : return Flm_gauss_sp(a,b, NULL, *pp);
1969 : }
1970 : }
1971 : GEN
1972 35 : FpM_gauss(GEN a, GEN b, GEN p)
1973 : {
1974 35 : pari_sp av = avma;
1975 : ulong pp;
1976 : GEN u;
1977 35 : if (lg(a) == 1 || lg(b)==1) return cgetg(1, t_MAT);
1978 35 : u = FpM_gauss_i(a, b, p, &pp);
1979 35 : if (!u) return gc_NULL(av);
1980 28 : switch(pp)
1981 : {
1982 28 : case 0: return gerepilecopy(av, u);
1983 0 : case 2: u = F2m_to_ZM(u); break;
1984 0 : default: u = Flm_to_ZM(u); break;
1985 : }
1986 0 : return gerepileupto(av, u);
1987 : }
1988 :
1989 : static GEN
1990 84 : F2xqM_gauss_gen(GEN a, GEN b, GEN T)
1991 : {
1992 : void *E;
1993 84 : const struct bb_field *S = get_F2xq_field(&E, T);
1994 84 : return gen_gauss(a, b, E, S, _F2xqM_mul);
1995 : }
1996 :
1997 : GEN
1998 21 : F2xqM_gauss(GEN a, GEN b, GEN T)
1999 : {
2000 21 : pari_sp av = avma;
2001 21 : long n = lg(a)-1;
2002 : GEN u;
2003 21 : if (!n || lg(b)==1) { set_avma(av); return cgetg(1, t_MAT); }
2004 21 : u = F2xqM_gauss_gen(a, b, T);
2005 21 : if (!u) return gc_NULL(av);
2006 14 : return gerepilecopy(av, u);
2007 : }
2008 :
2009 : static GEN
2010 91 : FlxqM_gauss_i(GEN a, GEN b, GEN T, ulong p) {
2011 : void *E;
2012 91 : const struct bb_field *S = get_Flxq_field(&E, T, p);
2013 91 : return gen_gauss(a, b, E, S, _FlxqM_mul);
2014 : }
2015 :
2016 : GEN
2017 21 : FlxqM_gauss(GEN a, GEN b, GEN T, ulong p)
2018 : {
2019 21 : pari_sp av = avma;
2020 21 : long n = lg(a)-1;
2021 : GEN u;
2022 21 : if (!n || lg(b)==1) { set_avma(av); return cgetg(1, t_MAT); }
2023 21 : u = FlxqM_gauss_i(a, b, T, p);
2024 21 : if (!u) return gc_NULL(av);
2025 14 : return gerepilecopy(av, u);
2026 : }
2027 :
2028 : static GEN
2029 133 : FqM_gauss_gen(GEN a, GEN b, GEN T, GEN p)
2030 : {
2031 : void *E;
2032 133 : const struct bb_field *S = get_Fq_field(&E,T,p);
2033 133 : return gen_gauss(a,b,E,S,_FqM_mul);
2034 : }
2035 : GEN
2036 21 : FqM_gauss(GEN a, GEN b, GEN T, GEN p)
2037 : {
2038 21 : pari_sp av = avma;
2039 : GEN u;
2040 : long n;
2041 21 : if (!T) return FpM_gauss(a,b,p);
2042 21 : n = lg(a)-1; if (!n || lg(b)==1) return cgetg(1, t_MAT);
2043 21 : u = FqM_gauss_gen(a,b,T,p);
2044 21 : if (!u) return gc_NULL(av);
2045 14 : return gerepilecopy(av, u);
2046 : }
2047 :
2048 : GEN
2049 14 : FpM_FpC_gauss(GEN a, GEN b, GEN p)
2050 : {
2051 14 : pari_sp av = avma;
2052 : ulong pp;
2053 : GEN u;
2054 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2055 14 : u = FpM_gauss_i(a, mkmat(b), p, &pp);
2056 14 : if (!u) return gc_NULL(av);
2057 14 : switch(pp)
2058 : {
2059 14 : case 0: return gerepilecopy(av, gel(u,1));
2060 0 : case 2: u = F2c_to_ZC(gel(u,1)); break;
2061 0 : default: u = Flc_to_ZC(gel(u,1)); break;
2062 : }
2063 0 : return gerepileupto(av, u);
2064 : }
2065 :
2066 : GEN
2067 28 : F2xqM_F2xqC_gauss(GEN a, GEN b, GEN T)
2068 : {
2069 28 : pari_sp av = avma;
2070 : GEN u;
2071 28 : if (lg(a) == 1) return cgetg(1, t_COL);
2072 28 : u = F2xqM_gauss_gen(a, mkmat(b), T);
2073 28 : if (!u) return gc_NULL(av);
2074 14 : return gerepilecopy(av, gel(u,1));
2075 : }
2076 :
2077 : GEN
2078 14 : FlxqM_FlxqC_gauss(GEN a, GEN b, GEN T, ulong p)
2079 : {
2080 14 : pari_sp av = avma;
2081 : GEN u;
2082 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2083 14 : u = FlxqM_gauss_i(a, mkmat(b), T, p);
2084 14 : if (!u) return gc_NULL(av);
2085 7 : return gerepilecopy(av, gel(u,1));
2086 : }
2087 :
2088 : GEN
2089 14 : FqM_FqC_gauss(GEN a, GEN b, GEN T, GEN p)
2090 : {
2091 14 : pari_sp av = avma;
2092 : GEN u;
2093 14 : if (!T) return FpM_FpC_gauss(a,b,p);
2094 14 : if (lg(a) == 1) return cgetg(1, t_COL);
2095 14 : u = FqM_gauss_gen(a,mkmat(b),T,p);
2096 14 : if (!u) return gc_NULL(av);
2097 7 : return gerepilecopy(av, gel(u,1));
2098 : }
2099 :
2100 : GEN
2101 311500 : FpM_inv(GEN a, GEN p)
2102 : {
2103 311500 : pari_sp av = avma;
2104 : ulong pp;
2105 : GEN u;
2106 311500 : if (lg(a) == 1) return cgetg(1, t_MAT);
2107 311500 : u = FpM_gauss_i(a, NULL, p, &pp);
2108 311489 : if (!u) return gc_NULL(av);
2109 311475 : switch(pp)
2110 : {
2111 3466 : case 0: return gerepilecopy(av, u);
2112 219026 : case 2: u = F2m_to_ZM(u); break;
2113 88983 : default: u = Flm_to_ZM(u); break;
2114 : }
2115 308013 : return gerepileupto(av, u);
2116 : }
2117 :
2118 : GEN
2119 35 : F2xqM_inv(GEN a, GEN T)
2120 : {
2121 35 : pari_sp av = avma;
2122 : GEN u;
2123 35 : if (lg(a) == 1) { set_avma(av); return cgetg(1, t_MAT); }
2124 35 : u = F2xqM_gauss_gen(a, matid_F2xqM(nbrows(a),T), T);
2125 35 : if (!u) return gc_NULL(av);
2126 28 : return gerepilecopy(av, u);
2127 : }
2128 :
2129 : GEN
2130 56 : FlxqM_inv(GEN a, GEN T, ulong p)
2131 : {
2132 56 : pari_sp av = avma;
2133 : GEN u;
2134 56 : if (lg(a) == 1) { set_avma(av); return cgetg(1, t_MAT); }
2135 56 : u = FlxqM_gauss_i(a, matid_FlxqM(nbrows(a),T,p), T,p);
2136 56 : if (!u) return gc_NULL(av);
2137 42 : return gerepilecopy(av, u);
2138 : }
2139 :
2140 : GEN
2141 98 : FqM_inv(GEN a, GEN T, GEN p)
2142 : {
2143 98 : pari_sp av = avma;
2144 : GEN u;
2145 98 : if (!T) return FpM_inv(a,p);
2146 98 : if (lg(a) == 1) return cgetg(1, t_MAT);
2147 98 : u = FqM_gauss_gen(a,matid(nbrows(a)),T,p);
2148 98 : if (!u) return gc_NULL(av);
2149 70 : return gerepilecopy(av, u);
2150 : }
2151 :
2152 : GEN
2153 283233 : FpM_intersect_i(GEN x, GEN y, GEN p)
2154 : {
2155 283233 : long j, lx = lg(x);
2156 : GEN z;
2157 :
2158 283233 : if (lx == 1 || lg(y) == 1) return cgetg(1,t_MAT);
2159 283233 : if (lgefint(p) == 3)
2160 : {
2161 283233 : ulong pp = p[2];
2162 283233 : return Flm_to_ZM(Flm_intersect_i(ZM_to_Flm(x,pp), ZM_to_Flm(y,pp), pp));
2163 : }
2164 0 : z = FpM_ker(shallowconcat(x,y), p);
2165 0 : for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
2166 0 : return FpM_mul(x,z,p);
2167 : }
2168 : GEN
2169 0 : FpM_intersect(GEN x, GEN y, GEN p)
2170 : {
2171 0 : pari_sp av = avma;
2172 : GEN z;
2173 0 : if (lgefint(p) == 3)
2174 : {
2175 0 : ulong pp = p[2];
2176 0 : z = Flm_image(Flm_intersect_i(ZM_to_Flm(x,pp), ZM_to_Flm(y,pp), pp), pp);
2177 : }
2178 : else
2179 0 : z = FpM_image(FpM_intersect_i(x,y,p), p);
2180 0 : return gerepileupto(av, z);
2181 : }
2182 :
2183 : static void
2184 244035 : init_suppl(GEN x)
2185 : {
2186 244035 : if (lg(x) == 1) pari_err_IMPL("suppl [empty matrix]");
2187 : /* HACK: avoid overwriting d from gauss_pivot() after set_avma(av) */
2188 244035 : (void)new_chunk(lgcols(x) * 2);
2189 244034 : }
2190 :
2191 : GEN
2192 242488 : FpM_suppl(GEN x, GEN p)
2193 : {
2194 : GEN d;
2195 : long r;
2196 242488 : init_suppl(x); d = FpM_gauss_pivot(x,p, &r);
2197 242486 : return get_suppl(x,d,nbrows(x),r,&col_ei);
2198 : }
2199 :
2200 : GEN
2201 14 : F2m_suppl(GEN x)
2202 : {
2203 : GEN d;
2204 : long r;
2205 14 : init_suppl(x); d = F2m_gauss_pivot(F2m_copy(x), &r);
2206 14 : return get_suppl(x,d,mael(x,1,1),r,&F2v_ei);
2207 : }
2208 :
2209 : GEN
2210 105 : Flm_suppl(GEN x, ulong p)
2211 : {
2212 : GEN d;
2213 : long r;
2214 105 : init_suppl(x); d = Flm_pivots(x, p, &r, 0);
2215 105 : return get_suppl(x,d,nbrows(x),r,&vecsmall_ei);
2216 : }
2217 :
2218 : GEN
2219 7 : F2xqM_suppl(GEN x, GEN T)
2220 : {
2221 : void *E;
2222 7 : const struct bb_field *S = get_F2xq_field(&E, T);
2223 7 : return gen_suppl(x, E, S, _F2xqM_mul);
2224 : }
2225 :
2226 : GEN
2227 14 : FlxqM_suppl(GEN x, GEN T, ulong p)
2228 : {
2229 : void *E;
2230 14 : const struct bb_field *S = get_Flxq_field(&E, T, p);
2231 14 : return gen_suppl(x, E, S, _FlxqM_mul);
2232 : }
2233 :
2234 : GEN
2235 4277 : FqM_suppl(GEN x, GEN T, GEN p)
2236 : {
2237 4277 : pari_sp av = avma;
2238 : GEN d;
2239 : long r;
2240 :
2241 4277 : if (!T) return FpM_suppl(x,p);
2242 1358 : init_suppl(x);
2243 1358 : d = FqM_gauss_pivot(x,T,p,&r);
2244 1358 : set_avma(av); return get_suppl(x,d,nbrows(x),r,&col_ei);
2245 : }
2246 :
2247 : static void
2248 173031 : init_indexrank(GEN x) {
2249 173031 : (void)new_chunk(3 + 2*lg(x)); /* HACK */
2250 173031 : }
2251 :
2252 : GEN
2253 30282 : FpM_indexrank(GEN x, GEN p) {
2254 30282 : pari_sp av = avma;
2255 : long r;
2256 : GEN d;
2257 30282 : init_indexrank(x);
2258 30282 : d = FpM_gauss_pivot(x,p,&r);
2259 30282 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2260 : }
2261 :
2262 : GEN
2263 52976 : Flm_indexrank(GEN x, ulong p) {
2264 52976 : pari_sp av = avma;
2265 : long r;
2266 : GEN d;
2267 52976 : init_indexrank(x);
2268 52976 : d = Flm_pivots(x, p, &r, 0);
2269 52976 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2270 : }
2271 :
2272 : GEN
2273 48 : F2m_indexrank(GEN x) {
2274 48 : pari_sp av = avma;
2275 : long r;
2276 : GEN d;
2277 48 : init_indexrank(x);
2278 48 : d = F2m_gauss_pivot(F2m_copy(x),&r);
2279 48 : set_avma(av); return indexrank0(lg(x)-1, r, d);
2280 : }
2281 :
2282 : GEN
2283 7 : F2xqM_indexrank(GEN x, GEN T) {
2284 7 : pari_sp av = avma;
2285 : long r;
2286 : GEN d;
2287 7 : init_indexrank(x);
2288 7 : d = F2xqM_gauss_pivot(x, T, &r);
2289 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2290 : }
2291 :
2292 : GEN
2293 7 : FlxqM_indexrank(GEN x, GEN T, ulong p) {
2294 7 : pari_sp av = avma;
2295 : long r;
2296 : GEN d;
2297 7 : init_indexrank(x);
2298 7 : d = FlxqM_gauss_pivot(x, T, p, &r);
2299 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2300 : }
2301 :
2302 : GEN
2303 7 : FqM_indexrank(GEN x, GEN T, GEN p) {
2304 7 : pari_sp av = avma;
2305 : long r;
2306 : GEN d;
2307 7 : init_indexrank(x);
2308 7 : d = FqM_gauss_pivot(x, T, p, &r);
2309 7 : set_avma(av); return indexrank0(lg(x) - 1, r, d);
2310 : }
2311 :
2312 : /*******************************************************************/
2313 : /* */
2314 : /* Solve A*X=B (Gauss pivot) */
2315 : /* */
2316 : /*******************************************************************/
2317 : /* x ~ 0 compared to reference y */
2318 : int
2319 1116592 : approx_0(GEN x, GEN y)
2320 : {
2321 1116592 : long tx = typ(x);
2322 1116592 : if (tx == t_COMPLEX)
2323 140 : return approx_0(gel(x,1), y) && approx_0(gel(x,2), y);
2324 1904908 : return gequal0(x) ||
2325 788455 : (tx == t_REAL && gexpo(y) - gexpo(x) > bit_prec(x));
2326 : }
2327 : /* x a column, x0 same column in the original input matrix (for reference),
2328 : * c list of pivots so far */
2329 : static long
2330 1183704 : gauss_get_pivot_max(GEN X, GEN X0, long ix, GEN c)
2331 : {
2332 1183704 : GEN p, r, x = gel(X,ix), x0 = gel(X0,ix);
2333 1183704 : long i, k = 0, ex = - (long)HIGHEXPOBIT, lx = lg(x);
2334 1183704 : if (c)
2335 : {
2336 312514 : for (i=1; i<lx; i++)
2337 178557 : if (!c[i])
2338 : {
2339 86180 : long e = gexpo(gel(x,i));
2340 86180 : if (e > ex) { ex = e; k = i; }
2341 : }
2342 : }
2343 : else
2344 : {
2345 3848096 : for (i=ix; i<lx; i++)
2346 : {
2347 2798322 : long e = gexpo(gel(x,i));
2348 2798349 : if (e > ex) { ex = e; k = i; }
2349 : }
2350 : }
2351 1183731 : if (!k) return lx;
2352 1116422 : p = gel(x,k);
2353 1116422 : r = gel(x0,k); if (isrationalzero(r)) r = x0;
2354 1116417 : return approx_0(p, r)? lx: k;
2355 : }
2356 : static long
2357 252357 : gauss_get_pivot_padic(GEN X, GEN p, long ix, GEN c)
2358 : {
2359 252357 : GEN x = gel(X, ix);
2360 252357 : long i, k = 0, ex = (long)HIGHVALPBIT, lx = lg(x);
2361 252357 : if (c)
2362 : {
2363 504 : for (i=1; i<lx; i++)
2364 378 : if (!c[i] && !gequal0(gel(x,i)))
2365 : {
2366 245 : long e = gvaluation(gel(x,i), p);
2367 245 : if (e < ex) { ex = e; k = i; }
2368 : }
2369 : }
2370 : else
2371 : {
2372 2311841 : for (i=ix; i<lx; i++)
2373 2059610 : if (!gequal0(gel(x,i)))
2374 : {
2375 1561196 : long e = gvaluation(gel(x,i), p);
2376 1561196 : if (e < ex) { ex = e; k = i; }
2377 : }
2378 : }
2379 252357 : return k? k: lx;
2380 : }
2381 : static long
2382 4193 : gauss_get_pivot_NZ(GEN X, GEN x0/*unused*/, long ix, GEN c)
2383 : {
2384 4193 : GEN x = gel(X, ix);
2385 4193 : long i, lx = lg(x);
2386 : (void)x0;
2387 4193 : if (c)
2388 : {
2389 11606 : for (i=1; i<lx; i++)
2390 10752 : if (!c[i] && !gequal0(gel(x,i))) return i;
2391 : }
2392 : else
2393 : {
2394 2380 : for (i=ix; i<lx; i++)
2395 2366 : if (!gequal0(gel(x,i))) return i;
2396 : }
2397 868 : return lx;
2398 : }
2399 :
2400 : /* Return pivot seeking function appropriate for the domain of the RgM x
2401 : * (first non zero pivot, maximal pivot...)
2402 : * x0 is a reference point used when guessing whether x[i,j] ~ 0
2403 : * (iff x[i,j] << x0[i,j]); typical case: mateigen, Gauss pivot on x - vp.Id,
2404 : * but use original x when deciding whether a prospective pivot is nonzero */
2405 : static pivot_fun
2406 410416 : get_pivot_fun(GEN x, GEN x0, GEN *data)
2407 : {
2408 410416 : long i, j, hx, lx = lg(x);
2409 410416 : int res = t_INT;
2410 410416 : GEN p = NULL;
2411 :
2412 410416 : *data = NULL;
2413 410416 : if (lx == 1) return &gauss_get_pivot_NZ;
2414 410381 : hx = lgcols(x);
2415 1907740 : for (j=1; j<lx; j++)
2416 : {
2417 1497402 : GEN xj = gel(x,j);
2418 10178668 : for (i=1; i<hx; i++)
2419 : {
2420 8681308 : GEN c = gel(xj,i);
2421 8681308 : switch(typ(c))
2422 : {
2423 3517504 : case t_REAL:
2424 3517504 : res = t_REAL;
2425 3517504 : break;
2426 364 : case t_COMPLEX:
2427 364 : if (typ(gel(c,1)) == t_REAL || typ(gel(c,2)) == t_REAL) res = t_REAL;
2428 364 : break;
2429 3610581 : case t_INT: case t_INTMOD: case t_FRAC: case t_FFELT: case t_QUAD:
2430 : case t_POLMOD: /* exact types */
2431 3610581 : break;
2432 1552817 : case t_PADIC:
2433 1552817 : p = gel(c,2);
2434 1552817 : res = t_PADIC;
2435 1552817 : break;
2436 42 : default: return &gauss_get_pivot_NZ;
2437 : }
2438 : }
2439 : }
2440 410338 : switch(res)
2441 : {
2442 380674 : case t_REAL: *data = x0; return &gauss_get_pivot_max;
2443 28399 : case t_PADIC: *data = p; return &gauss_get_pivot_padic;
2444 1265 : default: return &gauss_get_pivot_NZ;
2445 : }
2446 : }
2447 :
2448 : static GEN
2449 467849 : get_col(GEN a, GEN b, GEN p, long li)
2450 : {
2451 467849 : GEN u = cgetg(li+1,t_COL);
2452 : long i, j;
2453 :
2454 467849 : gel(u,li) = gdiv(gel(b,li), p);
2455 2079999 : for (i=li-1; i>0; i--)
2456 : {
2457 1612153 : pari_sp av = avma;
2458 1612153 : GEN m = gel(b,i);
2459 7909848 : for (j=i+1; j<=li; j++) m = gsub(m, gmul(gcoeff(a,i,j), gel(u,j)));
2460 1612133 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(a,i,i)));
2461 : }
2462 467846 : return u;
2463 : }
2464 :
2465 : /* bk -= m * bi */
2466 : static void
2467 7433105 : _submul(GEN b, long k, long i, GEN m)
2468 : {
2469 7433105 : gel(b,k) = gsub(gel(b,k), gmul(m, gel(b,i)));
2470 7433001 : }
2471 : static int
2472 1031893 : init_gauss(GEN a, GEN *b, long *aco, long *li, int *iscol)
2473 : {
2474 1031893 : *iscol = *b ? (typ(*b) == t_COL): 0;
2475 1031893 : *aco = lg(a) - 1;
2476 1031893 : if (!*aco) /* a empty */
2477 : {
2478 70 : if (*b && lg(*b) != 1) pari_err_DIM("gauss");
2479 70 : *li = 0; return 0;
2480 : }
2481 1031823 : *li = nbrows(a);
2482 1031822 : if (*li < *aco) pari_err_INV("gauss [no left inverse]", a);
2483 1031824 : if (*b)
2484 : {
2485 955960 : switch(typ(*b))
2486 : {
2487 123346 : case t_MAT:
2488 123346 : if (lg(*b) == 1) return 0;
2489 123346 : *b = RgM_shallowcopy(*b);
2490 123348 : break;
2491 832615 : case t_COL:
2492 832615 : *b = mkmat( leafcopy(*b) );
2493 832616 : break;
2494 0 : default: pari_err_TYPE("gauss",*b);
2495 : }
2496 955964 : if (nbrows(*b) != *li) pari_err_DIM("gauss");
2497 : }
2498 : else
2499 75864 : *b = matid(*li);
2500 1031824 : return 1;
2501 : }
2502 :
2503 : static GEN
2504 112 : RgM_inv_FpM(GEN a, GEN p)
2505 : {
2506 : ulong pp;
2507 112 : a = RgM_Fp_init(a, p, &pp);
2508 112 : switch(pp)
2509 : {
2510 35 : case 0:
2511 35 : a = FpM_inv(a,p);
2512 35 : if (a) a = FpM_to_mod(a, p);
2513 35 : break;
2514 35 : case 2:
2515 35 : a = F2m_inv(a);
2516 35 : if (a) a = F2m_to_mod(a);
2517 35 : break;
2518 42 : default:
2519 42 : a = Flm_inv_sp(a, NULL, pp);
2520 42 : if (a) a = Flm_to_mod(a, pp);
2521 : }
2522 112 : return a;
2523 : }
2524 :
2525 : static GEN
2526 42 : RgM_inv_FqM(GEN x, GEN pol, GEN p)
2527 : {
2528 42 : pari_sp av = avma;
2529 42 : GEN b, T = RgX_to_FpX(pol, p);
2530 42 : if (signe(T) == 0) pari_err_OP("^",x,gen_m1);
2531 42 : b = FqM_inv(RgM_to_FqM(x, T, p), T, p);
2532 42 : if (!b) return gc_NULL(av);
2533 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
2534 : }
2535 :
2536 : #define code(t1,t2) ((t1 << 6) | t2)
2537 : static GEN
2538 197331 : RgM_inv_fast(GEN x)
2539 : {
2540 : GEN p, pol;
2541 : long pa;
2542 197331 : long t = RgM_type(x, &p,&pol,&pa);
2543 197331 : switch(t)
2544 : {
2545 44443 : case t_INT: /* Fall back */
2546 44443 : case t_FRAC: return QM_inv(x);
2547 147 : case t_FFELT: return FFM_inv(x, pol);
2548 112 : case t_INTMOD: return RgM_inv_FpM(x, p);
2549 41 : case code(t_POLMOD, t_INTMOD):
2550 41 : return RgM_inv_FqM(x, pol, p);
2551 152588 : default: return gen_0;
2552 : }
2553 : }
2554 : #undef code
2555 :
2556 : static GEN
2557 49 : RgM_RgC_solve_FpC(GEN a, GEN b, GEN p)
2558 : {
2559 49 : pari_sp av = avma;
2560 : ulong pp;
2561 49 : a = RgM_Fp_init(a, p, &pp);
2562 49 : switch(pp)
2563 : {
2564 14 : case 0:
2565 14 : b = RgC_to_FpC(b, p);
2566 14 : a = FpM_FpC_gauss(a,b,p);
2567 14 : return a ? gerepileupto(av, FpC_to_mod(a, p)): NULL;
2568 14 : case 2:
2569 14 : b = RgV_to_F2v(b);
2570 14 : a = F2m_F2c_gauss(a,b);
2571 14 : return a ? gerepileupto(av, F2c_to_mod(a)): NULL;
2572 21 : default:
2573 21 : b = RgV_to_Flv(b, pp);
2574 21 : a = Flm_Flc_gauss(a, b, pp);
2575 21 : return a ? gerepileupto(av, Flc_to_mod(a, pp)): NULL;
2576 : }
2577 : }
2578 :
2579 : static GEN
2580 98 : RgM_solve_FpM(GEN a, GEN b, GEN p)
2581 : {
2582 98 : pari_sp av = avma;
2583 : ulong pp;
2584 98 : a = RgM_Fp_init(a, p, &pp);
2585 98 : switch(pp)
2586 : {
2587 35 : case 0:
2588 35 : b = RgM_to_FpM(b, p);
2589 35 : a = FpM_gauss(a,b,p);
2590 35 : return a ? gerepileupto(av, FpM_to_mod(a, p)): NULL;
2591 21 : case 2:
2592 21 : b = RgM_to_F2m(b);
2593 21 : a = F2m_gauss(a,b);
2594 21 : return a ? gerepileupto(av, F2m_to_mod(a)): NULL;
2595 42 : default:
2596 42 : b = RgM_to_Flm(b, pp);
2597 42 : a = Flm_gauss(a,b,pp);
2598 42 : return a ? gerepileupto(av, Flm_to_mod(a, pp)): NULL;
2599 : }
2600 : }
2601 :
2602 : /* Gaussan Elimination. If a is square, return a^(-1)*b;
2603 : * if a has more rows than columns and b is NULL, return c such that c a = Id.
2604 : * a is a (not necessarily square) matrix
2605 : * b is a matrix or column vector, NULL meaning: take the identity matrix,
2606 : * effectively returning the inverse of a
2607 : * If a and b are empty, the result is the empty matrix.
2608 : *
2609 : * li: number of rows of a and b
2610 : * aco: number of columns of a
2611 : * bco: number of columns of b (if matrix)
2612 : */
2613 : static GEN
2614 553845 : RgM_solve_basecase(GEN a, GEN b)
2615 : {
2616 553845 : pari_sp av = avma;
2617 : long i, j, k, li, bco, aco;
2618 : int iscol;
2619 : pivot_fun pivot;
2620 : GEN p, u, data;
2621 :
2622 553845 : set_avma(av);
2623 :
2624 553845 : if (lg(a)-1 == 2 && nbrows(a) == 2) {
2625 : /* 2x2 matrix, start by inverting a */
2626 259376 : GEN u = gcoeff(a,1,1), v = gcoeff(a,1,2);
2627 259376 : GEN w = gcoeff(a,2,1), x = gcoeff(a,2,2);
2628 259376 : GEN D = gsub(gmul(u,x), gmul(v,w)), ainv;
2629 259370 : if (gequal0(D)) return NULL;
2630 259370 : ainv = mkmat2(mkcol2(x, gneg(w)), mkcol2(gneg(v), u));
2631 259376 : ainv = gmul(ainv, ginv(D));
2632 259368 : if (b) ainv = gmul(ainv, b);
2633 259366 : return gerepileupto(av, ainv);
2634 : }
2635 :
2636 294468 : if (!init_gauss(a, &b, &aco, &li, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
2637 294472 : pivot = get_pivot_fun(a, a, &data);
2638 294472 : a = RgM_shallowcopy(a);
2639 294473 : bco = lg(b)-1;
2640 294473 : if(DEBUGLEVEL>4) err_printf("Entering gauss\n");
2641 :
2642 294473 : p = NULL; /* gcc -Wall */
2643 978078 : for (i=1; i<=aco; i++)
2644 : {
2645 : /* k is the line where we find the pivot */
2646 978073 : k = pivot(a, data, i, NULL);
2647 978085 : if (k > li) return NULL;
2648 978065 : if (k != i)
2649 : { /* exchange the lines s.t. k = i */
2650 841175 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
2651 641110 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
2652 : }
2653 978065 : p = gcoeff(a,i,i);
2654 978065 : if (i == aco) break;
2655 :
2656 2280033 : for (k=i+1; k<=li; k++)
2657 : {
2658 1596432 : GEN m = gcoeff(a,k,i);
2659 1596432 : if (!gequal0(m))
2660 : {
2661 1197073 : m = gdiv(m,p);
2662 5470101 : for (j=i+1; j<=aco; j++) _submul(gel(a,j),k,i,m);
2663 4357404 : for (j=1; j<=bco; j++) _submul(gel(b,j),k,i,m);
2664 : }
2665 : }
2666 683601 : if (gc_needed(av,1))
2667 : {
2668 12 : if(DEBUGMEM>1) pari_warn(warnmem,"gauss. i=%ld",i);
2669 12 : gerepileall(av,2, &a,&b);
2670 : }
2671 : }
2672 :
2673 294457 : if(DEBUGLEVEL>4) err_printf("Solving the triangular system\n");
2674 294457 : u = cgetg(bco+1,t_MAT);
2675 762296 : for (j=1; j<=bco; j++) gel(u,j) = get_col(a,gel(b,j),p,aco);
2676 294447 : return gerepilecopy(av, iscol? gel(u,1): u);
2677 : }
2678 :
2679 : static GEN
2680 367235 : RgM_RgC_solve_fast(GEN x, GEN y)
2681 : {
2682 : GEN p, pol;
2683 : long pa;
2684 367235 : long t = RgM_RgC_type(x, y, &p,&pol,&pa);
2685 367234 : switch(t)
2686 : {
2687 13545 : case t_INT: return ZM_gauss(x, y);
2688 42 : case t_FRAC: return QM_gauss(x, y);
2689 49 : case t_INTMOD: return RgM_RgC_solve_FpC(x, y, p);
2690 56 : case t_FFELT: return FFM_FFC_gauss(x, y, pol);
2691 353542 : default: return gen_0;
2692 : }
2693 : }
2694 :
2695 : static GEN
2696 47934 : RgM_solve_fast(GEN x, GEN y)
2697 : {
2698 : GEN p, pol;
2699 : long pa;
2700 47934 : long t = RgM_type2(x, y, &p,&pol,&pa);
2701 47934 : switch(t)
2702 : {
2703 42 : case t_INT: return ZM_gauss(x, y);
2704 14 : case t_FRAC: return QM_gauss(x, y);
2705 98 : case t_INTMOD: return RgM_solve_FpM(x, y, p);
2706 63 : case t_FFELT: return FFM_gauss(x, y, pol);
2707 47717 : default: return gen_0;
2708 : }
2709 : }
2710 :
2711 : GEN
2712 415170 : RgM_solve(GEN a, GEN b)
2713 : {
2714 415170 : pari_sp av = avma;
2715 : GEN u;
2716 415170 : if (!b) return RgM_inv(a);
2717 415170 : u = typ(b)==t_MAT ? RgM_solve_fast(a, b): RgM_RgC_solve_fast(a, b);
2718 415168 : if (!u) return gc_NULL(av);
2719 415070 : if (u != gen_0) return u;
2720 401259 : return RgM_solve_basecase(a, b);
2721 : }
2722 :
2723 : GEN
2724 197331 : RgM_inv(GEN a)
2725 : {
2726 197331 : GEN b = RgM_inv_fast(a);
2727 197318 : return b==gen_0? RgM_solve_basecase(a, NULL): b;
2728 : }
2729 :
2730 : /* assume dim A >= 1, A invertible + upper triangular */
2731 : static GEN
2732 830366 : RgM_inv_upper_ind(GEN A, long index)
2733 : {
2734 830366 : long n = lg(A)-1, i = index, j;
2735 830366 : GEN u = zerocol(n);
2736 830370 : gel(u,i) = ginv(gcoeff(A,i,i));
2737 2377062 : for (i--; i>0; i--)
2738 : {
2739 1546691 : pari_sp av = avma;
2740 1546691 : GEN m = gneg(gmul(gcoeff(A,i,i+1),gel(u,i+1))); /* j = i+1 */
2741 6532125 : for (j=i+2; j<=n; j++) m = gsub(m, gmul(gcoeff(A,i,j),gel(u,j)));
2742 1546684 : gel(u,i) = gerepileupto(av, gdiv(m, gcoeff(A,i,i)));
2743 : }
2744 830371 : return u;
2745 : }
2746 : GEN
2747 277659 : RgM_inv_upper(GEN A)
2748 : {
2749 : long i, l;
2750 277659 : GEN B = cgetg_copy(A, &l);
2751 1108021 : for (i = 1; i < l; i++) gel(B,i) = RgM_inv_upper_ind(A, i);
2752 277663 : return B;
2753 : }
2754 :
2755 : static GEN
2756 1734892 : split_realimag_col(GEN z, long r1, long r2)
2757 : {
2758 1734892 : long i, ru = r1+r2;
2759 1734892 : GEN x = cgetg(ru+r2+1,t_COL), y = x + r2;
2760 4652678 : for (i=1; i<=r1; i++) {
2761 2917783 : GEN a = gel(z,i);
2762 2917783 : if (typ(a) == t_COMPLEX) a = gel(a,1); /* paranoia: a should be real */
2763 2917783 : gel(x,i) = a;
2764 : }
2765 3396197 : for ( ; i<=ru; i++) {
2766 1661302 : GEN b, a = gel(z,i);
2767 1661302 : if (typ(a) == t_COMPLEX) { b = gel(a,2); a = gel(a,1); } else b = gen_0;
2768 1661302 : gel(x,i) = a;
2769 1661302 : gel(y,i) = b;
2770 : }
2771 1734895 : return x;
2772 : }
2773 : GEN
2774 938536 : split_realimag(GEN x, long r1, long r2)
2775 : {
2776 : long i,l; GEN y;
2777 938536 : if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
2778 461680 : y = cgetg_copy(x, &l);
2779 1719721 : for (i=1; i<l; i++) gel(y,i) = split_realimag_col(gel(x,i), r1, r2);
2780 461687 : return y;
2781 : }
2782 :
2783 : /* assume M = (r1+r2) x (r1+2r2) matrix and y compatible vector or matrix
2784 : * r1 first lines of M,y are real. Solve the system obtained by splitting
2785 : * real and imaginary parts. */
2786 : GEN
2787 400133 : RgM_solve_realimag(GEN M, GEN y)
2788 : {
2789 400133 : long l = lg(M), r2 = l - lgcols(M), r1 = l-1 - 2*r2;
2790 400133 : return RgM_solve(split_realimag(M, r1,r2),
2791 : split_realimag(y, r1,r2));
2792 : }
2793 :
2794 : GEN
2795 427 : gauss(GEN a, GEN b)
2796 : {
2797 : GEN z;
2798 427 : long t = typ(b);
2799 427 : if (typ(a)!=t_MAT) pari_err_TYPE("gauss",a);
2800 427 : if (t!=t_COL && t!=t_MAT) pari_err_TYPE("gauss",b);
2801 427 : z = RgM_solve(a,b);
2802 427 : if (!z) pari_err_INV("gauss",a);
2803 322 : return z;
2804 : }
2805 :
2806 : static GEN
2807 735172 : ZlM_gauss_ratlift(GEN a, GEN b, ulong p, long e, GEN C)
2808 : {
2809 735172 : pari_sp av = avma, av2;
2810 : GEN bb, xi, xb, pi, q, B, r;
2811 : long i, f, k;
2812 : ulong mask;
2813 735172 : if (!C) {
2814 0 : C = Flm_inv(ZM_to_Flm(a, p), p);
2815 0 : if (!C) pari_err_INV("ZlM_gauss", a);
2816 : }
2817 735172 : k = f = ZM_max_lg(a)-1;
2818 735176 : mask = quadratic_prec_mask((e+f-1)/f);
2819 735172 : pi = q = powuu(p, f);
2820 735147 : bb = b;
2821 735147 : C = ZpM_invlift(FpM_red(a, q), Flm_to_ZM(C), utoipos(p), f);
2822 735164 : av2 = avma;
2823 735164 : xb = xi = FpM_mul(C, b, q);
2824 770750 : for (i = f; i <= e; i+=f)
2825 : {
2826 136253 : if (i==k)
2827 : {
2828 131127 : k = (mask&1UL) ? 2*k-f: 2*k;
2829 131127 : mask >>= 1;
2830 131127 : B = sqrti(shifti(pi,-1));
2831 131127 : r = FpM_ratlift(xb, pi, B, B, NULL);
2832 131127 : if (r)
2833 : {
2834 108412 : GEN dr, nr = Q_remove_denom(r,&dr);
2835 108412 : if (ZM_equal(ZM_mul(a,nr), dr? ZM_Z_mul(b,dr): b))
2836 : {
2837 100670 : if (DEBUGLEVEL>=4)
2838 0 : err_printf("ZlM_gauss: early solution: %ld/%ld\n",i,e);
2839 100670 : return gerepilecopy(av, r);
2840 : }
2841 : }
2842 : }
2843 35583 : bb = ZM_Z_divexact(ZM_sub(bb, ZM_mul(a, xi)), q);
2844 35583 : if (gc_needed(av,2))
2845 : {
2846 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
2847 0 : gerepileall(av2,3, &pi,&bb,&xb);
2848 : }
2849 35583 : xi = FpM_mul(C, bb, q);
2850 35583 : xb = ZM_add(xb, ZM_Z_mul(xi, pi));
2851 35583 : pi = mulii(pi, q);
2852 : }
2853 634497 : B = sqrti(shifti(pi,-1));
2854 634495 : return gerepileupto(av, FpM_ratlift(xb, pi, B, B, NULL));
2855 : }
2856 :
2857 : /* Dixon p-adic lifting algorithm.
2858 : * Numer. Math. 40, 137-141 (1982), DOI: 10.1007/BF01459082 */
2859 : GEN
2860 737424 : ZM_gauss(GEN a, GEN b0)
2861 : {
2862 737424 : pari_sp av = avma, av2;
2863 : int iscol;
2864 : long n, ncol, i, m, elim;
2865 : ulong p;
2866 737424 : GEN C, delta, nb, nmin, res, b = b0;
2867 : forprime_t S;
2868 :
2869 737424 : if (!init_gauss(a, &b, &n, &ncol, &iscol)) return cgetg(1, iscol?t_COL:t_MAT);
2870 737353 : nb = gen_0; ncol = lg(b);
2871 1637806 : for (i = 1; i < ncol; i++)
2872 : {
2873 900474 : GEN ni = gnorml2(gel(b, i));
2874 900454 : if (cmpii(nb, ni) < 0) nb = ni;
2875 : }
2876 737332 : if (!signe(nb)) {set_avma(av); return iscol? zerocol(n): zeromat(n,lg(b)-1);}
2877 735169 : delta = gen_1; nmin = nb;
2878 3077889 : for (i = 1; i <= n; i++)
2879 : {
2880 2342793 : GEN ni = gnorml2(gel(a, i));
2881 2342896 : if (cmpii(ni, nmin) < 0)
2882 : {
2883 81826 : delta = mulii(delta, nmin); nmin = ni;
2884 : }
2885 : else
2886 2261064 : delta = mulii(delta, ni);
2887 : }
2888 735096 : if (!signe(nmin)) return NULL;
2889 735082 : elim = expi(delta)+1;
2890 735158 : av2 = avma;
2891 735158 : init_modular_big(&S);
2892 : for(;;)
2893 : {
2894 735169 : p = u_forprime_next(&S);
2895 735169 : C = Flm_inv_sp(ZM_to_Flm(a, p), NULL, p);
2896 735182 : if (C) break;
2897 7 : elim -= expu(p);
2898 7 : if (elim < 0) return NULL;
2899 0 : set_avma(av2);
2900 : }
2901 : /* N.B. Our delta/lambda are SQUARES of those in the paper
2902 : * log(delta lambda) / log p, where lambda is 3+sqrt(5) / 2,
2903 : * whose log is < 1, hence + 1 (to cater for rounding errors) */
2904 735175 : m = (long)ceil((dbllog2(delta)*M_LN2 + 1) / log((double)p));
2905 735172 : res = ZlM_gauss_ratlift(a, b, p, m, C);
2906 735169 : if (iscol) return gerepilecopy(av, gel(res, 1));
2907 82048 : return gerepileupto(av, res);
2908 : }
2909 :
2910 : /* #C = n, C[z[i]] = K[i], complete by 0s */
2911 : static GEN
2912 14 : RgC_inflate(GEN K, GEN z, long n)
2913 : {
2914 14 : GEN c = zerocol(n);
2915 14 : long j, l = lg(K);
2916 28 : for (j = 1; j < l; j++) gel(c, z[j]) = gel(K, j);
2917 14 : return c;
2918 : }
2919 : /* in place: C[i] *= cB / v[i] */
2920 : static void
2921 6048 : QC_normalize(GEN C, GEN v, GEN cB)
2922 : {
2923 6048 : long l = lg(C), i;
2924 45465 : for (i = 1; i < l; i++)
2925 : {
2926 39417 : GEN c = cB, k = gel(C,i), d = gel(v,i);
2927 39417 : if (d)
2928 : {
2929 23772 : if (isintzero(d)) { gel(C,i) = gen_0; continue; }
2930 23772 : c = div_content(c, d);
2931 : }
2932 39417 : gel(C,i) = c? gmul(k,c): k;
2933 : }
2934 6048 : }
2935 :
2936 : /* same as above, M rational; if flag = 1, call indexrank and return 1 sol */
2937 : GEN
2938 6041 : QM_gauss_i(GEN M, GEN B, long flag)
2939 : {
2940 6041 : pari_sp av = avma;
2941 : long i, l, n;
2942 6041 : int col = typ(B) == t_COL;
2943 6041 : GEN K, cB, N = cgetg_copy(M, &l), v = cgetg(l, t_VEC), z2 = NULL;
2944 :
2945 45485 : for (i = 1; i < l; i++)
2946 39444 : gel(N,i) = Q_primitive_part(gel(M,i), &gel(v,i));
2947 6041 : if (flag)
2948 : {
2949 301 : GEN z = ZM_indexrank(N), z1 = gel(z,1);
2950 301 : z2 = gel(z,2);
2951 301 : N = shallowmatextract(N, z1, z2);
2952 301 : B = col? vecpermute(B,z1): rowpermute(B,z1);
2953 301 : if (lg(z2) == l) z2 = NULL; else v = vecpermute(v, z2);
2954 : }
2955 6041 : B = Q_primitive_part(B, &cB);
2956 6041 : K = ZM_gauss(N, B); if (!K) return gc_NULL(av);
2957 6041 : n = l - 1;
2958 6041 : if (col)
2959 : {
2960 6027 : QC_normalize(K, v, cB);
2961 6027 : if (z2) K = RgC_inflate(K, z2, n);
2962 : }
2963 : else
2964 : {
2965 14 : long lK = lg(K);
2966 35 : for (i = 1; i < lK; i++)
2967 : {
2968 21 : QC_normalize(gel(K,i), v, cB);
2969 21 : if (z2) gel(K,i) = RgC_inflate(gel(K,i), z2, n);
2970 : }
2971 : }
2972 6041 : return gerepilecopy(av, K);
2973 : }
2974 : GEN
2975 5740 : QM_gauss(GEN M, GEN B) { return QM_gauss_i(M, B, 0); }
2976 :
2977 : static GEN
2978 397717 : ZM_inv_slice(GEN A, GEN P, GEN *mod)
2979 : {
2980 397717 : pari_sp av = avma;
2981 397717 : long i, n = lg(P)-1;
2982 : GEN H, T;
2983 397717 : if (n == 1)
2984 : {
2985 382998 : ulong p = uel(P,1);
2986 382998 : GEN Hp, a = ZM_to_Flm(A, p);
2987 383001 : Hp = Flm_adjoint(a, p);
2988 382999 : Hp = gerepileupto(av, Flm_to_ZM(Hp));
2989 383000 : *mod = utoipos(p); return Hp;
2990 : }
2991 14719 : T = ZV_producttree(P);
2992 14719 : A = ZM_nv_mod_tree(A, P, T);
2993 14719 : H = cgetg(n+1, t_VEC);
2994 67764 : for(i=1; i <= n; i++)
2995 53045 : gel(H,i) = Flm_adjoint(gel(A, i), uel(P,i));
2996 14719 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
2997 14719 : *mod = gmael(T, lg(T)-1, 1);
2998 14719 : gerepileall(av, 2, &H, mod);
2999 14719 : return H;
3000 : }
3001 :
3002 : static GEN
3003 357981 : RgM_true_Hadamard(GEN a)
3004 : {
3005 357981 : pari_sp av = avma;
3006 357981 : long n = lg(a)-1, i;
3007 : GEN B;
3008 357981 : if (n == 0) return gen_1;
3009 357981 : a = RgM_gtofp(a, LOWDEFAULTPREC);
3010 357980 : B = gnorml2(gel(a,1));
3011 1732925 : for (i = 2; i <= n; i++) B = gmul(B, gnorml2(gel(a,i)));
3012 357976 : return gerepileuptoint(av, ceil_safe(sqrtr(B)));
3013 : }
3014 :
3015 : GEN
3016 397718 : ZM_inv_worker(GEN P, GEN A)
3017 : {
3018 397718 : GEN V = cgetg(3, t_VEC);
3019 397718 : gel(V,1) = ZM_inv_slice(A, P, &gel(V,2));
3020 397715 : return V;
3021 : }
3022 :
3023 : static GEN
3024 40655 : ZM_inv0(GEN A, GEN *pden)
3025 : {
3026 40655 : if (pden) *pden = gen_1;
3027 40655 : (void)A; return cgetg(1, t_MAT);
3028 : }
3029 : static GEN
3030 65942 : ZM_inv1(GEN A, GEN *pden)
3031 : {
3032 65942 : GEN a = gcoeff(A,1,1);
3033 65942 : long s = signe(a);
3034 65942 : if (!s) return NULL;
3035 65942 : if (pden) *pden = absi(a);
3036 65942 : retmkmat(mkcol(s == 1? gen_1: gen_m1));
3037 : }
3038 : static GEN
3039 118239 : ZM_inv2(GEN A, GEN *pden)
3040 : {
3041 : GEN a, b, c, d, D, cA;
3042 : long s;
3043 118239 : A = Q_primitive_part(A, &cA);
3044 118235 : a = gcoeff(A,1,1); b = gcoeff(A,1,2);
3045 118235 : c = gcoeff(A,2,1); d = gcoeff(A,2,2);
3046 118235 : D = subii(mulii(a,d), mulii(b,c)); /* left on stack */
3047 118225 : s = signe(D);
3048 118225 : if (!s) return NULL;
3049 118225 : if (s < 0) D = negi(D);
3050 118231 : if (pden) *pden = mul_denom(D, cA);
3051 118232 : if (s > 0)
3052 78426 : retmkmat2(mkcol2(icopy(d), negi(c)), mkcol2(negi(b), icopy(a)));
3053 : else
3054 39806 : retmkmat2(mkcol2(negi(d), icopy(c)), mkcol2(icopy(b), negi(a)));
3055 : }
3056 :
3057 : /* to be used when denom(M^(-1)) << det(M) and a sharp multiple is
3058 : * not available. Return H primitive such that M*H = den*Id */
3059 : GEN
3060 0 : ZM_inv_ratlift(GEN M, GEN *pden)
3061 : {
3062 0 : pari_sp av2, av = avma;
3063 : GEN Hp, q, H;
3064 : ulong p;
3065 0 : long m = lg(M)-1;
3066 : forprime_t S;
3067 : pari_timer ti;
3068 :
3069 0 : if (m == 0) return ZM_inv0(M,pden);
3070 0 : if (m == 1 && nbrows(M)==1) return ZM_inv1(M,pden);
3071 0 : if (m == 2 && nbrows(M)==2) return ZM_inv2(M,pden);
3072 :
3073 0 : if (DEBUGLEVEL>5) timer_start(&ti);
3074 0 : init_modular_big(&S);
3075 0 : av2 = avma;
3076 0 : H = NULL;
3077 0 : while ((p = u_forprime_next(&S)))
3078 : {
3079 : GEN Mp, B, Hr;
3080 0 : Mp = ZM_to_Flm(M,p);
3081 0 : Hp = Flm_inv_sp(Mp, NULL, p);
3082 0 : if (!Hp) continue;
3083 0 : if (!H)
3084 : {
3085 0 : H = ZM_init_CRT(Hp, p);
3086 0 : q = utoipos(p);
3087 : }
3088 : else
3089 0 : ZM_incremental_CRT(&H, Hp, &q, p);
3090 0 : B = sqrti(shifti(q,-1));
3091 0 : Hr = FpM_ratlift(H,q,B,B,NULL);
3092 0 : if (DEBUGLEVEL>5)
3093 0 : timer_printf(&ti,"ZM_inv mod %lu (ratlift=%ld)", p,!!Hr);
3094 0 : if (Hr) {/* DONE ? */
3095 0 : GEN Hl = Q_remove_denom(Hr, pden);
3096 0 : if (ZM_isscalar(ZM_mul(Hl, M), *pden)) { H = Hl; break; }
3097 : }
3098 :
3099 0 : if (gc_needed(av,2))
3100 : {
3101 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZM_inv_ratlift");
3102 0 : gerepileall(av2, 2, &H, &q);
3103 : }
3104 : }
3105 0 : if (!*pden) *pden = gen_1;
3106 0 : gerepileall(av, 2, &H, pden);
3107 0 : return H;
3108 : }
3109 :
3110 : GEN
3111 74105 : FpM_ratlift_worker(GEN A, GEN mod, GEN B)
3112 : {
3113 : long l, i;
3114 74105 : GEN H = cgetg_copy(A, &l);
3115 150158 : for (i = 1; i < l; i++)
3116 : {
3117 76336 : GEN c = FpC_ratlift(gel(A,i), mod, B, B, NULL);
3118 76087 : gel(H,i) = c? c: gen_0;
3119 : }
3120 73822 : return H;
3121 : }
3122 : static int
3123 385094 : can_ratlift(GEN x, GEN mod, GEN B)
3124 : {
3125 385094 : pari_sp av = avma;
3126 : GEN a, b;
3127 385094 : return gc_bool(av, Fp_ratlift(x, mod, B, B, &a,&b));
3128 : }
3129 : static GEN
3130 445891 : FpM_ratlift_parallel(GEN A, GEN mod, GEN B)
3131 : {
3132 445891 : pari_sp av = avma;
3133 : GEN worker;
3134 445891 : long i, l = lg(A), m = mt_nbthreads();
3135 445893 : int test = !!B;
3136 :
3137 445893 : if (l == 1 || lgcols(A) == 1) return gcopy(A);
3138 445897 : if (!B) B = sqrti(shifti(mod,-1));
3139 445898 : if (m == 1 || l == 2 || lgcols(A) < 10)
3140 : {
3141 438972 : A = FpM_ratlift(A, mod, B, B, NULL);
3142 438972 : return A? A: gc_NULL(av);
3143 : }
3144 : /* test one coefficient first */
3145 6926 : if (test && !can_ratlift(gcoeff(A,1,1), mod, B)) return gc_NULL(av);
3146 6804 : worker = snm_closure(is_entry("_FpM_ratlift_worker"), mkvec2(mod,B));
3147 6804 : A = gen_parapply_slice(worker, A, m);
3148 75596 : for (i = 1; i < l; i++) if (typ(gel(A,i)) != t_COL) return gc_NULL(av);
3149 5898 : return A;
3150 : }
3151 :
3152 : static GEN
3153 378361 : ZM_adj_ratlift(GEN A, GEN H, GEN mod, GEN T)
3154 : {
3155 378361 : pari_sp av = avma;
3156 : GEN B, D, g;
3157 378361 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3158 378358 : if (T) D = mulii(T, D);
3159 378358 : g = gcdii(D, mod);
3160 378354 : if (!equali1(g))
3161 : {
3162 14 : mod = diviiexact(mod, g);
3163 14 : H = FpM_red(H, mod);
3164 : }
3165 378352 : D = Fp_inv(Fp_red(D, mod), mod);
3166 : /* test 1 coeff first */
3167 378353 : B = sqrti(shifti(mod,-1));
3168 378350 : if (!can_ratlift(Fp_mul(D, gcoeff(A,1,1), mod), mod, B)) return gc_NULL(av);
3169 365355 : H = FpM_Fp_mul(H, D, mod);
3170 365349 : H = FpM_ratlift_parallel(H, mod, B);
3171 365358 : return H? H: gc_NULL(av);
3172 : }
3173 :
3174 : /* if (T) return T A^(-1) in Mn(Q), else B in Mn(Z) such that A B = den*Id */
3175 : static GEN
3176 582822 : ZM_inv_i(GEN A, GEN *pden, GEN T)
3177 : {
3178 582822 : pari_sp av = avma;
3179 582822 : long m = lg(A)-1, n, k1 = 1, k2;
3180 582822 : GEN H = NULL, D, H1 = NULL, mod1 = NULL, worker;
3181 : ulong bnd, mask;
3182 : forprime_t S;
3183 : pari_timer ti;
3184 :
3185 582822 : if (m == 0) return ZM_inv0(A,pden);
3186 542167 : if (pden) *pden = gen_1;
3187 542167 : if (nbrows(A) < m) return NULL;
3188 542162 : if (m == 1 && nbrows(A)==1 && !T) return ZM_inv1(A,pden);
3189 476220 : if (m == 2 && nbrows(A)==2 && !T) return ZM_inv2(A,pden);
3190 :
3191 357980 : if (DEBUGLEVEL>=5) timer_start(&ti);
3192 357980 : init_modular_big(&S);
3193 357981 : bnd = expi(RgM_true_Hadamard(A));
3194 357979 : worker = snm_closure(is_entry("_ZM_inv_worker"), mkvec(A));
3195 357984 : gen_inccrt("ZM_inv_r", worker, NULL, k1, 0, &S, &H1, &mod1, nmV_chinese_center, FpM_center);
3196 357983 : n = (bnd+1)/expu(S.p)+1;
3197 357983 : if (DEBUGLEVEL>=5) timer_printf(&ti,"inv (%ld/%ld primes)", k1, n);
3198 357983 : mask = quadratic_prec_mask(n);
3199 357983 : for (k2 = 0;;)
3200 39465 : {
3201 : GEN Hr;
3202 397448 : if (k2 > 0)
3203 : {
3204 33984 : gen_inccrt("ZM_inv_r", worker, NULL, k2, 0, &S, &H1, &mod1,nmV_chinese_center,FpM_center);
3205 33984 : k1 += k2;
3206 33984 : if (DEBUGLEVEL>=5) timer_printf(&ti,"CRT (%ld/%ld primes)", k1, n);
3207 : }
3208 397448 : if (mask == 1) break;
3209 378361 : k2 = (mask&1UL) ? k1-1: k1;
3210 378361 : mask >>= 1;
3211 :
3212 378361 : Hr = ZM_adj_ratlift(A, H1, mod1, T);
3213 378358 : if (DEBUGLEVEL>=5) timer_printf(&ti,"ratlift (%ld/%ld primes)", k1, n);
3214 378358 : if (Hr) {/* DONE ? */
3215 341219 : GEN Hl = Q_primpart(Hr), R = ZM_mul(Hl, A), d = gcoeff(R,1,1);
3216 341217 : if (gsigne(d) < 0) { d = gneg(d); Hl = ZM_neg(Hl); }
3217 341216 : if (DEBUGLEVEL>=5) timer_printf(&ti,"mult (%ld/%ld primes)", k1, n);
3218 341216 : if (equali1(d))
3219 : {
3220 252199 : if (ZM_isidentity(R)) { H = Hl; break; }
3221 : }
3222 89019 : else if (ZM_isscalar(R, d))
3223 : {
3224 86695 : if (T) T = gdiv(T,d);
3225 83473 : else if (pden) *pden = d;
3226 86695 : H = Hl; break;
3227 : }
3228 : }
3229 : }
3230 357981 : if (!H)
3231 : {
3232 : GEN d;
3233 19087 : H = H1;
3234 19087 : D = ZMrow_ZC_mul(H, gel(A,1), 1);
3235 19087 : if (signe(D)==0) pari_err_INV("ZM_inv", A);
3236 19087 : if (T) T = gdiv(T, D);
3237 : else
3238 : {
3239 18165 : d = gcdii(Q_content_safe(H), D);
3240 18165 : if (signe(D) < 0) d = negi(d);
3241 18165 : if (!equali1(d))
3242 : {
3243 11235 : H = ZM_Z_divexact(H, d);
3244 11235 : D = diviiexact(D, d);
3245 : }
3246 18164 : if (pden) *pden = D;
3247 : }
3248 : }
3249 357980 : if (T && !isint1(T)) H = ZM_Q_mul(H, T);
3250 357980 : gerepileall(av, pden? 2: 1, &H, pden);
3251 357984 : return H;
3252 : }
3253 : GEN
3254 523431 : ZM_inv(GEN A, GEN *pden) { return ZM_inv_i(A, pden, NULL); }
3255 :
3256 : /* same as above, M rational */
3257 : GEN
3258 59388 : QM_inv(GEN M)
3259 : {
3260 59388 : pari_sp av = avma;
3261 : GEN den, dM, K;
3262 59388 : M = Q_remove_denom(M, &dM);
3263 59388 : K = ZM_inv_i(M, &den, dM);
3264 59387 : if (!K) return gc_NULL(av);
3265 59380 : if (den && !equali1(den)) K = ZM_Q_mul(K, ginv(den));
3266 59367 : return gerepileupto(av, K);
3267 : }
3268 :
3269 : static GEN
3270 105201 : ZM_ker_filter(GEN A, GEN P)
3271 : {
3272 105201 : long i, j, l = lg(A), n = 1, d = lg(gmael(A,1,1));
3273 105201 : GEN B, Q, D = gmael(A,1,2);
3274 215143 : for (i=2; i<l; i++)
3275 : {
3276 109942 : GEN Di = gmael(A,i,2);
3277 109942 : long di = lg(gmael(A,i,1));
3278 109942 : int c = vecsmall_lexcmp(D, Di);
3279 109942 : if (di==d && c==0) n++;
3280 45588 : else if (d > di || (di==d && c>0))
3281 37680 : { n = 1; d = di; D = Di; }
3282 : }
3283 105201 : B = cgetg(n+1, t_VEC);
3284 105202 : Q = cgetg(n+1, typ(P));
3285 320340 : for (i=1, j=1; i<l; i++)
3286 : {
3287 215140 : if (lg(gmael(A,i,1))==d && vecsmall_lexcmp(D, gmael(A,i,2))==0)
3288 : {
3289 169553 : gel(B,j) = gmael(A,i,1);
3290 169553 : Q[j] = P[i];
3291 169553 : j++;
3292 : }
3293 : }
3294 105200 : return mkvec3(B,Q,D);
3295 : }
3296 :
3297 : static GEN
3298 69539 : ZM_ker_chinese(GEN A, GEN P, GEN *mod)
3299 : {
3300 69539 : GEN BQD = ZM_ker_filter(A, P);
3301 69539 : return mkvec2(nmV_chinese_center(gel(BQD,1), gel(BQD,2), mod), gel(BQD,3));
3302 : }
3303 :
3304 : static GEN
3305 129250 : ZM_ker_slice(GEN A, GEN P, GEN *mod)
3306 : {
3307 129250 : pari_sp av = avma;
3308 129250 : long i, n = lg(P)-1;
3309 : GEN BQD, D, H, T, Q;
3310 129250 : if (n == 1)
3311 : {
3312 93585 : ulong p = uel(P,1);
3313 93585 : GEN K = Flm_ker_sp(ZM_to_Flm(A, p), p, 2);
3314 93585 : *mod = utoipos(p); return mkvec2(Flm_to_ZM(gel(K,1)), gel(K,2));
3315 : }
3316 35665 : T = ZV_producttree(P);
3317 35665 : A = ZM_nv_mod_tree(A, P, T);
3318 35664 : H = cgetg(n+1, t_VEC);
3319 111490 : for(i=1 ; i <= n; i++)
3320 75829 : gel(H,i) = Flm_ker_sp(gel(A, i), P[i], 2);
3321 35661 : BQD = ZM_ker_filter(H, P); Q = gel(BQD,2);
3322 35661 : if (lg(Q) != lg(P)) T = ZV_producttree(Q);
3323 35661 : H = nmV_chinese_center_tree_seq(gel(BQD,1), Q, T, ZV_chinesetree(Q,T));
3324 35664 : *mod = gmael(T, lg(T)-1, 1);
3325 35664 : D = gel(BQD, 3);
3326 35664 : gerepileall(av, 3, &H, &D, mod);
3327 35665 : return mkvec2(H,D);
3328 : }
3329 :
3330 : GEN
3331 129250 : ZM_ker_worker(GEN P, GEN A)
3332 : {
3333 129250 : GEN V = cgetg(3, t_VEC);
3334 129250 : gel(V,1) = ZM_ker_slice(A, P, &gel(V,2));
3335 129250 : return V;
3336 : }
3337 :
3338 : /* assume lg(A) > 1 */
3339 : static GEN
3340 62507 : ZM_ker_i(GEN A)
3341 : {
3342 : pari_sp av;
3343 62507 : long k, m = lg(A)-1;
3344 62507 : GEN HD = NULL, mod = gen_1, worker;
3345 : forprime_t S;
3346 :
3347 62507 : if (m >= 2*nbrows(A))
3348 : {
3349 3031 : GEN v = ZM_indexrank(A), y = gel(v,2), z = indexcompl(y, m);
3350 : GEN B, A1, A1i, d;
3351 3031 : A = rowpermute(A, gel(v,1)); /* same kernel */
3352 3031 : A1 = vecpermute(A, y); /* maximal rank submatrix */
3353 3031 : B = vecpermute(A, z);
3354 3031 : A1i = ZM_inv(A1, &d);
3355 3031 : if (!d) d = gen_1;
3356 3031 : B = vconcat(ZM_mul(ZM_neg(A1i), B), scalarmat_shallow(d, lg(B)-1));
3357 3031 : if (!gequal(y, identity_perm(lg(y)-1)))
3358 665 : B = rowpermute(B, perm_inv(shallowconcat(y,z)));
3359 3031 : return vec_Q_primpart(B);
3360 : }
3361 59476 : init_modular_big(&S);
3362 59476 : worker = snm_closure(is_entry("_ZM_ker_worker"), mkvec(A));
3363 59476 : av = avma;
3364 59476 : for (k = 1;; k <<= 1)
3365 65326 : {
3366 : pari_timer ti;
3367 : GEN H, Hr;
3368 124802 : gen_inccrt_i("ZM_ker", worker, NULL, (k+1)>>1, 0,
3369 : &S, &HD, &mod, ZM_ker_chinese, NULL);
3370 124802 : gerepileall(av, 2, &HD, &mod);
3371 140017 : H = gel(HD, 1); if (lg(H) == 1) return H;
3372 80541 : if (DEBUGLEVEL >= 4) timer_start(&ti);
3373 80541 : Hr = FpM_ratlift_parallel(H, mod, NULL);
3374 80541 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: ratlift (%ld)",!!Hr);
3375 80541 : if (Hr)
3376 : {
3377 : GEN MH;
3378 69811 : Hr = vec_Q_primpart(Hr);
3379 69811 : MH = ZM_mul(A, Hr);
3380 69811 : if (DEBUGLEVEL >= 4) timer_printf(&ti,"ZM_ker: QM_mul");
3381 69811 : if (ZM_equal0(MH)) return Hr;
3382 : }
3383 : }
3384 : }
3385 :
3386 : GEN
3387 47016 : ZM_ker(GEN M)
3388 : {
3389 47016 : pari_sp av = avma;
3390 47016 : long l = lg(M)-1;
3391 47016 : if (l==0) return cgetg(1, t_MAT);
3392 47016 : if (lgcols(M)==1) return matid(l);
3393 47016 : return gerepilecopy(av, ZM_ker_i(M));
3394 : }
3395 :
3396 : GEN
3397 16331 : QM_ker(GEN M)
3398 : {
3399 16331 : pari_sp av = avma;
3400 16331 : long l = lg(M)-1;
3401 16331 : if (l==0) return cgetg(1, t_MAT);
3402 16289 : if (lgcols(M)==1) return matid(l);
3403 15428 : return gerepilecopy(av, ZM_ker_i(row_Q_primpart(M)));
3404 : }
3405 :
3406 : /* x a ZM. Return a multiple of the determinant of the lattice generated by
3407 : * the columns of x. From Algorithm 2.2.6 in GTM138 */
3408 : GEN
3409 47549 : detint(GEN A)
3410 : {
3411 47549 : if (typ(A) != t_MAT) pari_err_TYPE("detint",A);
3412 47549 : RgM_check_ZM(A, "detint");
3413 47549 : return ZM_detmult(A);
3414 : }
3415 : GEN
3416 154655 : ZM_detmult(GEN A)
3417 : {
3418 154655 : pari_sp av1, av = avma;
3419 : GEN B, c, v, piv;
3420 154655 : long rg, i, j, k, m, n = lg(A) - 1;
3421 :
3422 154655 : if (!n) return gen_1;
3423 154655 : m = nbrows(A);
3424 154655 : if (n < m) return gen_0;
3425 154585 : c = zero_zv(m);
3426 154585 : av1 = avma;
3427 154585 : B = zeromatcopy(m,m);
3428 154585 : v = cgetg(m+1, t_COL);
3429 154585 : piv = gen_1; rg = 0;
3430 679438 : for (k=1; k<=n; k++)
3431 : {
3432 679424 : GEN pivprec = piv;
3433 679424 : long t = 0;
3434 4901587 : for (i=1; i<=m; i++)
3435 : {
3436 4222168 : pari_sp av2 = avma;
3437 : GEN vi;
3438 4222168 : if (c[i]) continue;
3439 :
3440 2451042 : vi = mulii(piv, gcoeff(A,i,k));
3441 22439561 : for (j=1; j<=m; j++)
3442 19988478 : if (c[j]) vi = addii(vi, mulii(gcoeff(B,j,i),gcoeff(A,j,k)));
3443 2451083 : if (!t && signe(vi)) t = i;
3444 2451083 : gel(v,i) = gerepileuptoint(av2, vi);
3445 : }
3446 679419 : if (!t) continue;
3447 : /* at this point c[t] = 0 */
3448 :
3449 679335 : if (++rg >= m) { /* full rank; mostly done */
3450 154571 : GEN det = gel(v,t); /* last on stack */
3451 154571 : if (++k > n)
3452 154446 : det = absi(det);
3453 : else
3454 : {
3455 : /* improve further; at this point c[i] is set for all i != t */
3456 125 : gcoeff(B,t,t) = piv; v = centermod(gel(B,t), det);
3457 418 : for ( ; k<=n; k++)
3458 293 : det = gcdii(det, ZV_dotproduct(v, gel(A,k)));
3459 : }
3460 154571 : return gerepileuptoint(av, det);
3461 : }
3462 :
3463 524764 : piv = gel(v,t);
3464 4067107 : for (i=1; i<=m; i++)
3465 : {
3466 : GEN mvi;
3467 3542343 : if (c[i] || i == t) continue;
3468 :
3469 1771170 : gcoeff(B,t,i) = mvi = negi(gel(v,i));
3470 17534203 : for (j=1; j<=m; j++)
3471 15763033 : if (c[j]) /* implies j != t */
3472 : {
3473 4073565 : pari_sp av2 = avma;
3474 4073565 : GEN z = addii(mulii(gcoeff(B,j,i), piv), mulii(gcoeff(B,j,t), mvi));
3475 4073563 : if (rg > 1) z = diviiexact(z, pivprec);
3476 4073563 : gcoeff(B,j,i) = gerepileuptoint(av2, z);
3477 : }
3478 : }
3479 524764 : c[t] = k;
3480 524764 : if (gc_needed(av,1))
3481 : {
3482 0 : if(DEBUGMEM>1) pari_warn(warnmem,"detint. k=%ld",k);
3483 0 : gerepileall(av1, 2, &piv,&B); v = zerovec(m);
3484 : }
3485 : }
3486 14 : return gc_const(av, gen_0);
3487 : }
3488 :
3489 : /* Reduce x modulo (invertible) y */
3490 : GEN
3491 13433 : closemodinvertible(GEN x, GEN y)
3492 : {
3493 13433 : return gmul(y, ground(RgM_solve(y,x)));
3494 : }
3495 : GEN
3496 7 : reducemodinvertible(GEN x, GEN y)
3497 : {
3498 7 : return gsub(x, closemodinvertible(x,y));
3499 : }
3500 : GEN
3501 0 : reducemodlll(GEN x,GEN y)
3502 : {
3503 0 : return reducemodinvertible(x, ZM_lll(y, 0.75, LLL_INPLACE));
3504 : }
3505 :
3506 : /*******************************************************************/
3507 : /* */
3508 : /* KERNEL of an m x n matrix */
3509 : /* return n - rk(x) linearly independent vectors */
3510 : /* */
3511 : /*******************************************************************/
3512 : static GEN
3513 28 : RgM_deplin_i(GEN x0)
3514 : {
3515 28 : pari_sp av = avma, av2;
3516 28 : long i, j, k, nl, nc = lg(x0)-1;
3517 : GEN D, x, y, c, l, d, ck;
3518 :
3519 28 : if (!nc) return NULL;
3520 28 : nl = nbrows(x0);
3521 28 : c = zero_zv(nl);
3522 28 : l = cgetg(nc+1, t_VECSMALL); /* not initialized */
3523 28 : av2 = avma;
3524 28 : x = RgM_shallowcopy(x0);
3525 28 : d = const_vec(nl, gen_1); /* pivot list */
3526 28 : ck = NULL; /* gcc -Wall */
3527 98 : for (k=1; k<=nc; k++)
3528 : {
3529 91 : ck = gel(x,k);
3530 196 : for (j=1; j<k; j++)
3531 : {
3532 105 : GEN cj = gel(x,j), piv = gel(d,j), q = gel(ck,l[j]);
3533 420 : for (i=1; i<=nl; i++)
3534 315 : if (i!=l[j]) gel(ck,i) = gsub(gmul(piv, gel(ck,i)), gmul(q, gel(cj,i)));
3535 : }
3536 :
3537 91 : i = gauss_get_pivot_NZ(x, NULL, k, c);
3538 91 : if (i > nl) break;
3539 70 : if (gc_needed(av,1))
3540 : {
3541 0 : if (DEBUGMEM>1) pari_warn(warnmem,"deplin k = %ld/%ld",k,nc);
3542 0 : gerepileall(av2, 2, &x, &d);
3543 0 : ck = gel(x,k);
3544 : }
3545 70 : gel(d,k) = gel(ck,i);
3546 70 : c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */
3547 : }
3548 28 : if (k > nc) return gc_NULL(av);
3549 21 : if (k == 1) { set_avma(av); return scalarcol_shallow(gen_1,nc); }
3550 21 : y = cgetg(nc+1,t_COL);
3551 21 : gel(y,1) = gcopy(gel(ck, l[1]));
3552 49 : for (D=gel(d,1),j=2; j<k; j++)
3553 : {
3554 28 : gel(y,j) = gmul(gel(ck, l[j]), D);
3555 28 : D = gmul(D, gel(d,j));
3556 : }
3557 21 : gel(y,j) = gneg(D);
3558 21 : for (j++; j<=nc; j++) gel(y,j) = gen_0;
3559 21 : y = primitive_part(y, &c);
3560 21 : return c? gerepileupto(av, y): gerepilecopy(av, y);
3561 : }
3562 : static GEN
3563 0 : RgV_deplin(GEN v)
3564 : {
3565 0 : pari_sp av = avma;
3566 0 : long n = lg(v)-1;
3567 0 : GEN y, p = NULL;
3568 0 : if (n <= 1)
3569 : {
3570 0 : if (n == 1 && gequal0(gel(v,1))) return mkcol(gen_1);
3571 0 : return cgetg(1, t_COL);
3572 : }
3573 0 : if (gequal0(gel(v,1))) return scalarcol_shallow(gen_1, n);
3574 0 : v = primpart(mkvec2(gel(v,1),gel(v,2)));
3575 0 : if (RgV_is_FpV(v, &p) && p) v = centerlift(v);
3576 0 : y = zerocol(n);
3577 0 : gel(y,1) = gneg(gel(v,2));
3578 0 : gel(y,2) = gcopy(gel(v,1));
3579 0 : return gerepileupto(av, y);
3580 :
3581 : }
3582 :
3583 : static GEN
3584 105 : RgM_deplin_FpM(GEN x, GEN p)
3585 : {
3586 105 : pari_sp av = avma;
3587 : ulong pp;
3588 105 : x = RgM_Fp_init(x, p, &pp);
3589 105 : switch(pp)
3590 : {
3591 35 : case 0:
3592 35 : x = FpM_ker_gen(x,p,1);
3593 35 : if (!x) return gc_NULL(av);
3594 21 : x = FpC_center(x,p,shifti(p,-1));
3595 21 : break;
3596 14 : case 2:
3597 14 : x = F2m_ker_sp(x,1);
3598 14 : if (!x) return gc_NULL(av);
3599 7 : x = F2c_to_ZC(x); break;
3600 56 : default:
3601 56 : x = Flm_ker_sp(x,pp,1);
3602 56 : if (!x) return gc_NULL(av);
3603 35 : x = Flv_center(x, pp, pp>>1);
3604 35 : x = zc_to_ZC(x);
3605 35 : break;
3606 : }
3607 63 : return gerepileupto(av, x);
3608 : }
3609 :
3610 : /* FIXME: implement direct modular ZM_deplin ? */
3611 : static GEN
3612 98 : QM_deplin(GEN M)
3613 : {
3614 98 : pari_sp av = avma;
3615 98 : long l = lg(M)-1;
3616 : GEN k;
3617 98 : if (l==0) return NULL;
3618 63 : if (lgcols(M)==1) return col_ei(l, 1);
3619 63 : k = ZM_ker_i(row_Q_primpart(M));
3620 63 : if (lg(k)== 1) return gc_NULL(av);
3621 49 : return gerepilecopy(av, gel(k,1));
3622 : }
3623 :
3624 : static GEN
3625 42 : RgM_deplin_FqM(GEN x, GEN pol, GEN p)
3626 : {
3627 42 : pari_sp av = avma;
3628 42 : GEN b, T = RgX_to_FpX(pol, p);
3629 42 : if (signe(T) == 0) pari_err_OP("deplin",x,pol);
3630 42 : b = FqM_deplin(RgM_to_FqM(x, T, p), T, p);
3631 42 : return gerepileupto(av, b);
3632 : }
3633 :
3634 : #define code(t1,t2) ((t1 << 6) | t2)
3635 : static GEN
3636 357 : RgM_deplin_fast(GEN x)
3637 : {
3638 : GEN p, pol;
3639 : long pa;
3640 357 : long t = RgM_type(x, &p,&pol,&pa);
3641 357 : switch(t)
3642 : {
3643 98 : case t_INT: /* fall through */
3644 98 : case t_FRAC: return QM_deplin(x);
3645 84 : case t_FFELT: return FFM_deplin(x, pol);
3646 105 : case t_INTMOD: return RgM_deplin_FpM(x, p);
3647 42 : case code(t_POLMOD, t_INTMOD):
3648 42 : return RgM_deplin_FqM(x, pol, p);
3649 28 : default: return gen_0;
3650 : }
3651 : }
3652 : #undef code
3653 :
3654 : static GEN
3655 357 : RgM_deplin(GEN x)
3656 : {
3657 357 : GEN z = RgM_deplin_fast(x);
3658 357 : if (z!= gen_0) return z;
3659 28 : return RgM_deplin_i(x);
3660 : }
3661 :
3662 : GEN
3663 357 : deplin(GEN x)
3664 : {
3665 357 : switch(typ(x))
3666 : {
3667 357 : case t_MAT:
3668 : {
3669 357 : GEN z = RgM_deplin(x);
3670 357 : if (z) return z;
3671 140 : return cgetg(1, t_COL);
3672 : }
3673 0 : case t_VEC: return RgV_deplin(x);
3674 0 : default: pari_err_TYPE("deplin",x);
3675 : }
3676 : return NULL;/*LCOV_EXCL_LINE*/
3677 : }
3678 :
3679 : /*******************************************************************/
3680 : /* */
3681 : /* GAUSS REDUCTION OF MATRICES (m lines x n cols) */
3682 : /* (kernel, image, complementary image, rank) */
3683 : /* */
3684 : /*******************************************************************/
3685 : /* return the transform of x under a standard Gauss pivot.
3686 : * x0 is a reference point when guessing whether x[i,j] ~ 0
3687 : * (iff x[i,j] << x0[i,j])
3688 : * Set r = dim ker(x). d[k] contains the index of the first nonzero pivot
3689 : * in column k */
3690 : static GEN
3691 1042 : gauss_pivot_ker(GEN x, GEN x0, GEN *dd, long *rr)
3692 : {
3693 : GEN c, d, p, data;
3694 : pari_sp av;
3695 : long i, j, k, r, t, n, m;
3696 : pivot_fun pivot;
3697 :
3698 1042 : n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
3699 1042 : m=nbrows(x); r=0;
3700 1042 : pivot = get_pivot_fun(x, x0, &data);
3701 1042 : x = RgM_shallowcopy(x);
3702 1042 : c = zero_zv(m);
3703 1042 : d = cgetg(n+1,t_VECSMALL);
3704 1042 : av=avma;
3705 6204 : for (k=1; k<=n; k++)
3706 : {
3707 5162 : j = pivot(x, data, k, c);
3708 5162 : if (j > m)
3709 : {
3710 1099 : r++; d[k]=0;
3711 5026 : for(j=1; j<k; j++)
3712 3927 : if (d[j]) gcoeff(x,d[j],k) = gclone(gcoeff(x,d[j],k));
3713 : }
3714 : else
3715 : { /* pivot for column k on row j */
3716 4063 : c[j]=k; d[k]=j; p = gdiv(gen_m1,gcoeff(x,j,k));
3717 4063 : gcoeff(x,j,k) = gen_m1;
3718 : /* x[j,] /= - x[j,k] */
3719 21837 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3720 38774 : for (t=1; t<=m; t++)
3721 34711 : if (t!=j)
3722 : { /* x[t,] -= 1 / x[j,k] x[j,] */
3723 30648 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3724 30648 : if (gequal0(p)) continue;
3725 81156 : for (i=k+1; i<=n; i++)
3726 65411 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
3727 15745 : if (gc_needed(av,1)) gerepile_gauss_ker(x,k,t,av);
3728 : }
3729 : }
3730 : }
3731 1042 : *dd=d; *rr=r; return x;
3732 : }
3733 :
3734 : /* r = dim ker(x).
3735 : * Returns d:
3736 : * d[k] != 0 contains the index of a nonzero pivot in column k
3737 : * d[k] == 0 if column k is a linear combination of the (k-1) first ones */
3738 : GEN
3739 160955 : RgM_pivots(GEN x0, GEN data, long *rr, pivot_fun pivot)
3740 : {
3741 : GEN x, c, d, p;
3742 160955 : long i, j, k, r, t, m, n = lg(x0)-1;
3743 : pari_sp av;
3744 :
3745 160955 : if (RgM_is_ZM(x0)) return ZM_pivots(x0, rr);
3746 138725 : if (!n) { *rr = 0; return NULL; }
3747 :
3748 138725 : d = cgetg(n+1, t_VECSMALL);
3749 138724 : x = RgM_shallowcopy(x0);
3750 138723 : m = nbrows(x); r = 0;
3751 138723 : c = zero_zv(m);
3752 138789 : av = avma;
3753 3344699 : for (k=1; k<=n; k++)
3754 : {
3755 3205976 : j = pivot(x, data, k, c);
3756 3205992 : if (j > m) { r++; d[k] = 0; }
3757 : else
3758 : {
3759 280105 : c[j] = k; d[k] = j; p = gdiv(gen_m1, gcoeff(x,j,k));
3760 10022490 : for (i=k+1; i<=n; i++) gcoeff(x,j,i) = gmul(p,gcoeff(x,j,i));
3761 :
3762 1134245 : for (t=1; t<=m; t++)
3763 854222 : if (!c[t]) /* no pivot on that line yet */
3764 : {
3765 334649 : p = gcoeff(x,t,k); gcoeff(x,t,k) = gen_0;
3766 24487177 : for (i=k+1; i<=n; i++)
3767 24152533 : gcoeff(x,t,i) = gadd(gcoeff(x,t,i), gmul(p, gcoeff(x,j,i)));
3768 334644 : if (gc_needed(av,1)) gerepile_gauss(x,k,t,av,j,c);
3769 : }
3770 10303056 : for (i=k; i<=n; i++) gcoeff(x,j,i) = gen_0; /* dummy */
3771 : }
3772 : }
3773 138723 : *rr = r; return gc_const((pari_sp)d, d);
3774 : }
3775 :
3776 : static long
3777 369139 : ZM_count_0_cols(GEN M)
3778 : {
3779 369139 : long i, l = lg(M), n = 0;
3780 1910127 : for (i = 1; i < l; i++)
3781 1540986 : if (ZV_equal0(gel(M,i))) n++;
3782 369141 : return n;
3783 : }
3784 :
3785 : static void indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol);
3786 : /* As RgM_pivots, integer entries. Set *rr = dim Ker M0 */
3787 : GEN
3788 395463 : ZM_pivots(GEN M0, long *rr)
3789 : {
3790 395463 : GEN d, dbest = NULL;
3791 : long m, mm, n, nn, i, imax, rmin, rbest, zc;
3792 395463 : int beenthere = 0;
3793 395463 : pari_sp av, av0 = avma;
3794 : forprime_t S;
3795 :
3796 395463 : rbest = n = lg(M0)-1;
3797 395463 : if (n == 0) { *rr = 0; return NULL; }
3798 369139 : zc = ZM_count_0_cols(M0);
3799 369141 : if (n == zc) { *rr = zc; return zero_zv(n); }
3800 :
3801 368679 : m = nbrows(M0);
3802 368676 : rmin = maxss(zc, n-m);
3803 368676 : init_modular_small(&S);
3804 368675 : if (n <= m) { nn = n; mm = m; } else { nn = m; mm = n; }
3805 368675 : imax = (nn < 16)? 1: (nn < 64)? 2: 3; /* heuristic */
3806 :
3807 : for(;;)
3808 0 : {
3809 : GEN row, col, M, KM, IM, RHS, X, cX;
3810 : long rk;
3811 394927 : for (av = avma, i = 0;; set_avma(av), i++)
3812 26254 : {
3813 394927 : ulong p = u_forprime_next(&S);
3814 : long rp;
3815 394927 : if (!p) pari_err_OVERFLOW("ZM_pivots [ran out of primes]");
3816 394927 : d = Flm_pivots(ZM_to_Flm(M0, p), p, &rp, 1);
3817 394932 : if (rp == rmin) { rbest = rp; goto END; } /* maximal rank, return */
3818 50980 : if (rp < rbest) { /* save best r so far */
3819 24728 : rbest = rp;
3820 24728 : guncloneNULL(dbest);
3821 24728 : dbest = gclone(d);
3822 24728 : if (beenthere) break;
3823 : }
3824 50980 : if (!beenthere && i >= imax) break;
3825 : }
3826 24726 : beenthere = 1;
3827 : /* Dubious case: there is (probably) a non trivial kernel */
3828 24726 : indexrank_all(m,n, rbest, dbest, &row, &col);
3829 24726 : M = rowpermute(vecpermute(M0, col), row);
3830 24725 : rk = n - rbest; /* (probable) dimension of image */
3831 24725 : if (n > m) M = shallowtrans(M);
3832 24725 : IM = vecslice(M,1,rk);
3833 24725 : KM = vecslice(M,rk+1, nn);
3834 24725 : M = rowslice(IM, 1,rk); /* square maximal rank */
3835 24725 : X = ZM_gauss(M, rowslice(KM, 1,rk));
3836 24726 : RHS = rowslice(KM,rk+1,mm);
3837 24726 : M = rowslice(IM,rk+1,mm);
3838 24726 : X = Q_remove_denom(X, &cX);
3839 24726 : if (cX) RHS = ZM_Z_mul(RHS, cX);
3840 24726 : if (ZM_equal(ZM_mul(M, X), RHS)) { d = vecsmall_copy(dbest); goto END; }
3841 0 : set_avma(av);
3842 : }
3843 368677 : END:
3844 368677 : *rr = rbest; guncloneNULL(dbest);
3845 368677 : return gerepileuptoleaf(av0, d);
3846 : }
3847 :
3848 : /* set *pr = dim Ker x */
3849 : static GEN
3850 58575 : gauss_pivot(GEN x, long *pr) {
3851 : GEN data;
3852 58575 : pivot_fun pivot = get_pivot_fun(x, x, &data);
3853 58575 : return RgM_pivots(x, data, pr, pivot);
3854 : }
3855 :
3856 : /* compute ker(x), x0 is a reference point when guessing whether x[i,j] ~ 0
3857 : * (iff x[i,j] << x0[i,j]) */
3858 : static GEN
3859 1042 : ker_aux(GEN x, GEN x0)
3860 : {
3861 1042 : pari_sp av = avma;
3862 : GEN d,y;
3863 : long i,j,k,r,n;
3864 :
3865 1042 : x = gauss_pivot_ker(x,x0,&d,&r);
3866 1042 : if (!r) { set_avma(av); return cgetg(1,t_MAT); }
3867 1001 : n = lg(x)-1; y=cgetg(r+1,t_MAT);
3868 2100 : for (j=k=1; j<=r; j++,k++)
3869 : {
3870 1099 : GEN p = cgetg(n+1,t_COL);
3871 :
3872 4662 : gel(y,j) = p; while (d[k]) k++;
3873 5026 : for (i=1; i<k; i++)
3874 3927 : if (d[i])
3875 : {
3876 3745 : GEN p1=gcoeff(x,d[i],k);
3877 3745 : gel(p,i) = gcopy(p1); gunclone(p1);
3878 : }
3879 : else
3880 182 : gel(p,i) = gen_0;
3881 1869 : gel(p,k) = gen_1; for (i=k+1; i<=n; i++) gel(p,i) = gen_0;
3882 : }
3883 1001 : return gerepileupto(av,y);
3884 : }
3885 :
3886 : static GEN
3887 77 : RgM_ker_FpM(GEN x, GEN p)
3888 : {
3889 77 : pari_sp av = avma;
3890 : ulong pp;
3891 77 : x = RgM_Fp_init(x, p, &pp);
3892 77 : switch(pp)
3893 : {
3894 35 : case 0: x = FpM_to_mod(FpM_ker_gen(x,p,0),p); break;
3895 7 : case 2: x = F2m_to_mod(F2m_ker_sp(x,0)); break;
3896 35 : default:x = Flm_to_mod(Flm_ker_sp(x,pp,0), pp); break;
3897 : }
3898 77 : return gerepileupto(av, x);
3899 : }
3900 :
3901 : static GEN
3902 91 : RgM_ker_FqM(GEN x, GEN pol, GEN p)
3903 : {
3904 91 : pari_sp av = avma;
3905 91 : GEN b, T = RgX_to_FpX(pol, p);
3906 91 : if (signe(T) == 0) pari_err_OP("ker",x,pol);
3907 84 : b = FqM_ker(RgM_to_FqM(x, T, p), T, p);
3908 84 : return gerepileupto(av, FqM_to_mod(b, T, p));
3909 : }
3910 :
3911 : #define code(t1,t2) ((t1 << 6) | t2)
3912 : static GEN
3913 8358 : RgM_ker_fast(GEN x)
3914 : {
3915 : GEN p, pol;
3916 : long pa;
3917 8358 : long t = RgM_type(x, &p,&pol,&pa);
3918 8358 : switch(t)
3919 : {
3920 7357 : case t_INT: /* fall through */
3921 7357 : case t_FRAC: return QM_ker(x);
3922 77 : case t_FFELT: return FFM_ker(x, pol);
3923 77 : case t_INTMOD: return RgM_ker_FpM(x, p);
3924 91 : case code(t_POLMOD, t_INTMOD):
3925 91 : return RgM_ker_FqM(x, pol, p);
3926 756 : default: return NULL;
3927 : }
3928 : }
3929 : #undef code
3930 :
3931 : GEN
3932 8358 : ker(GEN x)
3933 : {
3934 8358 : GEN b = RgM_ker_fast(x);
3935 8351 : if (b) return b;
3936 756 : return ker_aux(x,x);
3937 : }
3938 :
3939 : GEN
3940 46214 : matker0(GEN x,long flag)
3941 : {
3942 46214 : if (typ(x)!=t_MAT) pari_err_TYPE("matker",x);
3943 46214 : if (!flag) return ker(x);
3944 45934 : RgM_check_ZM(x, "matker");
3945 45934 : return ZM_ker(x);
3946 : }
3947 :
3948 : static GEN
3949 63 : RgM_image_FpM(GEN x, GEN p)
3950 : {
3951 63 : pari_sp av = avma;
3952 : ulong pp;
3953 63 : x = RgM_Fp_init(x, p, &pp);
3954 63 : switch(pp)
3955 : {
3956 28 : case 0: x = FpM_to_mod(FpM_image(x,p),p); break;
3957 7 : case 2: x = F2m_to_mod(F2m_image(x)); break;
3958 28 : default:x = Flm_to_mod(Flm_image(x,pp), pp); break;
3959 : }
3960 63 : return gerepileupto(av, x);
3961 : }
3962 :
3963 : static GEN
3964 35 : RgM_image_FqM(GEN x, GEN pol, GEN p)
3965 : {
3966 35 : pari_sp av = avma;
3967 35 : GEN b, T = RgX_to_FpX(pol, p);
3968 35 : if (signe(T) == 0) pari_err_OP("image",x,pol);
3969 28 : b = FqM_image(RgM_to_FqM(x, T, p), T, p);
3970 28 : return gerepileupto(av, FqM_to_mod(b, T, p));
3971 : }
3972 :
3973 : GEN
3974 5971 : QM_image_shallow(GEN A)
3975 : {
3976 5971 : A = vec_Q_primpart(A);
3977 5971 : return vecpermute(A, ZM_indeximage(A));
3978 : }
3979 : GEN
3980 5201 : QM_image(GEN A)
3981 : {
3982 5201 : pari_sp av = avma;
3983 5201 : return gerepilecopy(av, QM_image_shallow(A));
3984 : }
3985 :
3986 : #define code(t1,t2) ((t1 << 6) | t2)
3987 : static GEN
3988 5362 : RgM_image_fast(GEN x)
3989 : {
3990 : GEN p, pol;
3991 : long pa;
3992 5362 : long t = RgM_type(x, &p,&pol,&pa);
3993 5362 : switch(t)
3994 : {
3995 5201 : case t_INT: /* fall through */
3996 5201 : case t_FRAC: return QM_image(x);
3997 49 : case t_FFELT: return FFM_image(x, pol);
3998 63 : case t_INTMOD: return RgM_image_FpM(x, p);
3999 35 : case code(t_POLMOD, t_INTMOD):
4000 35 : return RgM_image_FqM(x, pol, p);
4001 14 : default: return NULL;
4002 : }
4003 : }
4004 : #undef code
4005 :
4006 : GEN
4007 5362 : image(GEN x)
4008 : {
4009 : GEN d, M;
4010 : long r;
4011 :
4012 5362 : if (typ(x)!=t_MAT) pari_err_TYPE("matimage",x);
4013 5362 : M = RgM_image_fast(x);
4014 5355 : if (M) return M;
4015 14 : d = gauss_pivot(x,&r); /* d left on stack for efficiency */
4016 14 : return image_from_pivot(x,d,r);
4017 : }
4018 :
4019 : static GEN
4020 84 : imagecompl_aux(GEN x, GEN(*PIVOT)(GEN,long*))
4021 : {
4022 84 : pari_sp av = avma;
4023 : GEN d,y;
4024 : long j,i,r;
4025 :
4026 84 : if (typ(x)!=t_MAT) pari_err_TYPE("imagecompl",x);
4027 84 : (void)new_chunk(lg(x) * 4 + 1); /* HACK */
4028 84 : d = PIVOT(x,&r); /* if (!d) then r = 0 */
4029 84 : set_avma(av); y = cgetg(r+1,t_VECSMALL);
4030 126 : for (i=j=1; j<=r; i++)
4031 42 : if (!d[i]) y[j++] = i;
4032 84 : return y;
4033 : }
4034 : GEN
4035 84 : imagecompl(GEN x) { return imagecompl_aux(x, &gauss_pivot); }
4036 : GEN
4037 0 : ZM_imagecompl(GEN x) { return imagecompl_aux(x, &ZM_pivots); }
4038 :
4039 : static GEN
4040 28 : RgM_RgC_invimage_FpC(GEN A, GEN y, GEN p)
4041 : {
4042 28 : pari_sp av = avma;
4043 : ulong pp;
4044 : GEN x;
4045 28 : A = RgM_Fp_init(A,p,&pp);
4046 28 : switch(pp)
4047 : {
4048 7 : case 0:
4049 7 : y = RgC_to_FpC(y,p);
4050 7 : x = FpM_FpC_invimage(A, y, p);
4051 7 : return x ? gerepileupto(av, FpC_to_mod(x,p)): NULL;
4052 7 : case 2:
4053 7 : y = RgV_to_F2v(y);
4054 7 : x = F2m_F2c_invimage(A, y);
4055 7 : return x ? gerepileupto(av, F2c_to_mod(x)): NULL;
4056 14 : default:
4057 14 : y = RgV_to_Flv(y,pp);
4058 14 : x = Flm_Flc_invimage(A, y, pp);
4059 14 : return x ? gerepileupto(av, Flc_to_mod(x,pp)): NULL;
4060 : }
4061 : }
4062 :
4063 : static GEN
4064 2072 : RgM_RgC_invimage_fast(GEN x, GEN y)
4065 : {
4066 : GEN p, pol;
4067 : long pa;
4068 2072 : long t = RgM_RgC_type(x, y, &p,&pol,&pa);
4069 2072 : switch(t)
4070 : {
4071 28 : case t_INTMOD: return RgM_RgC_invimage_FpC(x, y, p);
4072 63 : case t_FFELT: return FFM_FFC_invimage(x, y, pol);
4073 1981 : default: return gen_0;
4074 : }
4075 : }
4076 :
4077 : GEN
4078 2177 : RgM_RgC_invimage(GEN A, GEN y)
4079 : {
4080 2177 : pari_sp av = avma;
4081 2177 : long i, l = lg(A);
4082 : GEN M, x, t;
4083 2177 : if (l==1) return NULL;
4084 2072 : if (lg(y) != lgcols(A)) pari_err_DIM("inverseimage");
4085 2072 : M = RgM_RgC_invimage_fast(A, y);
4086 2072 : if (!M) return gc_NULL(av);
4087 2051 : if (M != gen_0) return M;
4088 1981 : M = ker(shallowconcat(A, y));
4089 1981 : i = lg(M)-1;
4090 1981 : if (!i) return gc_NULL(av);
4091 :
4092 1722 : x = gel(M,i); t = gel(x,l);
4093 1722 : if (gequal0(t)) return gc_NULL(av);
4094 :
4095 1687 : t = gneg_i(t); setlg(x,l);
4096 1687 : return gerepileupto(av, RgC_Rg_div(x, t));
4097 : }
4098 :
4099 : /* Return X such that m X = v (t_COL or t_MAT), resp. an empty t_COL / t_MAT
4100 : * if no solution exist */
4101 : GEN
4102 2359 : inverseimage(GEN m, GEN v)
4103 : {
4104 : GEN y;
4105 2359 : if (typ(m)!=t_MAT) pari_err_TYPE("inverseimage",m);
4106 2359 : switch(typ(v))
4107 : {
4108 2121 : case t_COL:
4109 2121 : y = RgM_RgC_invimage(m,v);
4110 2121 : return y? y: cgetg(1,t_COL);
4111 238 : case t_MAT:
4112 238 : y = RgM_invimage(m, v);
4113 238 : return y? y: cgetg(1,t_MAT);
4114 : }
4115 0 : pari_err_TYPE("inverseimage",v);
4116 : return NULL;/*LCOV_EXCL_LINE*/
4117 : }
4118 :
4119 : static GEN
4120 84 : RgM_invimage_FpM(GEN A, GEN B, GEN p)
4121 : {
4122 84 : pari_sp av = avma;
4123 : ulong pp;
4124 : GEN x;
4125 84 : A = RgM_Fp_init(A,p,&pp);
4126 84 : switch(pp)
4127 : {
4128 35 : case 0:
4129 35 : B = RgM_to_FpM(B,p);
4130 35 : x = FpM_invimage_gen(A, B, p);
4131 35 : return x ? gerepileupto(av, FpM_to_mod(x, p)): x;
4132 7 : case 2:
4133 7 : B = RgM_to_F2m(B);
4134 7 : x = F2m_invimage_i(A, B);
4135 7 : return x ? gerepileupto(av, F2m_to_mod(x)): x;
4136 42 : default:
4137 42 : B = RgM_to_Flm(B,pp);
4138 42 : x = Flm_invimage_i(A, B, pp);
4139 42 : return x ? gerepileupto(av, Flm_to_mod(x, pp)): x;
4140 : }
4141 : }
4142 :
4143 : static GEN
4144 252 : RgM_invimage_fast(GEN x, GEN y)
4145 : {
4146 : GEN p, pol;
4147 : long pa;
4148 252 : long t = RgM_type2(x, y, &p,&pol,&pa);
4149 252 : switch(t)
4150 : {
4151 84 : case t_INTMOD: return RgM_invimage_FpM(x, y, p);
4152 105 : case t_FFELT: return FFM_invimage(x, y, pol);
4153 63 : default: return gen_0;
4154 : }
4155 : }
4156 :
4157 : /* find Z such that A Z = B. Return NULL if no solution */
4158 : GEN
4159 252 : RgM_invimage(GEN A, GEN B)
4160 : {
4161 252 : pari_sp av = avma;
4162 : GEN d, x, X, Y;
4163 252 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
4164 252 : X = RgM_invimage_fast(A, B);
4165 252 : if (!X) return gc_NULL(av);
4166 140 : if (X != gen_0) return X;
4167 63 : x = ker(shallowconcat(RgM_neg(A), B));
4168 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
4169 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
4170 : * Y has at least nB columns and full rank */
4171 63 : nY = lg(x)-1;
4172 63 : if (nY < nB) return gc_NULL(av);
4173 49 : Y = rowslice(x, nA+1, nA+nB); /* nB rows */
4174 49 : d = cgetg(nB+1, t_VECSMALL);
4175 441 : for (i = nB, j = nY; i >= 1; i--, j--)
4176 : {
4177 546 : for (; j>=1; j--)
4178 532 : if (!gequal0(gcoeff(Y,i,j))) { d[i] = j; break; }
4179 406 : if (!j) return gc_NULL(av);
4180 : }
4181 : /* reduce to the case Y square, upper triangular with 1s on diagonal */
4182 35 : Y = vecpermute(Y, d);
4183 35 : x = vecpermute(x, d);
4184 35 : X = rowslice(x, 1, nA);
4185 35 : return gerepileupto(av, RgM_mul(X, RgM_inv_upper(Y)));
4186 : }
4187 :
4188 : static GEN
4189 70 : RgM_suppl_FpM(GEN x, GEN p)
4190 : {
4191 70 : pari_sp av = avma;
4192 : ulong pp;
4193 70 : x = RgM_Fp_init(x, p, &pp);
4194 70 : switch(pp)
4195 : {
4196 21 : case 0: x = FpM_to_mod(FpM_suppl(x,p), p); break;
4197 14 : case 2: x = F2m_to_mod(F2m_suppl(x)); break;
4198 35 : default:x = Flm_to_mod(Flm_suppl(x,pp), pp); break;
4199 : }
4200 70 : return gerepileupto(av, x);
4201 : }
4202 :
4203 : static GEN
4204 175 : RgM_suppl_fast(GEN x)
4205 : {
4206 : GEN p, pol;
4207 : long pa;
4208 175 : long t = RgM_type(x,&p,&pol,&pa);
4209 175 : switch(t)
4210 : {
4211 70 : case t_INTMOD: return RgM_suppl_FpM(x, p);
4212 35 : case t_FFELT: return FFM_suppl(x, pol);
4213 70 : default: return NULL;
4214 : }
4215 : }
4216 :
4217 : /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
4218 : * whose first k columns are given by x. If rank(x) < k, undefined result. */
4219 : GEN
4220 175 : suppl(GEN x)
4221 : {
4222 175 : pari_sp av = avma;
4223 : GEN d, M;
4224 : long r;
4225 175 : if (typ(x)!=t_MAT) pari_err_TYPE("suppl",x);
4226 175 : M = RgM_suppl_fast(x);
4227 175 : if (M) return M;
4228 70 : init_suppl(x);
4229 70 : d = gauss_pivot(x,&r);
4230 70 : set_avma(av); return get_suppl(x,d,nbrows(x),r,&col_ei);
4231 : }
4232 :
4233 : GEN
4234 7 : image2(GEN x)
4235 : {
4236 7 : pari_sp av = avma;
4237 : long k, n, i;
4238 : GEN A, B;
4239 :
4240 7 : if (typ(x)!=t_MAT) pari_err_TYPE("image2",x);
4241 7 : if (lg(x) == 1) return cgetg(1,t_MAT);
4242 7 : A = ker(x); k = lg(A)-1;
4243 7 : if (!k) { set_avma(av); return gcopy(x); }
4244 7 : A = suppl(A); n = lg(A)-1;
4245 7 : B = cgetg(n-k+1, t_MAT);
4246 21 : for (i = k+1; i <= n; i++) gel(B,i-k) = RgM_RgC_mul(x, gel(A,i));
4247 7 : return gerepileupto(av, B);
4248 : }
4249 :
4250 : GEN
4251 210 : matimage0(GEN x,long flag)
4252 : {
4253 210 : switch(flag)
4254 : {
4255 203 : case 0: return image(x);
4256 7 : case 1: return image2(x);
4257 0 : default: pari_err_FLAG("matimage");
4258 : }
4259 : return NULL; /* LCOV_EXCL_LINE */
4260 : }
4261 :
4262 : static long
4263 126 : RgM_rank_FpM(GEN x, GEN p)
4264 : {
4265 126 : pari_sp av = avma;
4266 : ulong pp;
4267 : long r;
4268 126 : x = RgM_Fp_init(x,p,&pp);
4269 126 : switch(pp)
4270 : {
4271 28 : case 0: r = FpM_rank(x,p); break;
4272 63 : case 2: r = F2m_rank(x); break;
4273 35 : default:r = Flm_rank(x,pp); break;
4274 : }
4275 126 : return gc_long(av, r);
4276 : }
4277 :
4278 : static long
4279 49 : RgM_rank_FqM(GEN x, GEN pol, GEN p)
4280 : {
4281 49 : pari_sp av = avma;
4282 : long r;
4283 49 : GEN T = RgX_to_FpX(pol, p);
4284 49 : if (signe(T) == 0) pari_err_OP("rank",x,pol);
4285 42 : r = FqM_rank(RgM_to_FqM(x, T, p), T, p);
4286 42 : return gc_long(av,r);
4287 : }
4288 :
4289 : #define code(t1,t2) ((t1 << 6) | t2)
4290 : static long
4291 294 : RgM_rank_fast(GEN x)
4292 : {
4293 : GEN p, pol;
4294 : long pa;
4295 294 : long t = RgM_type(x,&p,&pol,&pa);
4296 294 : switch(t)
4297 : {
4298 42 : case t_INT: return ZM_rank(x);
4299 0 : case t_FRAC: return QM_rank(x);
4300 126 : case t_INTMOD: return RgM_rank_FpM(x, p);
4301 70 : case t_FFELT: return FFM_rank(x, pol);
4302 49 : case code(t_POLMOD, t_INTMOD):
4303 49 : return RgM_rank_FqM(x, pol, p);
4304 7 : default: return -1;
4305 : }
4306 : }
4307 : #undef code
4308 :
4309 : long
4310 294 : rank(GEN x)
4311 : {
4312 294 : pari_sp av = avma;
4313 : long r;
4314 :
4315 294 : if (typ(x)!=t_MAT) pari_err_TYPE("rank",x);
4316 294 : r = RgM_rank_fast(x);
4317 287 : if (r >= 0) return r;
4318 7 : (void)gauss_pivot(x, &r);
4319 7 : return gc_long(av, lg(x)-1 - r);
4320 : }
4321 :
4322 : /* d a t_VECSMALL of integers in 1..n. Return the vector of the d[i]
4323 : * followed by the missing indices */
4324 : static GEN
4325 49451 : perm_complete(GEN d, long n)
4326 : {
4327 49451 : GEN y = cgetg(n+1, t_VECSMALL);
4328 49451 : long i, j = 1, k = n, l = lg(d);
4329 49451 : pari_sp av = avma;
4330 49451 : char *T = stack_calloc(n+1);
4331 238495 : for (i = 1; i < l; i++) T[d[i]] = 1;
4332 469122 : for (i = 1; i <= n; i++)
4333 419673 : if (T[i]) y[j++] = i; else y[k--] = i;
4334 49449 : return gc_const(av, y);
4335 : }
4336 :
4337 : /* n = dim x, r = dim Ker(x), d from gauss_pivot */
4338 : static GEN
4339 5971 : indeximage0(long n, long r, GEN d)
4340 : {
4341 : long i, j;
4342 : GEN v;
4343 :
4344 5971 : r = n - r; /* now r = dim Im(x) */
4345 5971 : v = cgetg(r+1,t_VECSMALL);
4346 32571 : if (d) for (i=j=1; j<=n; j++)
4347 26600 : if (d[j]) v[i++] = j;
4348 5971 : return v;
4349 : }
4350 : /* x an m x n t_MAT, n > 0, r = dim Ker(x), d from gauss_pivot */
4351 : static void
4352 24726 : indexrank_all(long m, long n, long r, GEN d, GEN *prow, GEN *pcol)
4353 : {
4354 24726 : GEN IR = indexrank0(n, r, d);
4355 24726 : *prow = perm_complete(gel(IR,1), m);
4356 24726 : *pcol = perm_complete(gel(IR,2), n);
4357 24726 : }
4358 :
4359 : static GEN
4360 28 : RgM_indexrank_FpM(GEN x, GEN p)
4361 : {
4362 28 : pari_sp av = avma;
4363 : ulong pp;
4364 : GEN r;
4365 28 : x = RgM_Fp_init(x,p,&pp);
4366 28 : switch(pp)
4367 : {
4368 7 : case 0: r = FpM_indexrank(x,p); break;
4369 7 : case 2: r = F2m_indexrank(x); break;
4370 14 : default: r = Flm_indexrank(x,pp); break;
4371 : }
4372 28 : return gerepileupto(av, r);
4373 : }
4374 :
4375 : static GEN
4376 0 : RgM_indexrank_FqM(GEN x, GEN pol, GEN p)
4377 : {
4378 0 : pari_sp av = avma;
4379 0 : GEN r, T = RgX_to_FpX(pol, p);
4380 0 : if (signe(T) == 0) pari_err_OP("indexrank",x,pol);
4381 0 : r = FqM_indexrank(RgM_to_FqM(x, T, p), T, p);
4382 0 : return gerepileupto(av, r);
4383 : }
4384 :
4385 : #define code(t1,t2) ((t1 << 6) | t2)
4386 : static GEN
4387 59987 : RgM_indexrank_fast(GEN x)
4388 : {
4389 : GEN p, pol;
4390 : long pa;
4391 59987 : long t = RgM_type(x,&p,&pol,&pa);
4392 59989 : switch(t)
4393 : {
4394 392 : case t_INT: return ZM_indexrank(x);
4395 1148 : case t_FRAC: return QM_indexrank(x);
4396 28 : case t_INTMOD: return RgM_indexrank_FpM(x, p);
4397 21 : case t_FFELT: return FFM_indexrank(x, pol);
4398 0 : case code(t_POLMOD, t_INTMOD):
4399 0 : return RgM_indexrank_FqM(x, pol, p);
4400 58400 : default: return NULL;
4401 : }
4402 : }
4403 : #undef code
4404 :
4405 : GEN
4406 59987 : indexrank(GEN x)
4407 : {
4408 : pari_sp av;
4409 : long r;
4410 : GEN d;
4411 59987 : if (typ(x)!=t_MAT) pari_err_TYPE("indexrank",x);
4412 59987 : d = RgM_indexrank_fast(x);
4413 59989 : if (d) return d;
4414 58400 : av = avma;
4415 58400 : init_indexrank(x);
4416 58400 : d = gauss_pivot(x, &r);
4417 58398 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4418 : }
4419 :
4420 : GEN
4421 5971 : ZM_indeximage(GEN x) {
4422 5971 : pari_sp av = avma;
4423 : long r;
4424 : GEN d;
4425 5971 : init_indexrank(x);
4426 5971 : d = ZM_pivots(x,&r);
4427 5971 : set_avma(av); return indeximage0(lg(x)-1, r, d);
4428 : }
4429 : long
4430 99084 : ZM_rank(GEN x) {
4431 99084 : pari_sp av = avma;
4432 : long r;
4433 99084 : (void)ZM_pivots(x,&r);
4434 99082 : return gc_long(av, lg(x)-1-r);
4435 : }
4436 : GEN
4437 25333 : ZM_indexrank(GEN x) {
4438 25333 : pari_sp av = avma;
4439 : long r;
4440 : GEN d;
4441 25333 : init_indexrank(x);
4442 25333 : d = ZM_pivots(x,&r);
4443 25333 : set_avma(av); return indexrank0(lg(x)-1, r, d);
4444 : }
4445 :
4446 : long
4447 0 : QM_rank(GEN x)
4448 : {
4449 0 : pari_sp av = avma;
4450 0 : long r = ZM_rank(Q_primpart(x));
4451 0 : set_avma(av);
4452 0 : return r;
4453 : }
4454 :
4455 : GEN
4456 1148 : QM_indexrank(GEN x)
4457 : {
4458 1148 : pari_sp av = avma;
4459 1148 : GEN r = ZM_indexrank(Q_primpart(x));
4460 1148 : return gerepileupto(av, r);
4461 : }
4462 :
4463 : /*******************************************************************/
4464 : /* */
4465 : /* ZabM */
4466 : /* */
4467 : /*******************************************************************/
4468 :
4469 : static GEN
4470 1276 : FpXM_ratlift(GEN a, GEN q)
4471 : {
4472 : GEN B, y;
4473 1276 : long i, j, l = lg(a), n;
4474 1276 : B = sqrti(shifti(q,-1));
4475 1276 : y = cgetg(l, t_MAT);
4476 1276 : if (l==1) return y;
4477 1276 : n = lgcols(a);
4478 3059 : for (i=1; i<l; i++)
4479 : {
4480 2404 : GEN yi = cgetg(n, t_COL);
4481 32311 : for (j=1; j<n; j++)
4482 : {
4483 30528 : GEN v = FpX_ratlift(gmael(a,i,j), q, B, B, NULL);
4484 30528 : if (!v) return NULL;
4485 29907 : gel(yi, j) = RgX_renormalize(v);
4486 : }
4487 1783 : gel(y,i) = yi;
4488 : }
4489 655 : return y;
4490 : }
4491 :
4492 : static GEN
4493 4449 : FlmV_recover_pre(GEN a, GEN M, ulong p, ulong pi, long sv)
4494 : {
4495 4449 : GEN a1 = gel(a,1);
4496 4449 : long i, j, k, l = lg(a1), n, lM = lg(M);
4497 4449 : GEN v = cgetg(lM, t_VECSMALL);
4498 4449 : GEN y = cgetg(l, t_MAT);
4499 4449 : if (l==1) return y;
4500 4449 : n = lgcols(a1);
4501 22380 : for (i=1; i<l; i++)
4502 : {
4503 17932 : GEN yi = cgetg(n, t_COL);
4504 347069 : for (j=1; j<n; j++)
4505 : {
4506 4673155 : for (k=1; k<lM; k++) uel(v,k) = umael(gel(a,k),i,j);
4507 329138 : gel(yi, j) = Flm_Flc_mul_pre_Flx(M, v, p, pi, sv);
4508 : }
4509 17931 : gel(y,i) = yi;
4510 : }
4511 4448 : return y;
4512 : }
4513 :
4514 : static GEN
4515 0 : FlkM_inv(GEN M, GEN P, ulong p)
4516 : {
4517 0 : ulong pi = get_Fl_red(p);
4518 0 : GEN R = Flx_roots(P, p);
4519 0 : long l = lg(R), i;
4520 0 : GEN W = Flv_invVandermonde(R, 1UL, p);
4521 0 : GEN V = cgetg(l, t_VEC);
4522 0 : for(i=1; i<l; i++)
4523 : {
4524 0 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, pi);
4525 0 : GEN H = Flm_inv_sp(FlxM_eval_powers_pre(M, pows, p, pi), NULL, p);
4526 0 : if (!H) return NULL;
4527 0 : gel(V, i) = H;
4528 : }
4529 0 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4530 : }
4531 :
4532 : static GEN
4533 3173 : FlkM_adjoint(GEN M, GEN P, ulong p)
4534 : {
4535 3173 : ulong pi = get_Fl_red(p);
4536 3173 : GEN R = Flx_roots(P, p);
4537 3173 : long l = lg(R), i;
4538 3173 : GEN W = Flv_invVandermonde(R, 1UL, p);
4539 3173 : GEN V = cgetg(l, t_VEC);
4540 15436 : for(i=1; i<l; i++)
4541 : {
4542 12263 : GEN pows = Fl_powers_pre(uel(R,i), degpol(P), p, pi);
4543 12262 : gel(V, i) = Flm_adjoint(FlxM_eval_powers_pre(M, pows, p, pi), p);
4544 : }
4545 3173 : return FlmV_recover_pre(V, W, p, pi, P[1]);
4546 : }
4547 :
4548 : static GEN
4549 1950 : ZabM_inv_slice(GEN A, GEN Q, GEN P, GEN *mod)
4550 : {
4551 1950 : pari_sp av = avma;
4552 1950 : long i, n = lg(P)-1, w = varn(Q);
4553 : GEN H, T;
4554 1950 : if (n == 1)
4555 : {
4556 1520 : ulong p = uel(P,1);
4557 1520 : GEN Qp = ZX_to_Flx(Q, p);
4558 1520 : GEN Ap = ZXM_to_FlxM(A, p, get_Flx_var(Qp));
4559 1520 : GEN Hp = FlkM_adjoint(Ap, Qp, p);
4560 1520 : Hp = gerepileupto(av, FlxM_to_ZXM(Hp));
4561 1520 : *mod = utoipos(p); return Hp;
4562 : }
4563 430 : T = ZV_producttree(P);
4564 430 : A = ZXM_nv_mod_tree(A, P, T, w);
4565 430 : Q = ZX_nv_mod_tree(Q, P, T);
4566 430 : H = cgetg(n+1, t_VEC);
4567 2083 : for(i=1; i <= n; i++)
4568 : {
4569 1653 : ulong p = P[i];
4570 1653 : GEN a = gel(A,i), q = gel(Q, i);
4571 1653 : gel(H,i) = FlkM_adjoint(a, q, p);
4572 : }
4573 430 : H = nxMV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P,T));
4574 430 : *mod = gmael(T, lg(T)-1, 1);
4575 430 : gerepileall(av, 2, &H, mod);
4576 430 : return H;
4577 : }
4578 :
4579 : GEN
4580 1950 : ZabM_inv_worker(GEN P, GEN A, GEN Q)
4581 : {
4582 1950 : GEN V = cgetg(3, t_VEC);
4583 1950 : gel(V,1) = ZabM_inv_slice(A, Q, P, &gel(V,2));
4584 1950 : return V;
4585 : }
4586 :
4587 : static GEN
4588 5411 : vecnorml1(GEN a)
4589 : {
4590 : long i, l;
4591 5411 : GEN g = cgetg_copy(a, &l);
4592 60046 : for (i=1; i<l; i++)
4593 54635 : gel(g, i) = gnorml1_fake(gel(a,i));
4594 5411 : return g;
4595 : }
4596 :
4597 : static GEN
4598 1792 : ZabM_true_Hadamard(GEN a)
4599 : {
4600 1792 : pari_sp av = avma;
4601 1792 : long n = lg(a)-1, i;
4602 : GEN B;
4603 1792 : if (n == 0) return gen_1;
4604 1792 : if (n == 1) return gnorml1_fake(gcoeff(a,1,1));
4605 1148 : B = gen_1;
4606 6559 : for (i = 1; i <= n; i++)
4607 5411 : B = gmul(B, gnorml2(RgC_gtofp(vecnorml1(gel(a,i)),DEFAULTPREC)));
4608 1148 : return gerepileuptoint(av, ceil_safe(sqrtr_abs(B)));
4609 : }
4610 :
4611 : GEN
4612 1792 : ZabM_inv(GEN A, GEN Q, long n, GEN *pt_den)
4613 : {
4614 1792 : pari_sp av = avma;
4615 : forprime_t S;
4616 : GEN bnd, H, D, d, mod, worker;
4617 1792 : if (lg(A) == 1)
4618 : {
4619 0 : if (pt_den) *pt_den = gen_1;
4620 0 : return cgetg(1, t_MAT);
4621 : }
4622 1792 : bnd = ZabM_true_Hadamard(A);
4623 1792 : worker = snm_closure(is_entry("_ZabM_inv_worker"), mkvec2(A, Q));
4624 1792 : u_forprime_arith_init(&S, HIGHBIT+1, ULONG_MAX, 1, n);
4625 1792 : H = gen_crt("ZabM_inv", worker, &S, NULL, expi(bnd), 0, &mod,
4626 : nxMV_chinese_center, FpXM_center);
4627 1792 : D = RgMrow_RgC_mul(H, gel(A,1), 1);
4628 1792 : D = ZX_rem(D, Q);
4629 1792 : d = Z_content(mkvec2(H, D));
4630 1792 : if (d)
4631 : {
4632 511 : D = ZX_Z_divexact(D, d);
4633 511 : H = Q_div_to_int(H, d);
4634 : }
4635 1792 : if (pt_den)
4636 : {
4637 1792 : gerepileall(av, 2, &H, &D);
4638 1792 : *pt_den = D; return H;
4639 : }
4640 0 : return gerepileupto(av, H);
4641 : }
4642 :
4643 : GEN
4644 0 : ZabM_inv_ratlift(GEN M, GEN P, long n, GEN *pden)
4645 : {
4646 0 : pari_sp av2, av = avma;
4647 : GEN q, H;
4648 0 : ulong m = LONG_MAX>>1;
4649 0 : ulong p= 1 + m - (m % n);
4650 0 : long lM = lg(M);
4651 0 : if (lM == 1) { *pden = gen_1; return cgetg(1,t_MAT); }
4652 :
4653 0 : av2 = avma;
4654 0 : H = NULL;
4655 : for(;;)
4656 0 : {
4657 : GEN Hp, Pp, Mp, Hr;
4658 0 : do p += n; while(!uisprime(p));
4659 0 : Pp = ZX_to_Flx(P, p);
4660 0 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4661 0 : Hp = FlkM_inv(Mp, Pp, p);
4662 0 : if (!Hp) continue;
4663 0 : if (!H)
4664 : {
4665 0 : H = ZXM_init_CRT(Hp, degpol(P)-1, p);
4666 0 : q = utoipos(p);
4667 : }
4668 : else
4669 0 : ZXM_incremental_CRT(&H, Hp, &q, p);
4670 0 : Hr = FpXM_ratlift(H, q);
4671 0 : if (DEBUGLEVEL>5) err_printf("ZabM_inv mod %ld (ratlift=%ld)\n", p,!!Hr);
4672 0 : if (Hr) {/* DONE ? */
4673 0 : GEN Hl = Q_remove_denom(Hr, pden);
4674 0 : GEN MH = ZXQM_mul(Hl, M, P);
4675 0 : if (*pden)
4676 0 : { if (RgM_isscalar(MH, *pden)) { H = Hl; break; }}
4677 : else
4678 0 : { if (RgM_isidentity(MH)) { H = Hl; *pden = gen_1; break; } }
4679 : }
4680 :
4681 0 : if (gc_needed(av,2))
4682 : {
4683 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_inv");
4684 0 : gerepileall(av2, 2, &H, &q);
4685 : }
4686 : }
4687 0 : gerepileall(av, 2, &H, pden);
4688 0 : return H;
4689 : }
4690 :
4691 : static GEN
4692 1276 : FlkM_ker(GEN M, GEN P, ulong p)
4693 : {
4694 1276 : ulong pi = get_Fl_red(p);
4695 1276 : GEN R = Flx_roots(P, p);
4696 1276 : long l = lg(R), i, dP = degpol(P), r;
4697 : GEN M1, K, D;
4698 1276 : GEN W = Flv_invVandermonde(R, 1UL, p);
4699 1276 : GEN V = cgetg(l, t_VEC);
4700 1276 : M1 = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,1), dP, p, pi), p, pi);
4701 1276 : K = Flm_ker_sp(M1, p, 2);
4702 1276 : r = lg(gel(K,1)); D = gel(K,2);
4703 1276 : gel(V, 1) = gel(K,1);
4704 2652 : for(i=2; i<l; i++)
4705 : {
4706 1376 : GEN Mi = FlxM_eval_powers_pre(M, Fl_powers_pre(uel(R,i), dP, p, pi), p, pi);
4707 1376 : GEN K = Flm_ker_sp(Mi, p, 2);
4708 1376 : if (lg(gel(K,1)) != r || !zv_equal(D, gel(K,2))) return NULL;
4709 1376 : gel(V, i) = gel(K,1);
4710 : }
4711 1276 : return mkvec2(FlmV_recover_pre(V, W, p, pi, P[1]), D);
4712 : }
4713 :
4714 : static int
4715 655 : ZabM_ker_check(GEN M, GEN H, ulong p, GEN P, long n)
4716 : {
4717 : GEN pow;
4718 655 : long j, l = lg(H);
4719 : ulong pi, r;
4720 3899 : do p += n; while(!uisprime(p));
4721 655 : pi = get_Fl_red(p);
4722 655 : P = ZX_to_Flx(P, p);
4723 655 : r = Flx_oneroot(P, p);
4724 655 : pow = Fl_powers_pre(r, degpol(P),p,pi);
4725 655 : M = ZXM_to_FlxM(M, p, P[1]); M = FlxM_eval_powers_pre(M, pow, p, pi);
4726 655 : H = ZXM_to_FlxM(H, p, P[1]); H = FlxM_eval_powers_pre(H, pow, p, pi);
4727 2178 : for (j = 1; j < l; j++)
4728 1555 : if (!zv_equal0(Flm_Flc_mul_pre(M, gel(H,j), p, pi))) return 0;
4729 623 : return 1;
4730 : }
4731 :
4732 : GEN
4733 623 : ZabM_ker(GEN M, GEN P, long n)
4734 : {
4735 623 : pari_sp av = avma;
4736 : pari_timer ti;
4737 623 : GEN q, H = NULL, D = NULL;
4738 623 : ulong m = LONG_MAX>>1;
4739 623 : ulong p = 1 + m - (m % n);
4740 :
4741 623 : if (DEBUGLEVEL>5) timer_start(&ti);
4742 : for(;;)
4743 653 : {
4744 : GEN Kp, Hp, Dp, Pp, Mp, Hr;
4745 22341 : do p += n; while(!uisprime(p));
4746 1276 : Pp = ZX_to_Flx(P, p);
4747 1276 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4748 1276 : Kp = FlkM_ker(Mp, Pp, p);
4749 1276 : if (!Kp) continue;
4750 1276 : Hp = gel(Kp,1); Dp = gel(Kp,2);
4751 1276 : if (H && (lg(Hp)>lg(H) || (lg(Hp)==lg(H) && vecsmall_lexcmp(Dp,D)>0))) continue;
4752 1276 : if (!H || (lg(Hp)<lg(H) || vecsmall_lexcmp(Dp,D)<0))
4753 : {
4754 623 : H = ZXM_init_CRT(Hp, degpol(P)-1, p); D = Dp;
4755 623 : q = utoipos(p);
4756 : }
4757 : else
4758 653 : ZXM_incremental_CRT(&H, Hp, &q, p);
4759 1276 : Hr = FpXM_ratlift(H, q);
4760 1276 : if (DEBUGLEVEL>5) timer_printf(&ti,"ZabM_ker mod %ld (ratlift=%ld)", p,!!Hr);
4761 1276 : if (Hr) {/* DONE ? */
4762 655 : GEN Hl = vec_Q_primpart(Hr);
4763 655 : if (ZabM_ker_check(M, Hl, p, P, n)) { H = Hl; break; }
4764 : }
4765 :
4766 653 : if (gc_needed(av,2))
4767 : {
4768 23 : if (DEBUGMEM>1) pari_warn(warnmem,"ZabM_ker");
4769 23 : gerepileall(av, 3, &H, &D, &q);
4770 : }
4771 : }
4772 623 : return gerepilecopy(av, H);
4773 : }
4774 :
4775 : GEN
4776 2352 : ZabM_indexrank(GEN M, GEN P, long n)
4777 : {
4778 2352 : pari_sp av = avma;
4779 2352 : ulong m = LONG_MAX>>1;
4780 2352 : ulong p = 1+m-(m%n), D = degpol(P);
4781 2352 : long lM = lg(M), lmax = 0, c = 0;
4782 : GEN v;
4783 : for(;;)
4784 735 : {
4785 : GEN R, Pp, Mp, K;
4786 : ulong pi;
4787 : long l;
4788 60720 : do p += n; while (!uisprime(p));
4789 3087 : pi = get_Fl_red(p);
4790 3087 : Pp = ZX_to_Flx(P, p);
4791 3087 : R = Flx_roots(Pp, p);
4792 3087 : Mp = ZXM_to_FlxM(M, p, get_Flx_var(Pp));
4793 3087 : K = FlxM_eval_powers_pre(Mp, Fl_powers_pre(uel(R,1), D,p,pi), p,pi);
4794 3087 : v = Flm_indexrank(K, p);
4795 3087 : l = lg(gel(v,2));
4796 3087 : if (l == lM) break;
4797 980 : if (lmax >= 0 && l > lmax) { lmax = l; c = 0; } else c++;
4798 980 : if (c > 2)
4799 : { /* probably not maximal rank, expensive check */
4800 245 : lM -= lg(ZabM_ker(M, P, n))-1; /* actual rank (+1) */
4801 245 : if (lmax == lM) break;
4802 0 : lmax = -1; /* disable check */
4803 : }
4804 : }
4805 2352 : return gerepileupto(av, v);
4806 : }
4807 :
4808 : #if 0
4809 : GEN
4810 : ZabM_gauss(GEN M, GEN P, long n, GEN *den)
4811 : {
4812 : pari_sp av = avma;
4813 : GEN v, S, W;
4814 : v = ZabM_indexrank(M, P, n);
4815 : S = shallowmatextract(M,gel(v,1),gel(v,2));
4816 : W = ZabM_inv(S, P, n, den);
4817 : gerepileall(av,2,&W,den);
4818 : return W;
4819 : }
4820 : #endif
4821 :
4822 : GEN
4823 140 : ZabM_pseudoinv(GEN M, GEN P, long n, GEN *pv, GEN *den)
4824 : {
4825 140 : GEN v = ZabM_indexrank(M, P, n);
4826 140 : if (pv) *pv = v;
4827 140 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4828 140 : return ZabM_inv(M, P, n, den);
4829 : }
4830 : GEN
4831 4844 : ZM_pseudoinv(GEN M, GEN *pv, GEN *den)
4832 : {
4833 4844 : GEN v = ZM_indexrank(M);
4834 4844 : if (pv) *pv = v;
4835 4844 : M = shallowmatextract(M,gel(v,1),gel(v,2));
4836 4844 : return ZM_inv(M, den);
4837 : }
4838 :
4839 : /*******************************************************************/
4840 : /* */
4841 : /* Structured Elimination */
4842 : /* */
4843 : /*******************************************************************/
4844 :
4845 : static void
4846 108051 : rem_col(GEN c, long i, GEN iscol, GEN Wrow, long *rcol, long *rrow)
4847 : {
4848 108051 : long lc = lg(c), k;
4849 108051 : iscol[i] = 0; (*rcol)--;
4850 1010133 : for (k = 1; k < lc; ++k)
4851 : {
4852 902082 : Wrow[c[k]]--;
4853 902082 : if (Wrow[c[k]]==0) (*rrow)--;
4854 : }
4855 108051 : }
4856 :
4857 : static void
4858 6261 : rem_singleton(GEN M, GEN iscol, GEN Wrow, long idx, long *rcol, long *rrow)
4859 : {
4860 : long i, j;
4861 6261 : long nbcol = lg(iscol)-1, last;
4862 : do
4863 : {
4864 8293 : last = 0;
4865 19310206 : for (i = 1; i <= nbcol; ++i)
4866 19301913 : if (iscol[i])
4867 : {
4868 10012830 : GEN c = idx ? gmael(M, i, idx): gel(M,i);
4869 10012830 : long lc = lg(c);
4870 95090091 : for (j = 1; j < lc; ++j)
4871 85095956 : if (Wrow[c[j]] == 1)
4872 : {
4873 18695 : rem_col(c, i, iscol, Wrow, rcol, rrow);
4874 18695 : last=1; break;
4875 : }
4876 : }
4877 8293 : } while (last);
4878 6261 : }
4879 :
4880 : static GEN
4881 6101 : fill_wcol(GEN M, GEN iscol, GEN Wrow, long *w, GEN wcol)
4882 : {
4883 6101 : long nbcol = lg(iscol)-1;
4884 : long i, j, m, last;
4885 : GEN per;
4886 15290 : for (m = 2, last=0; !last ; m++)
4887 : {
4888 22446705 : for (i = 1; i <= nbcol; ++i)
4889 : {
4890 22437516 : wcol[i] = 0;
4891 22437516 : if (iscol[i])
4892 : {
4893 11520312 : GEN c = gmael(M, i, 1);
4894 11520312 : long lc = lg(c);
4895 107248231 : for (j = 1; j < lc; ++j)
4896 95727919 : if (Wrow[c[j]] == m) { wcol[i]++; last = 1; }
4897 : }
4898 : }
4899 : }
4900 6101 : per = vecsmall_indexsort(wcol);
4901 6101 : *w = wcol[per[nbcol]];
4902 6101 : return per;
4903 : }
4904 :
4905 : /* M is a RgMs with nbrow rows, A a list of row indices.
4906 : Eliminate rows of M with a single entry that do not belong to A,
4907 : and the corresponding columns. Also eliminate columns until #colums=#rows.
4908 : Return pcol and prow:
4909 : pcol is a map from the new columns indices to the old one.
4910 : prow is a map from the old rows indices to the new one (0 if removed).
4911 : */
4912 :
4913 : void
4914 119 : RgMs_structelim_col(GEN M, long nbcol, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4915 : {
4916 : long i,j,k;
4917 119 : long lA = lg(A);
4918 119 : GEN prow = cgetg(nbrow+1, t_VECSMALL);
4919 119 : GEN pcol = zero_zv(nbcol);
4920 119 : pari_sp av = avma;
4921 119 : long rcol = nbcol, rrow = 0, imin = nbcol - usqrt(nbcol);
4922 119 : GEN iscol = const_vecsmall(nbcol, 1);
4923 119 : GEN Wrow = zero_zv(nbrow);
4924 119 : GEN wcol = cgetg(nbcol+1, t_VECSMALL);
4925 119 : pari_sp av2=avma;
4926 126763 : for (i = 1; i <= nbcol; ++i)
4927 : {
4928 126644 : GEN F = gmael(M, i, 1);
4929 126644 : long l = lg(F)-1;
4930 1115940 : for (j = 1; j <= l; ++j)
4931 989296 : Wrow[F[j]]++;
4932 : }
4933 119 : for (j = 1; j < lA; ++j)
4934 : {
4935 0 : if (Wrow[A[j]] == 0) { *p_col=NULL; return; }
4936 0 : Wrow[A[j]] = -1;
4937 : }
4938 235298 : for (i = 1; i <= nbrow; ++i)
4939 235179 : if (Wrow[i])
4940 66651 : rrow++;
4941 119 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow);
4942 119 : if (rcol<rrow) pari_err_BUG("RgMs_structelim, rcol<rrow");
4943 6220 : for (; rcol>rrow;)
4944 : {
4945 : long w;
4946 6101 : GEN per = fill_wcol(M, iscol, Wrow, &w, wcol);
4947 95457 : for (i = nbcol; i>=imin && wcol[per[i]]>=w && rcol>rrow; i--)
4948 89356 : rem_col(gmael(M, per[i], 1), per[i], iscol, Wrow, &rcol, &rrow);
4949 6101 : rem_singleton(M, iscol, Wrow, 1, &rcol, &rrow);
4950 6101 : set_avma(av2);
4951 : }
4952 126763 : for (j = 1, i = 1; i <= nbcol; ++i)
4953 126644 : if (iscol[i])
4954 24929 : pcol[j++] = i;
4955 119 : setlg(pcol,j);
4956 235298 : for (k = 1, i = 1; i <= nbrow; ++i)
4957 235179 : prow[i] = Wrow[i] ? k++: 0;
4958 119 : set_avma(av);
4959 119 : *p_col = pcol; *p_row = prow;
4960 : }
4961 :
4962 : void
4963 0 : RgMs_structelim(GEN M, long nbrow, GEN A, GEN *p_col, GEN *p_row)
4964 : {
4965 0 : RgMs_structelim_col(M, lg(M)-1, nbrow, A, p_col, p_row);
4966 0 : }
4967 :
4968 : GEN
4969 41 : F2Ms_colelim(GEN M, long nbrow)
4970 : {
4971 : long i,j;
4972 41 : long nbcol = lg(M)-1;
4973 41 : GEN pcol = zero_zv(nbcol);
4974 41 : pari_sp av = avma;
4975 41 : long rcol = nbcol, rrow = 0;
4976 41 : GEN iscol = const_vecsmall(nbcol, 1);
4977 41 : GEN Wrow = zero_zv(nbrow);
4978 54612 : for (i = 1; i <= nbcol; ++i)
4979 : {
4980 54571 : GEN F = gel(M, i);
4981 54571 : long l = lg(F)-1;
4982 984416 : for (j = 1; j <= l; ++j)
4983 929845 : Wrow[F[j]]++;
4984 : }
4985 41 : rem_singleton(M, iscol, Wrow, 0, &rcol, &rrow);
4986 54612 : for (j = 1, i = 1; i <= nbcol; ++i)
4987 54571 : if (iscol[i])
4988 48235 : pcol[j++] = i;
4989 41 : fixlg(pcol,j);
4990 41 : set_avma(av);
4991 41 : return pcol;
4992 : }
4993 :
4994 : /*******************************************************************/
4995 : /* */
4996 : /* EIGENVECTORS */
4997 : /* (independent eigenvectors, sorted by increasing eigenvalue) */
4998 : /* */
4999 : /*******************************************************************/
5000 : /* assume x is square of dimension > 0 */
5001 : static int
5002 34 : RgM_is_symmetric_cx(GEN x, long bit)
5003 : {
5004 34 : pari_sp av = avma;
5005 34 : long i, j, l = lg(x);
5006 201 : for (i = 1; i < l; i++)
5007 670 : for (j = 1; j < i; j++)
5008 : {
5009 503 : GEN a = gcoeff(x,i,j), b = gcoeff(x,j,i), c = gsub(a,b);
5010 503 : if (!gequal0(c) && gexpo(c) - gexpo(a) > -bit) return gc_long(av,0);
5011 : }
5012 21 : return gc_long(av,1);
5013 : }
5014 : static GEN
5015 34 : eigen_err(int exact, GEN x, long flag, long prec)
5016 : {
5017 34 : pari_sp av = avma;
5018 34 : if (RgM_is_symmetric_cx(x, prec2nbits(prec) - 10))
5019 : { /* approximately symmetric: recover */
5020 21 : x = jacobi(x, prec); if (flag) return x;
5021 14 : return gerepilecopy(av, gel(x,2));
5022 : }
5023 13 : if (exact)
5024 : {
5025 6 : GEN y = mateigen(x, flag, precdbl(prec));
5026 6 : return gerepilecopy(av, gprec_wtrunc(y, prec));
5027 : }
5028 7 : pari_err_PREC("mateigen");
5029 : return NULL; /* LCOV_EXCL_LINE */
5030 : }
5031 : GEN
5032 104 : mateigen(GEN x, long flag, long prec)
5033 : {
5034 : GEN y, R, T;
5035 104 : long k, l, ex, n = lg(x);
5036 : int exact;
5037 104 : pari_sp av = avma;
5038 :
5039 104 : if (typ(x)!=t_MAT) pari_err_TYPE("eigen",x);
5040 104 : if (n != 1 && n != lgcols(x)) pari_err_DIM("eigen");
5041 104 : if (flag < 0 || flag > 1) pari_err_FLAG("mateigen");
5042 104 : if (n == 1)
5043 : {
5044 14 : if (flag) retmkvec2(cgetg(1,t_VEC), cgetg(1,t_MAT));
5045 7 : return cgetg(1,t_VEC);
5046 : }
5047 90 : if (n == 2)
5048 : {
5049 14 : if (flag) retmkvec2(mkveccopy(gcoeff(x,1,1)), matid(1));
5050 7 : return matid(1);
5051 : }
5052 :
5053 76 : ex = 16 - prec2nbits(prec);
5054 76 : T = charpoly(x,0);
5055 76 : exact = RgX_is_QX(T);
5056 76 : if (exact)
5057 : {
5058 41 : T = ZX_radical( Q_primpart(T) );
5059 41 : R = nfrootsQ(T);
5060 41 : if (lg(R)-1 < degpol(T))
5061 : { /* add missing complex roots */
5062 27 : GEN r = cleanroots(RgX_div(T, roots_to_pol(R, 0)), prec);
5063 27 : settyp(r, t_VEC);
5064 27 : R = shallowconcat(R, r);
5065 : }
5066 : }
5067 : else
5068 : {
5069 35 : GEN r1, v = vectrunc_init(lg(T));
5070 : long e;
5071 35 : R = cleanroots(T,prec);
5072 35 : r1 = NULL;
5073 231 : for (k = 1; k < lg(R); k++)
5074 : {
5075 196 : GEN r2 = gel(R,k), r = grndtoi(r2, &e);
5076 196 : if (e < ex) r2 = r;
5077 196 : if (r1)
5078 : {
5079 161 : r = gsub(r1,r2);
5080 161 : if (gequal0(r) || gexpo(r) < ex) continue;
5081 : }
5082 154 : vectrunc_append(v, r2);
5083 154 : r1 = r2;
5084 : }
5085 35 : R = v;
5086 : }
5087 : /* R = distinct complex roots of charpoly(x) */
5088 76 : l = lg(R); y = cgetg(l, t_VEC);
5089 328 : for (k = 1; k < l; k++)
5090 : {
5091 286 : GEN F = ker_aux(RgM_Rg_sub_shallow(x, gel(R,k)), x);
5092 286 : long d = lg(F)-1;
5093 286 : if (!d) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5094 252 : gel(y,k) = F;
5095 252 : if (flag) gel(R,k) = const_vec(d, gel(R,k));
5096 : }
5097 42 : y = shallowconcat1(y);
5098 42 : if (lg(y) > n) { set_avma(av); return eigen_err(exact, x, flag, prec); }
5099 : /* lg(y) < n if x is not diagonalizable */
5100 42 : if (flag) y = mkvec2(shallowconcat1(R), y);
5101 42 : return gerepilecopy(av,y);
5102 : }
5103 : GEN
5104 0 : eigen(GEN x, long prec) { return mateigen(x, 0, prec); }
5105 :
5106 : /*******************************************************************/
5107 : /* */
5108 : /* DETERMINANT */
5109 : /* */
5110 : /*******************************************************************/
5111 :
5112 : GEN
5113 4074 : det0(GEN a,long flag)
5114 : {
5115 4074 : switch(flag)
5116 : {
5117 4060 : case 0: return det(a);
5118 14 : case 1: return det2(a);
5119 0 : default: pari_err_FLAG("matdet");
5120 : }
5121 : return NULL; /* LCOV_EXCL_LINE */
5122 : }
5123 :
5124 : /* M a 2x2 matrix, returns det(M) */
5125 : static GEN
5126 51977 : RgM_det2(GEN M)
5127 : {
5128 51977 : pari_sp av = avma;
5129 51977 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5130 51977 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5131 51977 : return gerepileupto(av, gsub(gmul(a,d), gmul(b,c)));
5132 : }
5133 : /* M a 2x2 ZM, returns det(M) */
5134 : static GEN
5135 6118 : ZM_det2(GEN M)
5136 : {
5137 6118 : pari_sp av = avma;
5138 6118 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2);
5139 6118 : GEN c = gcoeff(M,2,1), d = gcoeff(M,2,2);
5140 6118 : return gerepileuptoint(av, subii(mulii(a,d), mulii(b, c)));
5141 : }
5142 : /* M a 3x3 ZM, return det(M) */
5143 : static GEN
5144 39326 : ZM_det3(GEN M)
5145 : {
5146 39326 : pari_sp av = avma;
5147 39326 : GEN a = gcoeff(M,1,1), b = gcoeff(M,1,2), c = gcoeff(M,1,3);
5148 39326 : GEN d = gcoeff(M,2,1), e = gcoeff(M,2,2), f = gcoeff(M,2,3);
5149 39326 : GEN g = gcoeff(M,3,1), h = gcoeff(M,3,2), i = gcoeff(M,3,3);
5150 39326 : GEN t, D = signe(i)? mulii(subii(mulii(a,e), mulii(b,d)), i): gen_0;
5151 39326 : if (signe(g))
5152 : {
5153 24941 : t = mulii(subii(mulii(b,f), mulii(c,e)), g);
5154 24941 : D = addii(D, t);
5155 : }
5156 39326 : if (signe(h))
5157 : {
5158 30317 : t = mulii(subii(mulii(c,d), mulii(a,f)), h);
5159 30317 : D = addii(D, t);
5160 : }
5161 39326 : return gerepileuptoint(av, D);
5162 : }
5163 :
5164 : static GEN
5165 56091 : det_simple_gauss(GEN a, GEN data, pivot_fun pivot)
5166 : {
5167 56091 : pari_sp av = avma;
5168 56091 : long i,j,k, s = 1, nbco = lg(a)-1;
5169 56091 : GEN p, x = gen_1;
5170 :
5171 56091 : a = RgM_shallowcopy(a);
5172 380975 : for (i=1; i<nbco; i++)
5173 : {
5174 324892 : k = pivot(a, data, i, NULL);
5175 324894 : if (k > nbco) return gerepilecopy(av, gcoeff(a,i,i));
5176 324887 : if (k != i)
5177 : { /* exchange the lines s.t. k = i */
5178 1498956 : for (j=i; j<=nbco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
5179 144507 : s = -s;
5180 : }
5181 324887 : p = gcoeff(a,i,i);
5182 :
5183 324887 : x = gmul(x,p);
5184 2289877 : for (k=i+1; k<=nbco; k++)
5185 : {
5186 1964991 : GEN m = gcoeff(a,i,k);
5187 1964991 : if (gequal0(m)) continue;
5188 :
5189 1363432 : m = gdiv(m,p);
5190 13398712 : for (j=i+1; j<=nbco; j++)
5191 12035283 : gcoeff(a,j,k) = gsub(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
5192 : }
5193 324886 : if (gc_needed(av,2))
5194 : {
5195 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5196 0 : gerepileall(av,2, &a,&x);
5197 : }
5198 : }
5199 56083 : if (s < 0) x = gneg_i(x);
5200 56083 : return gerepileupto(av, gmul(x, gcoeff(a,nbco,nbco)));
5201 : }
5202 :
5203 : GEN
5204 115282 : det2(GEN a)
5205 : {
5206 : GEN data;
5207 : pivot_fun pivot;
5208 115282 : long n = lg(a)-1;
5209 115282 : if (typ(a)!=t_MAT) pari_err_TYPE("det2",a);
5210 115282 : if (!n) return gen_1;
5211 115282 : if (n != nbrows(a)) pari_err_DIM("det2");
5212 115282 : if (n == 1) return gcopy(gcoeff(a,1,1));
5213 76904 : if (n == 2) return RgM_det2(a);
5214 27741 : pivot = get_pivot_fun(a, a, &data);
5215 27741 : return det_simple_gauss(a, data, pivot);
5216 : }
5217 :
5218 : /* Assumes a a square t_MAT of dimension n > 0. Returns det(a) using
5219 : * Gauss-Bareiss. */
5220 : static GEN
5221 448 : det_bareiss(GEN a)
5222 : {
5223 448 : pari_sp av = avma;
5224 448 : long nbco = lg(a)-1,i,j,k,s = 1;
5225 : GEN p, pprec;
5226 :
5227 448 : a = RgM_shallowcopy(a);
5228 1274 : for (pprec=gen_1,i=1; i<nbco; i++,pprec=p)
5229 : {
5230 826 : int diveuc = (gequal1(pprec)==0);
5231 : GEN ci;
5232 :
5233 826 : p = gcoeff(a,i,i);
5234 826 : if (gequal0(p))
5235 : {
5236 0 : k=i+1; while (k<=nbco && gequal0(gcoeff(a,i,k))) k++;
5237 0 : if (k>nbco) return gerepilecopy(av, p);
5238 0 : swap(gel(a,k), gel(a,i)); s = -s;
5239 0 : p = gcoeff(a,i,i);
5240 : }
5241 826 : ci = gel(a,i);
5242 2184 : for (k=i+1; k<=nbco; k++)
5243 : {
5244 1358 : GEN ck = gel(a,k), m = gel(ck,i);
5245 1358 : if (gequal0(m))
5246 : {
5247 0 : if (gequal1(p))
5248 : {
5249 0 : if (diveuc)
5250 0 : gel(a,k) = gdiv(gel(a,k), pprec);
5251 : }
5252 : else
5253 0 : for (j=i+1; j<=nbco; j++)
5254 : {
5255 0 : GEN p1 = gmul(p, gel(ck,j));
5256 0 : if (diveuc) p1 = gdiv(p1,pprec);
5257 0 : gel(ck,j) = p1;
5258 : }
5259 : }
5260 : else
5261 4088 : for (j=i+1; j<=nbco; j++)
5262 : {
5263 2730 : pari_sp av2 = avma;
5264 2730 : GEN p1 = gsub(gmul(p,gel(ck,j)), gmul(m,gel(ci,j)));
5265 2730 : if (diveuc) p1 = gdiv(p1,pprec);
5266 2730 : gel(ck,j) = gerepileupto(av2, p1);
5267 : }
5268 1358 : if (gc_needed(av,2))
5269 : {
5270 0 : if(DEBUGMEM>1) pari_warn(warnmem,"det. col = %ld",i);
5271 0 : gerepileall(av,2, &a,&pprec);
5272 0 : ci = gel(a,i);
5273 0 : p = gcoeff(a,i,i);
5274 : }
5275 : }
5276 : }
5277 448 : p = gcoeff(a,nbco,nbco);
5278 448 : p = (s < 0)? gneg(p): gcopy(p);
5279 448 : return gerepileupto(av, p);
5280 : }
5281 :
5282 : /* count nonzero entries in col j, at most 'max' of them.
5283 : * Return their indices */
5284 : static GEN
5285 1400 : col_count_non_zero(GEN a, long j, long max)
5286 : {
5287 1400 : GEN v = cgetg(max+1, t_VECSMALL);
5288 1400 : GEN c = gel(a,j);
5289 1400 : long i, l = lg(a), k = 1;
5290 5278 : for (i = 1; i < l; i++)
5291 5040 : if (!gequal0(gel(c,i)))
5292 : {
5293 4788 : if (k > max) return NULL; /* fail */
5294 3626 : v[k++] = i;
5295 : }
5296 238 : setlg(v, k); return v;
5297 : }
5298 : /* count nonzero entries in row i, at most 'max' of them.
5299 : * Return their indices */
5300 : static GEN
5301 1386 : row_count_non_zero(GEN a, long i, long max)
5302 : {
5303 1386 : GEN v = cgetg(max+1, t_VECSMALL);
5304 1386 : long j, l = lg(a), k = 1;
5305 5222 : for (j = 1; j < l; j++)
5306 4998 : if (!gequal0(gcoeff(a,i,j)))
5307 : {
5308 4774 : if (k > max) return NULL; /* fail */
5309 3612 : v[k++] = j;
5310 : }
5311 224 : setlg(v, k); return v;
5312 : }
5313 :
5314 : static GEN det_develop(GEN a, long max, double bound);
5315 : /* (-1)^(i+j) a[i,j] * det RgM_minor(a,i,j) */
5316 : static GEN
5317 406 : coeff_det(GEN a, long i, long j, long max, double bound)
5318 : {
5319 406 : GEN c = gcoeff(a, i, j);
5320 406 : c = gmul(c, det_develop(RgM_minor(a, i,j), max, bound));
5321 406 : if (odd(i+j)) c = gneg(c);
5322 406 : return c;
5323 : }
5324 : /* a square t_MAT, 'bound' a rough upper bound for the number of
5325 : * multiplications we are willing to pay while developing rows/columns before
5326 : * switching to Gaussian elimination */
5327 : static GEN
5328 644 : det_develop(GEN M, long max, double bound)
5329 : {
5330 644 : pari_sp av = avma;
5331 644 : long i,j, n = lg(M)-1, lbest = max+2, best_col = 0, best_row = 0;
5332 644 : GEN best = NULL;
5333 :
5334 644 : if (bound < 1.) return det_bareiss(M); /* too costly now */
5335 :
5336 420 : switch(n)
5337 : {
5338 0 : case 0: return gen_1;
5339 0 : case 1: return gcopy(gcoeff(M,1,1));
5340 14 : case 2: return RgM_det2(M);
5341 : }
5342 406 : if (max > ((n+2)>>1)) max = (n+2)>>1;
5343 1792 : for (j = 1; j <= n; j++)
5344 : {
5345 1400 : pari_sp av2 = avma;
5346 1400 : GEN v = col_count_non_zero(M, j, max);
5347 : long lv;
5348 1400 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5349 182 : if (lv == 1) { set_avma(av); return gen_0; }
5350 182 : if (lv == 2) {
5351 14 : set_avma(av);
5352 14 : return gerepileupto(av, coeff_det(M,v[1],j,max,bound));
5353 : }
5354 168 : best = v; lbest = lv; best_col = j;
5355 : }
5356 1778 : for (i = 1; i <= n; i++)
5357 : {
5358 1386 : pari_sp av2 = avma;
5359 1386 : GEN v = row_count_non_zero(M, i, max);
5360 : long lv;
5361 1386 : if (!v || (lv = lg(v)) >= lbest) { set_avma(av2); continue; }
5362 0 : if (lv == 1) { set_avma(av); return gen_0; }
5363 0 : if (lv == 2) {
5364 0 : set_avma(av);
5365 0 : return gerepileupto(av, coeff_det(M,i,v[1],max,bound));
5366 : }
5367 0 : best = v; lbest = lv; best_row = i;
5368 : }
5369 392 : if (best_row)
5370 : {
5371 0 : double d = lbest-1;
5372 0 : GEN s = NULL;
5373 : long k;
5374 0 : bound /= d*d*d;
5375 0 : for (k = 1; k < lbest; k++)
5376 : {
5377 0 : GEN c = coeff_det(M, best_row, best[k], max, bound);
5378 0 : s = s? gadd(s, c): c;
5379 : }
5380 0 : return gerepileupto(av, s);
5381 : }
5382 392 : if (best_col)
5383 : {
5384 168 : double d = lbest-1;
5385 168 : GEN s = NULL;
5386 : long k;
5387 168 : bound /= d*d*d;
5388 560 : for (k = 1; k < lbest; k++)
5389 : {
5390 392 : GEN c = coeff_det(M, best[k], best_col, max, bound);
5391 392 : s = s? gadd(s, c): c;
5392 : }
5393 168 : return gerepileupto(av, s);
5394 : }
5395 224 : return det_bareiss(M);
5396 : }
5397 :
5398 : /* area of parallelogram bounded by (v1,v2) */
5399 : static GEN
5400 61411 : parallelogramarea(GEN v1, GEN v2)
5401 61411 : { return gsub(gmul(gnorml2(v1), gnorml2(v2)), gsqr(RgV_dotproduct(v1, v2))); }
5402 :
5403 : /* Square of Hadamard bound for det(a), a square matrix.
5404 : * Slight improvement: instead of using the column norms, use the area of
5405 : * the parallelogram formed by pairs of consecutive vectors */
5406 : GEN
5407 19215 : RgM_Hadamard(GEN a)
5408 : {
5409 19215 : pari_sp av = avma;
5410 19215 : long n = lg(a)-1, i;
5411 : GEN B;
5412 19215 : if (n == 0) return gen_1;
5413 19215 : if (n == 1) return gsqr(gcoeff(a,1,1));
5414 19215 : a = RgM_gtofp(a, LOWDEFAULTPREC);
5415 19215 : B = gen_1;
5416 80626 : for (i = 1; i <= n/2; i++)
5417 61411 : B = gmul(B, parallelogramarea(gel(a,2*i-1), gel(a,2*i)));
5418 19215 : if (odd(n)) B = gmul(B, gnorml2(gel(a, n)));
5419 19215 : return gerepileuptoint(av, ceil_safe(B));
5420 : }
5421 :
5422 : /* If B=NULL, assume B=A' */
5423 : static GEN
5424 20563 : ZM_det_slice(GEN A, GEN P, GEN *mod)
5425 : {
5426 20563 : pari_sp av = avma;
5427 20563 : long i, n = lg(P)-1;
5428 : GEN H, T;
5429 20563 : if (n == 1)
5430 : {
5431 0 : ulong Hp, p = uel(P,1);
5432 0 : GEN a = ZM_to_Flm(A, p);
5433 0 : Hp = Flm_det_sp(a, p);
5434 0 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
5435 : }
5436 20563 : T = ZV_producttree(P);
5437 20563 : A = ZM_nv_mod_tree(A, P, T);
5438 20563 : H = cgetg(n+1, t_VECSMALL);
5439 85543 : for(i=1; i <= n; i++)
5440 : {
5441 64980 : ulong p = P[i];
5442 64980 : GEN a = gel(A,i);
5443 64980 : H[i] = Flm_det_sp(a, p);
5444 : }
5445 20563 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
5446 20562 : *mod = gmael(T, lg(T)-1, 1);
5447 20562 : gerepileall(av, 2, &H, mod); return H;
5448 : }
5449 :
5450 : GEN
5451 20563 : ZM_det_worker(GEN P, GEN A)
5452 : {
5453 20563 : GEN V = cgetg(3, t_VEC);
5454 20563 : gel(V,1) = ZM_det_slice(A, P, &gel(V,2));
5455 20563 : return V;
5456 : }
5457 :
5458 : GEN
5459 65975 : ZM_det(GEN M)
5460 : {
5461 65975 : const long DIXON_THRESHOLD = 40;
5462 : pari_sp av, av2;
5463 65975 : long i, n = lg(M)-1;
5464 : ulong p, Dp;
5465 : forprime_t S;
5466 : pari_timer ti;
5467 : GEN H, D, mod, h, q, v, worker;
5468 : #ifdef LONG_IS_64BIT
5469 56550 : const ulong PMAX = 18446744073709551557UL;
5470 : #else
5471 9425 : const ulong PMAX = 4294967291UL;
5472 : #endif
5473 :
5474 65975 : switch(n)
5475 : {
5476 7 : case 0: return gen_1;
5477 1309 : case 1: return icopy(gcoeff(M,1,1));
5478 6118 : case 2: return ZM_det2(M);
5479 39326 : case 3: return ZM_det3(M);
5480 : }
5481 19215 : if (DEBUGLEVEL>=4) timer_start(&ti);
5482 19215 : av = avma; h = RgM_Hadamard(M); /* |D| <= sqrt(h) */
5483 19215 : if (!signe(h)) { set_avma(av); return gen_0; }
5484 19215 : h = sqrti(h);
5485 19215 : if (lgefint(h) == 3 && (ulong)h[2] <= (PMAX >> 1))
5486 : { /* h < p/2 => direct result */
5487 6571 : p = PMAX;
5488 6571 : Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5489 6571 : set_avma(av);
5490 6571 : if (!Dp) return gen_0;
5491 6571 : return (Dp <= (p>>1))? utoipos(Dp): utoineg(p - Dp);
5492 : }
5493 12644 : q = gen_1; Dp = 1;
5494 12644 : init_modular_big(&S);
5495 12644 : p = 0; /* -Wall */
5496 12644 : while (cmpii(q, h) <= 0 && (p = u_forprime_next(&S)))
5497 : {
5498 12644 : av2 = avma; Dp = Flm_det_sp(ZM_to_Flm(M, p), p);
5499 12644 : set_avma(av2);
5500 12644 : if (Dp) break;
5501 0 : q = muliu(q, p);
5502 : }
5503 12644 : if (!p) pari_err_OVERFLOW("ZM_det [ran out of primes]");
5504 12644 : if (!Dp) { set_avma(av); return gen_0; }
5505 12644 : if (mt_nbthreads() > 1 || n <= DIXON_THRESHOLD)
5506 12644 : D = q; /* never competitive when bound is sharp even with 2 threads */
5507 : else
5508 : {
5509 0 : av2 = avma;
5510 0 : v = cgetg(n+1, t_COL);
5511 0 : gel(v, 1) = gen_1; /* ensure content(v) = 1 */
5512 0 : for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
5513 0 : D = Q_denom(ZM_gauss(M, v));
5514 0 : if (expi(D) < expi(h) >> 1)
5515 : { /* First try unlucky, try once more */
5516 0 : for (i = 2; i <= n; i++) gel(v, i) = stoi(random_Fl(15) - 7);
5517 0 : D = lcmii(D, Q_denom(ZM_gauss(M, v)));
5518 : }
5519 0 : D = gerepileuptoint(av2, D);
5520 0 : if (q != gen_1) D = lcmii(D, q);
5521 : }
5522 12644 : if (DEBUGLEVEL >=4)
5523 0 : timer_printf(&ti,"ZM_det: Dixon %ld/%ld bits",expi(D),expi(h));
5524 : /* determinant is a multiple of D */
5525 12644 : if (is_pm1(D)) D = NULL;
5526 12644 : if (D) h = diviiexact(h, D);
5527 12644 : worker = snm_closure(is_entry("_ZM_det_worker"), mkvec(M));
5528 12644 : H = gen_crt("ZM_det", worker, &S, D, expi(h)+1, 0, &mod,
5529 : ZV_chinese, NULL);
5530 12644 : if (D) H = Fp_div(H, D, mod);
5531 12644 : H = Fp_center(H, mod, shifti(mod,-1));
5532 12644 : if (D) H = mulii(H, D);
5533 12644 : return gerepileuptoint(av, H);
5534 : }
5535 :
5536 : static GEN
5537 1519 : RgM_det_FpM(GEN a, GEN p)
5538 : {
5539 1519 : pari_sp av = avma;
5540 : ulong pp, d;
5541 1519 : a = RgM_Fp_init(a,p,&pp);
5542 1519 : switch(pp)
5543 : {
5544 70 : case 0: return gerepileupto(av, Fp_to_mod(FpM_det(a,p),p)); break;
5545 14 : case 2: d = F2m_det_sp(a); break;
5546 1435 : default:d = Flm_det_sp(a, pp); break;
5547 : }
5548 1449 : set_avma(av); return mkintmodu(d, pp);
5549 : }
5550 :
5551 : static GEN
5552 42 : RgM_det_FqM(GEN x, GEN pol, GEN p)
5553 : {
5554 42 : pari_sp av = avma;
5555 42 : GEN b, T = RgX_to_FpX(pol, p);
5556 42 : if (signe(T) == 0) pari_err_OP("%",x,pol);
5557 42 : b = FqM_det(RgM_to_FqM(x, T, p), T, p);
5558 42 : if (!b) return gc_NULL(av);
5559 42 : return gerepilecopy(av, mkpolmod(FpX_to_mod(b, p), FpX_to_mod(T, p)));
5560 : }
5561 :
5562 : #define code(t1,t2) ((t1 << 6) | t2)
5563 : static GEN
5564 30583 : RgM_det_fast(GEN x)
5565 : {
5566 : GEN p, pol;
5567 : long pa;
5568 30583 : long t = RgM_type(x, &p,&pol,&pa);
5569 30583 : switch(t)
5570 : {
5571 175 : case t_INT: return ZM_det(x);
5572 196 : case t_FRAC: return QM_det(x);
5573 63 : case t_FFELT: return FFM_det(x, pol);
5574 1519 : case t_INTMOD: return RgM_det_FpM(x, p);
5575 42 : case code(t_POLMOD, t_INTMOD):
5576 42 : return RgM_det_FqM(x, pol, p);
5577 28588 : default: return NULL;
5578 : }
5579 : }
5580 : #undef code
5581 :
5582 : static long
5583 238 : det_init_max(long n)
5584 : {
5585 238 : if (n > 100) return 0;
5586 238 : if (n > 50) return 1;
5587 238 : if (n > 30) return 4;
5588 238 : return 7;
5589 : }
5590 :
5591 : GEN
5592 278474 : det(GEN a)
5593 : {
5594 278474 : long n = lg(a)-1;
5595 : double B;
5596 : GEN data, b;
5597 : pivot_fun pivot;
5598 :
5599 278474 : if (typ(a)!=t_MAT) pari_err_TYPE("det",a);
5600 278474 : if (!n) return gen_1;
5601 278432 : if (n != nbrows(a)) pari_err_DIM("det");
5602 278425 : if (n == 1) return gcopy(gcoeff(a,1,1));
5603 33383 : if (n == 2) return RgM_det2(a);
5604 30583 : b = RgM_det_fast(a);
5605 30583 : if (b) return b;
5606 28588 : pivot = get_pivot_fun(a, a, &data);
5607 28588 : if (pivot != gauss_get_pivot_NZ) return det_simple_gauss(a, data, pivot);
5608 238 : B = (double)n;
5609 238 : return det_develop(a, det_init_max(n), B*B*B);
5610 : }
5611 :
5612 : GEN
5613 196 : QM_det(GEN M)
5614 : {
5615 196 : pari_sp av = avma;
5616 196 : GEN cM, pM = Q_primitive_part(M, &cM);
5617 196 : GEN b = ZM_det(pM);
5618 196 : if (cM) b = gmul(b, gpowgs(cM, lg(M)-1));
5619 196 : return gerepileupto(av, b);
5620 : }
|
