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