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