1 /* Copyright (C) 2000-2019 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 #include "pari.h"
16 #include "paripriv.h"
17
18 GEN
Flv_to_ZV(GEN x)19 Flv_to_ZV(GEN x)
20 { pari_APPLY_type(t_VEC, utoi(x[i])) }
21
22 GEN
Flc_to_ZC(GEN x)23 Flc_to_ZC(GEN x)
24 { pari_APPLY_type(t_COL, utoi(x[i])) }
25
26 GEN
Flm_to_ZM(GEN x)27 Flm_to_ZM(GEN x)
28 { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) }
29
30 GEN
Flc_to_ZC_inplace(GEN z)31 Flc_to_ZC_inplace(GEN z)
32 {
33 long i, l = lg(z);
34 for (i=1; i<l; i++) gel(z,i) = utoi(z[i]);
35 settyp(z, t_COL);
36 return z;
37 }
38
39 GEN
Flm_to_ZM_inplace(GEN z)40 Flm_to_ZM_inplace(GEN z)
41 {
42 long i, l = lg(z);
43 for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
44 return z;
45 }
46
47 static GEN
Flm_solve_upper_1(GEN U,GEN B,ulong p,ulong pi)48 Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
49 { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
50
51 static GEN
Flm_rsolve_upper_2(GEN U,GEN B,ulong p,ulong pi)52 Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
53 {
54 ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
55 ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
56 ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
57 ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
58 GEN B1 = rowslice(B, 1, 1);
59 GEN B2 = rowslice(B, 2, 2);
60 GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
61 GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
62 ainv, p, pi);
63 return vconcat(X1, X2);
64 }
65
66 /* solve U*X = B, U upper triangular and invertible */
67 static GEN
Flm_rsolve_upper_pre(GEN U,GEN B,ulong p,ulong pi)68 Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
69 {
70 long n = lg(U) - 1, n1;
71 GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
72 pari_sp av = avma;
73
74 if (n == 0) return B;
75 if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
76 if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
77 n1 = (n + 1)/2;
78 U2 = vecslice(U, n1 + 1, n);
79 U22 = rowslice(U2, n1 + 1, n);
80 B2 = rowslice(B, n1 + 1, n);
81 X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
82 U12 = rowslice(U2, 1, n1);
83 B1 = rowslice(B, 1, n1);
84 B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
85 if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
86 U11 = matslice(U, 1,n1, 1,n1);
87 X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
88 X = vconcat(X1, X2);
89 if (gc_needed(av, 1)) X = gerepilecopy(av, X);
90 return X;
91 }
92
93 static GEN
Flm_lsolve_upper_2(GEN U,GEN B,ulong p,ulong pi)94 Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
95 {
96 ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
97 ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
98 ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
99 ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
100 GEN B1 = vecslice(B, 1, 1);
101 GEN B2 = vecslice(B, 2, 2);
102 GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
103 GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
104 dinv, p, pi);
105 return shallowconcat(X1, X2);
106 }
107
108 /* solve X*U = B, U upper triangular and invertible */
109 static GEN
Flm_lsolve_upper_pre(GEN U,GEN B,ulong p,ulong pi)110 Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
111 {
112 long n = lg(U) - 1, n1;
113 GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
114 pari_sp av = avma;
115
116 if (n == 0) return B;
117 if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
118 if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
119 n1 = (n + 1)/2;
120 U2 = vecslice(U, n1 + 1, n);
121 U11 = matslice(U, 1,n1, 1,n1);
122 U12 = rowslice(U2, 1, n1);
123 U22 = rowslice(U2, n1 + 1, n);
124 B1 = vecslice(B, 1, n1);
125 B2 = vecslice(B, n1 + 1, n);
126 X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
127 B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
128 if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
129 X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
130 X = shallowconcat(X1, X2);
131 if (gc_needed(av, 1)) X = gerepilecopy(av, X);
132 return X;
133 }
134
135 static GEN
Flm_rsolve_lower_unit_2(GEN L,GEN A,ulong p,ulong pi)136 Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
137 {
138 GEN X1 = rowslice(A, 1, 1);
139 GEN X2 = Flm_sub(rowslice(A, 2, 2),
140 Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
141 return vconcat(X1, X2);
142 }
143
144 /* solve L*X = A, L lower triangular with ones on the diagonal
145 * (at least as many rows as columns) */
146 static GEN
Flm_rsolve_lower_unit_pre(GEN L,GEN A,ulong p,ulong pi)147 Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
148 {
149 long m = lg(L) - 1, m1, n;
150 GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
151 pari_sp av = avma;
152
153 if (m == 0) return zero_Flm(0, lg(A) - 1);
154 if (m == 1) return rowslice(A, 1, 1);
155 if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
156 m1 = (m + 1)/2;
157 n = nbrows(L);
158 L1 = vecslice(L, 1, m1);
159 L11 = rowslice(L1, 1, m1);
160 L21 = rowslice(L1, m1 + 1, n);
161 A1 = rowslice(A, 1, m1);
162 X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
163 A2 = rowslice(A, m1 + 1, n);
164 A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
165 if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
166 L22 = matslice(L, m1+1,n, m1+1,m);
167 X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
168 X = vconcat(X1, X2);
169 if (gc_needed(av, 1)) X = gerepilecopy(av, X);
170 return X;
171 }
172
173 static GEN
Flm_lsolve_lower_unit_2(GEN L,GEN A,ulong p,ulong pi)174 Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
175 {
176 GEN X2 = vecslice(A, 2, 2);
177 GEN X1 = Flm_sub(vecslice(A, 1, 1),
178 Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
179 return shallowconcat(X1, X2);
180 }
181
182 /* solve L*X = A, L square lower triangular with ones on the diagonal */
183 static GEN
Flm_lsolve_lower_unit_pre(GEN L,GEN A,ulong p,ulong pi)184 Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
185 {
186 long m = lg(L) - 1, m1;
187 GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
188 pari_sp av = avma;
189
190 if (m <= 1) return A;
191 if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
192 m1 = (m + 1)/2;
193 L2 = vecslice(L, m1 + 1, m);
194 L22 = rowslice(L2, m1 + 1, m);
195 A2 = vecslice(A, m1 + 1, m);
196 X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
197 if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
198 L1 = vecslice(L, 1, m1);
199 L21 = rowslice(L1, m1 + 1, m);
200 A1 = vecslice(A, 1, m1);
201 A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
202 L11 = rowslice(L1, 1, m1);
203 if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
204 X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
205 X = shallowconcat(X1, X2);
206 if (gc_needed(av, 1)) X = gerepilecopy(av, X);
207 return X;
208 }
209
210 /* destroy A */
211 static long
Flm_CUP_basecase(GEN A,GEN * R,GEN * C,GEN * U,GEN * P,ulong p,ulong pi)212 Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
213 {
214 long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
215 ulong u, v;
216
217 if (P) *P = identity_perm(n);
218 *R = cgetg(m + 1, t_VECSMALL);
219 for (j = 1, pr = 0; j <= n; j++)
220 {
221 for (pr++, pc = 0; pr <= m; pr++)
222 {
223 for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
224 if (pc) break;
225 }
226 if (!pc) break;
227 (*R)[j] = pr;
228 if (pc != j)
229 {
230 swap(gel(A, j), gel(A, pc));
231 if (P) lswap((*P)[j], (*P)[pc]);
232 }
233 u = Fl_inv(ucoeff(A, pr, j), p);
234 for (i = pr + 1; i <= m; i++)
235 {
236 v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
237 ucoeff(A, i, j) = v;
238 v = Fl_neg(v, p);
239 for (k = j + 1; k <= n; k++)
240 ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
241 ucoeff(A, pr, k), v, p, pi);
242 }
243 }
244 setlg(*R, j);
245 *C = vecslice(A, 1, j - 1);
246 if (U) *U = rowpermute(A, *R);
247 return j - 1;
248 }
249
250 static const long Flm_CUP_LIMIT = 8;
251
252 static long
Flm_CUP_pre(GEN A,GEN * R,GEN * C,GEN * U,GEN * P,ulong p,ulong pi)253 Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
254 {
255 long m = nbrows(A), m1, n = lg(A) - 1, i, r1, r2, r;
256 GEN R1, C1, U1, P1, R2, C2, U2, P2;
257 GEN A1, A2, B2, C21, U11, U12, T21, T22;
258 pari_sp av = avma;
259
260 if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
261 /* destroy A; not called at the outermost recursion level */
262 return Flm_CUP_basecase(A, R, C, U, P, p, pi);
263 m1 = (minss(m, n) + 1)/2;
264 A1 = rowslice(A, 1, m1);
265 A2 = rowslice(A, m1 + 1, m);
266 r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
267 if (r1 == 0)
268 {
269 r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
270 *R = cgetg(r2 + 1, t_VECSMALL);
271 for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
272 *C = vconcat(zero_Flm(m1, r2), C2);
273 *U = U2;
274 *P = P2;
275 r = r2;
276 }
277 else
278 {
279 U11 = vecslice(U1, 1, r1);
280 U12 = vecslice(U1, r1 + 1, n);
281 T21 = vecslicepermute(A2, P1, 1, r1);
282 T22 = vecslicepermute(A2, P1, r1 + 1, n);
283 C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
284 if (gc_needed(av, 1))
285 gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
286 B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
287 r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
288 r = r1 + r2;
289 *R = cgetg(r + 1, t_VECSMALL);
290 for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
291 for (; i <= r; i++) (*R)[i] = R2[i - r1] + m1;
292 *C = shallowconcat(vconcat(C1, C21),
293 vconcat(zero_Flm(m1, r2), C2));
294 *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
295 vconcat(vecpermute(U12, P2), U2));
296 *P = cgetg(n + 1, t_VECSMALL);
297 for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
298 for (; i <= n; i++) (*P)[i] = P1[P2[i - r1] + r1];
299 }
300 if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
301 return r;
302 }
303
304 static long
Flm_echelon_gauss(GEN A,GEN * R,GEN * C,ulong p,ulong pi)305 Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
306 { return Flm_CUP_basecase(A, R, C, NULL, NULL, p, pi); }
307
308 /* complement of a strictly increasing subsequence of (1, 2, ..., n) */
309 static GEN
indexcompl(GEN v,long n)310 indexcompl(GEN v, long n)
311 {
312 long i, j, k, m = lg(v) - 1;
313 GEN w = cgetg(n - m + 1, t_VECSMALL);
314 for (i = j = k = 1; i <= n; i++)
315 if (j <= m && v[j] == i) j++; else w[k++] = i;
316 return w;
317 }
318
319 /* column echelon form */
320 static long
Flm_echelon_pre(GEN A,GEN * R,GEN * C,ulong p,ulong pi)321 Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
322 {
323 long j, j1, j2, m = nbrows(A), n = lg(A) - 1, n1, r, r1, r2;
324 GEN A1, A2, R1, R1c, C1, R2, C2;
325 GEN A12, A22, B2, C11, C21, M12;
326 pari_sp av = avma;
327
328 if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
329 return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
330
331 n1 = (n + 1)/2;
332 A1 = vecslice(A, 1, n1);
333 A2 = vecslice(A, n1 + 1, n);
334 r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
335 if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
336 if (r1 == m) { *R = R1; *C = C1; return r1; }
337
338 R1c = indexcompl(R1, m);
339 C11 = rowpermute(C1, R1);
340 C21 = rowpermute(C1, R1c);
341 A12 = rowpermute(A2, R1);
342 A22 = rowpermute(A2, R1c);
343 M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
344 B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
345 r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
346 if (!r2) { *R = R1; *C = C1; r = r1; }
347 else
348 {
349 R2 = perm_mul(R1c, R2);
350 C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
351 perm_inv(vecsmall_concat(R1, R1c)));
352 r = r1 + r2;
353 *R = cgetg(r + 1, t_VECSMALL);
354 *C = cgetg(r + 1, t_MAT);
355 for (j = j1 = j2 = 1; j <= r; j++)
356 if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
357 {
358 gel(*C, j) = gel(C1, j1);
359 (*R)[j] = R1[j1++];
360 }
361 else
362 {
363 gel(*C, j) = gel(C2, j2);
364 (*R)[j] = R2[j2++];
365 }
366 }
367 if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
368 return r;
369 }
370
371 static void /* assume m < p */
_Fl_addmul(GEN b,long k,long i,ulong m,ulong p,ulong pi)372 _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
373 {
374 uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
375 }
376 static void /* same m = 1 */
_Fl_add(GEN b,long k,long i,ulong p)377 _Fl_add(GEN b, long k, long i, ulong p)
378 {
379 uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
380 }
381 static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
_Fl_addmul_OK(GEN b,long k,long i,ulong m,ulong p)382 _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
383 {
384 uel(b,k) += m * uel(b,i);
385 if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
386 }
387 static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
_Fl_add_OK(GEN b,long k,long i,ulong p)388 _Fl_add_OK(GEN b, long k, long i, ulong p)
389 {
390 uel(b,k) += uel(b,i);
391 if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
392 }
393
394 /* assume 0 <= a[i,j] < p */
395 static GEN
Fl_get_col_OK(GEN a,GEN b,long li,ulong p)396 Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
397 {
398 GEN u = cgetg(li+1,t_VECSMALL);
399 ulong m = uel(b,li) % p;
400 long i,j;
401
402 uel(u,li) = (m * ucoeff(a,li,li)) % p;
403 for (i = li-1; i > 0; i--)
404 {
405 m = p - uel(b,i)%p;
406 for (j = i+1; j <= li; j++) {
407 if (m & HIGHBIT) m %= p;
408 m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
409 }
410 m %= p;
411 if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
412 uel(u,i) = m;
413 }
414 return u;
415 }
416 static GEN
Fl_get_col(GEN a,GEN b,long li,ulong p)417 Fl_get_col(GEN a, GEN b, long li, ulong p)
418 {
419 GEN u = cgetg(li+1,t_VECSMALL);
420 ulong m = uel(b,li) % p;
421 long i,j;
422
423 uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
424 for (i=li-1; i>0; i--)
425 {
426 m = b[i]%p;
427 for (j = i+1; j <= li; j++)
428 m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
429 if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
430 uel(u,i) = m;
431 }
432 return u;
433 }
434
435 static GEN
Flm_ker_gauss_OK(GEN x,ulong p,long deplin)436 Flm_ker_gauss_OK(GEN x, ulong p, long deplin)
437 {
438 GEN y, c, d;
439 long i, j, k, r, t, m, n;
440 ulong a;
441
442 n = lg(x)-1;
443 m=nbrows(x); r=0;
444
445 c = zero_zv(m);
446 d = cgetg(n+1, t_VECSMALL);
447 a = 0; /* for gcc -Wall */
448 for (k=1; k<=n; k++)
449 {
450 for (j=1; j<=m; j++)
451 if (!c[j])
452 {
453 a = ucoeff(x,j,k) % p;
454 if (a) break;
455 }
456 if (j > m)
457 {
458 if (deplin==1) {
459 c = cgetg(n+1, t_VECSMALL);
460 for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
461 c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
462 return c;
463 }
464 r++; d[k] = 0;
465 }
466 else
467 {
468 ulong piv = p - Fl_inv(a, p); /* -1/a */
469 c[j] = k; d[k] = j;
470 ucoeff(x,j,k) = p-1;
471 if (piv != 1)
472 for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
473 for (t=1; t<=m; t++)
474 {
475 if (t == j) continue;
476
477 piv = ( ucoeff(x,t,k) %= p );
478 if (!piv) continue;
479 if (piv == 1)
480 for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
481 else
482 for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
483 }
484 }
485 }
486 if (deplin==1) return NULL;
487
488 y = cgetg(r+1, t_MAT);
489 for (j=k=1; j<=r; j++,k++)
490 {
491 GEN C = cgetg(n+1, t_VECSMALL);
492
493 gel(y,j) = C; while (d[k]) k++;
494 for (i=1; i<k; i++)
495 if (d[i])
496 uel(C,i) = ucoeff(x,d[i],k) % p;
497 else
498 uel(C,i) = 0UL;
499 uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
500 }
501 if (deplin == 2) {
502 GEN pc = cgetg(n - r + 1, t_VECSMALL); /* indices of pivot columns */
503 for (i = j = 1; j <= n; j++)
504 if (d[j]) pc[i++] = j;
505 return mkvec2(y, pc);
506 }
507 return y;
508 }
509
510 /* in place, destroy x */
511 static GEN
Flm_ker_gauss(GEN x,ulong p,long deplin)512 Flm_ker_gauss(GEN x, ulong p, long deplin)
513 {
514 GEN y, c, d;
515 long i, j, k, r, t, m, n;
516 ulong a, pi;
517 n = lg(x)-1;
518 if (!n) return cgetg(1,t_MAT);
519 if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
520 pi = get_Fl_red(p);
521
522 m=nbrows(x); r=0;
523
524 c = zero_zv(m);
525 d = cgetg(n+1, t_VECSMALL);
526 a = 0; /* for gcc -Wall */
527 for (k=1; k<=n; k++)
528 {
529 for (j=1; j<=m; j++)
530 if (!c[j])
531 {
532 a = ucoeff(x,j,k);
533 if (a) break;
534 }
535 if (j > m)
536 {
537 if (deplin==1) {
538 c = cgetg(n+1, t_VECSMALL);
539 for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
540 c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
541 return c;
542 }
543 r++; d[k] = 0;
544 }
545 else
546 {
547 ulong piv = p - Fl_inv(a, p); /* -1/a */
548 c[j] = k; d[k] = j;
549 ucoeff(x,j,k) = p-1;
550 if (piv != 1)
551 for (i=k+1; i<=n; i++)
552 ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
553 for (t=1; t<=m; t++)
554 {
555 if (t == j) continue;
556
557 piv = ucoeff(x,t,k);
558 if (!piv) continue;
559 if (piv == 1)
560 for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
561 else
562 for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
563 }
564 }
565 }
566 if (deplin==1) return NULL;
567
568 y = cgetg(r+1, t_MAT);
569 for (j=k=1; j<=r; j++,k++)
570 {
571 GEN C = cgetg(n+1, t_VECSMALL);
572
573 gel(y,j) = C; while (d[k]) k++;
574 for (i=1; i<k; i++)
575 if (d[i])
576 uel(C,i) = ucoeff(x,d[i],k);
577 else
578 uel(C,i) = 0UL;
579 uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
580 }
581 if (deplin == 2) {
582 GEN pc = cgetg(n - r + 1, t_VECSMALL); /* indices of pivot columns */
583 for (i = j = 1; j <= n; j++)
584 if (d[j]) pc[i++] = j;
585 return mkvec2(y, pc);
586 }
587 return y;
588 }
589
590 GEN
Flm_intersect(GEN x,GEN y,ulong p)591 Flm_intersect(GEN x, GEN y, ulong p)
592 {
593 pari_sp av = avma;
594 long j, lx = lg(x);
595 GEN z;
596
597 if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
598 z = Flm_ker(shallowconcat(x,y), p);
599 for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
600 return gerepileupto(av, Flm_mul(x,z,p));
601 }
602
603 static GEN
Flm_ker_echelon(GEN x,ulong p,long pivots)604 Flm_ker_echelon(GEN x, ulong p, long pivots) {
605 pari_sp av = avma;
606 GEN R, Rc, C, C1, C2, S, K;
607 long n = lg(x) - 1, r;
608 ulong pi = get_Fl_red(p);
609 r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
610 Rc = indexcompl(R, n);
611 C1 = rowpermute(C, R);
612 C2 = rowpermute(C, Rc);
613 S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
614 K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
615 perm_inv(vecsmall_concat(R, Rc)));
616 K = Flm_transpose(K);
617 if (pivots)
618 return gerepilecopy(av, mkvec2(K, R));
619 return gerepileupto(av, K);
620 }
621
622 static GEN
Flm_deplin_echelon(GEN x,ulong p)623 Flm_deplin_echelon(GEN x, ulong p) {
624 pari_sp av = avma;
625 GEN R, Rc, C, C1, C2, s, v;
626 long i, n = lg(x) - 1, r;
627 ulong pi = get_Fl_red(p);
628 r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
629 if (r == n) return gc_NULL(av);
630 Rc = indexcompl(R, n);
631 i = Rc[1];
632 C1 = rowpermute(C, R);
633 C2 = rowslice(C, i, i);
634 s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
635 v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
636 perm_inv(vecsmall_concat(R, Rc)));
637 return gerepileuptoleaf(av, v);
638 }
639
640 static GEN
Flm_ker_i(GEN x,ulong p,long deplin,long inplace)641 Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
642 if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
643 switch(deplin) {
644 case 0: return Flm_ker_echelon(x, p, 0);
645 case 1: return Flm_deplin_echelon(x, p);
646 case 2: return Flm_ker_echelon(x, p, 1);
647 }
648 return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
649 }
650
651 GEN
Flm_ker_sp(GEN x,ulong p,long deplin)652 Flm_ker_sp(GEN x, ulong p, long deplin) {
653 return Flm_ker_i(x, p, deplin, 1);
654 }
655
656 GEN
Flm_ker(GEN x,ulong p)657 Flm_ker(GEN x, ulong p) {
658 return Flm_ker_i(x, p, 0, 0);
659 }
660
661 GEN
Flm_deplin(GEN x,ulong p)662 Flm_deplin(GEN x, ulong p) {
663 return Flm_ker_i(x, p, 1, 0);
664 }
665
666 /* in place, destroy a, SMALL_ULONG(p) is TRUE */
667 static ulong
Flm_det_gauss_OK(GEN a,long nbco,ulong p)668 Flm_det_gauss_OK(GEN a, long nbco, ulong p)
669 {
670 long i,j,k, s = 1;
671 ulong q, x = 1;
672
673 for (i=1; i<nbco; i++)
674 {
675 for(k=i; k<=nbco; k++)
676 {
677 ulong c = ucoeff(a,k,i) % p;
678 ucoeff(a,k,i) = c;
679 if (c) break;
680 }
681 for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
682 if (k > nbco) return ucoeff(a,i,i);
683 if (k != i)
684 { /* exchange the lines s.t. k = i */
685 for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
686 s = -s;
687 }
688 q = ucoeff(a,i,i);
689
690 if (x & HIGHMASK) x %= p;
691 x *= q;
692 q = Fl_inv(q,p);
693 for (k=i+1; k<=nbco; k++)
694 {
695 ulong m = ucoeff(a,i,k) % p;
696 if (!m) continue;
697
698 m = p - ((m*q)%p);
699 for (j=i+1; j<=nbco; j++)
700 {
701 ulong c = ucoeff(a,j,k);
702 if (c & HIGHMASK) c %= p;
703 ucoeff(a,j,k) = c + m*ucoeff(a,j,i);
704 }
705 }
706 }
707 if (x & HIGHMASK) x %= p;
708 q = ucoeff(a,nbco,nbco);
709 if (q & HIGHMASK) q %= p;
710 x = (x*q) % p;
711 if (s < 0 && x) x = p - x;
712 return x;
713 }
714
715 /* in place, destroy a */
716 static ulong
Flm_det_gauss(GEN a,ulong p)717 Flm_det_gauss(GEN a, ulong p)
718 {
719 long i,j,k, s = 1, nbco = lg(a)-1;
720 ulong pi, q, x = 1;
721
722 if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
723 pi = get_Fl_red(p);
724 for (i=1; i<nbco; i++)
725 {
726 for(k=i; k<=nbco; k++)
727 if (ucoeff(a,k,i)) break;
728 if (k > nbco) return ucoeff(a,i,i);
729 if (k != i)
730 { /* exchange the lines s.t. k = i */
731 for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
732 s = -s;
733 }
734 q = ucoeff(a,i,i);
735
736 x = Fl_mul_pre(x, q, p, pi);
737 q = Fl_inv(q,p);
738 for (k=i+1; k<=nbco; k++)
739 {
740 ulong m = ucoeff(a,i,k);
741 if (!m) continue;
742
743 m = Fl_mul_pre(m, q, p, pi);
744 for (j=i+1; j<=nbco; j++)
745 ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
746 }
747 }
748 if (s < 0) x = Fl_neg(x, p);
749 return Fl_mul(x, ucoeff(a,nbco,nbco), p);
750 }
751
752 static ulong
Flm_det_CUP(GEN a,ulong p)753 Flm_det_CUP(GEN a, ulong p) {
754 GEN R, C, U, P;
755 long i, n = lg(a) - 1, r;
756 ulong d;
757 ulong pi = get_Fl_red(p);
758 r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
759 if (r < n)
760 d = 0;
761 else {
762 d = perm_sign(P) == 1? 1: p-1;
763 for (i = 1; i <= n; i++)
764 d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
765 }
766 return d;
767 }
768
769 static ulong
Flm_det_i(GEN x,ulong p,long inplace)770 Flm_det_i(GEN x, ulong p, long inplace) {
771 pari_sp av = avma;
772 ulong d;
773 if (lg(x) - 1 >= Flm_CUP_LIMIT)
774 d = Flm_det_CUP(x, p);
775 else
776 d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
777 return gc_ulong(av, d);
778 }
779
780 ulong
Flm_det_sp(GEN x,ulong p)781 Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); }
782 ulong
Flm_det(GEN x,ulong p)783 Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); }
784
785 /* Destroy x */
786 static GEN
Flm_gauss_pivot(GEN x,ulong p,long * rr)787 Flm_gauss_pivot(GEN x, ulong p, long *rr)
788 {
789 GEN c,d;
790 long i,j,k,r,t,n,m;
791
792 n=lg(x)-1; if (!n) { *rr=0; return NULL; }
793
794 m=nbrows(x); r=0;
795 d=cgetg(n+1,t_VECSMALL);
796 c = zero_zv(m);
797 for (k=1; k<=n; k++)
798 {
799 for (j=1; j<=m; j++)
800 if (!c[j])
801 {
802 ucoeff(x,j,k) %= p;
803 if (ucoeff(x,j,k)) break;
804 }
805 if (j>m) { r++; d[k]=0; }
806 else
807 {
808 ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
809 c[j]=k; d[k]=j;
810 for (i=k+1; i<=n; i++)
811 ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
812 for (t=1; t<=m; t++)
813 if (!c[t]) /* no pivot on that line yet */
814 {
815 piv = ucoeff(x,t,k);
816 if (piv)
817 {
818 ucoeff(x,t,k) = 0;
819 for (i=k+1; i<=n; i++)
820 ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
821 Fl_mul(piv,ucoeff(x,j,i),p),p);
822 }
823 }
824 for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
825 }
826 }
827 *rr = r; return gc_const((pari_sp)d, d);
828 }
829
830 static GEN
Flm_pivots_CUP(GEN x,ulong p,long * rr)831 Flm_pivots_CUP(GEN x, ulong p, long *rr)
832 {
833 long i, n = lg(x) - 1, r;
834 GEN R, C, U, P, d = zero_zv(n);
835 ulong pi = get_Fl_red(p);
836 r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
837 for(i = 1; i <= r; i++)
838 d[P[i]] = R[i];
839 *rr = n - r; return gc_const((pari_sp)d, d);
840 }
841
842 GEN
Flm_pivots(GEN x,ulong p,long * rr,long inplace)843 Flm_pivots(GEN x, ulong p, long *rr, long inplace)
844 {
845 if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
846 return Flm_pivots_CUP(x, p, rr);
847 return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
848 }
849
850 long
Flm_rank(GEN x,ulong p)851 Flm_rank(GEN x, ulong p)
852 {
853 pari_sp av = avma;
854 long r;
855 if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
856 GEN R, C;
857 ulong pi = get_Fl_red(p);
858 return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
859 }
860 (void) Flm_pivots(x, p, &r, 0);
861 return gc_long(av, lg(x)-1 - r);
862 }
863
864 /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
865 * reduced mod p */
866 static GEN
Flm_inv_upper_1_ind(GEN A,long index,ulong p)867 Flm_inv_upper_1_ind(GEN A, long index, ulong p)
868 {
869 long n = lg(A)-1, i = index, j;
870 GEN u = const_vecsmall(n, 0);
871 u[i] = 1;
872 if (SMALL_ULONG(p))
873 for (i--; i>0; i--)
874 {
875 ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
876 for (j=i+2; j<=n; j++)
877 {
878 if (m & HIGHMASK) m %= p;
879 m += ucoeff(A,i,j) * uel(u,j);
880 }
881 u[i] = Fl_neg(m % p, p);
882 }
883 else
884 for (i--; i>0; i--)
885 {
886 ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
887 for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
888 u[i] = Fl_neg(m, p);
889 }
890 return u;
891 }
892 static GEN
Flm_inv_upper_1(GEN A,ulong p)893 Flm_inv_upper_1(GEN A, ulong p)
894 {
895 long i, l;
896 GEN B = cgetg_copy(A, &l);
897 for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
898 return B;
899 }
900 /* destroy a, b */
901 static GEN
Flm_gauss_sp_OK(GEN a,GEN b,ulong * detp,ulong p)902 Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
903 {
904 long i, j, k, li, bco, aco = lg(a)-1, s = 1;
905 ulong det = 1;
906 GEN u;
907
908 li = nbrows(a);
909 bco = lg(b)-1;
910 for (i=1; i<=aco; i++)
911 {
912 ulong invpiv;
913 /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
914 for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
915 for (k = i; k <= li; k++)
916 {
917 ulong piv = ( ucoeff(a,k,i) %= p );
918 if (piv)
919 {
920 ucoeff(a,k,i) = Fl_inv(piv, p);
921 if (detp)
922 {
923 if (det & HIGHMASK) det %= p;
924 det *= piv;
925 }
926 break;
927 }
928 }
929 /* found a pivot on line k */
930 if (k > li) return NULL;
931 if (k != i)
932 { /* swap lines so that k = i */
933 s = -s;
934 for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
935 for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
936 }
937 if (i == aco) break;
938
939 invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
940 for (k=i+1; k<=li; k++)
941 {
942 ulong m = ( ucoeff(a,k,i) %= p );
943 if (!m) continue;
944
945 m = Fl_mul(m, invpiv, p);
946 if (m == 1) {
947 for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
948 for (j=1; j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
949 } else {
950 for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
951 for (j=1; j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
952 }
953 }
954 }
955 if (detp)
956 {
957 det %= p;
958 if (s < 0 && det) det = p - det;
959 *detp = det;
960 }
961 u = cgetg(bco+1,t_MAT);
962 for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
963 return u;
964 }
965
966 /* destroy a, b */
967 static GEN
Flm_gauss_sp_i(GEN a,GEN b,ulong * detp,ulong p)968 Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
969 {
970 long i, j, k, li, bco, aco = lg(a)-1, s = 1;
971 ulong det = 1;
972 GEN u;
973 ulong pi;
974 if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
975 if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
976 pi = get_Fl_red(p);
977 li = nbrows(a);
978 bco = lg(b)-1;
979 for (i=1; i<=aco; i++)
980 {
981 ulong invpiv;
982 /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
983 for (k = i; k <= li; k++)
984 {
985 ulong piv = ucoeff(a,k,i);
986 if (piv)
987 {
988 ucoeff(a,k,i) = Fl_inv(piv, p);
989 if (detp) det = Fl_mul_pre(det, piv, p, pi);
990 break;
991 }
992 }
993 /* found a pivot on line k */
994 if (k > li) return NULL;
995 if (k != i)
996 { /* swap lines so that k = i */
997 s = -s;
998 for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
999 for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
1000 }
1001 if (i == aco) break;
1002
1003 invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
1004 for (k=i+1; k<=li; k++)
1005 {
1006 ulong m = ucoeff(a,k,i);
1007 if (!m) continue;
1008
1009 m = Fl_mul_pre(m, invpiv, p, pi);
1010 if (m == 1) {
1011 for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
1012 for (j=1; j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
1013 } else {
1014 for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
1015 for (j=1; j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
1016 }
1017 }
1018 }
1019 if (detp)
1020 {
1021 if (s < 0 && det) det = p - det;
1022 *detp = det;
1023 }
1024 u = cgetg(bco+1,t_MAT);
1025 for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
1026 return u;
1027 }
1028
1029 static GEN
Flm_gauss_from_CUP(GEN b,GEN R,GEN C,GEN U,GEN P,ulong p,ulong pi,ulong * detp)1030 Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
1031 {
1032 GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
1033 GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
1034 if (detp)
1035 {
1036 ulong d = perm_sign(P) == 1? 1: p-1;
1037 long i, r = lg(R);
1038 for (i = 1; i < r; i++)
1039 d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
1040 *detp = d;
1041 }
1042 return X;
1043 }
1044
1045 static GEN
Flm_gauss_CUP(GEN a,GEN b,ulong * detp,ulong p)1046 Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
1047 GEN R, C, U, P;
1048 long n = lg(a) - 1, r;
1049 ulong pi = get_Fl_red(p);
1050 if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
1051 return NULL;
1052 return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
1053 }
1054
1055 GEN
Flm_gauss_sp(GEN a,GEN b,ulong * detp,ulong p)1056 Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
1057 pari_sp av = avma;
1058 GEN x;
1059 if (lg(a) - 1 >= Flm_CUP_LIMIT)
1060 x = Flm_gauss_CUP(a, b, detp, p);
1061 else
1062 x = Flm_gauss_sp_i(a, b, detp, p);
1063 if (!x) return gc_NULL(av);
1064 return gerepileupto(av, x);
1065 }
1066
1067 GEN
Flm_gauss(GEN a,GEN b,ulong p)1068 Flm_gauss(GEN a, GEN b, ulong p) {
1069 pari_sp av = avma;
1070 GEN x;
1071 if (lg(a) - 1 >= Flm_CUP_LIMIT)
1072 x = Flm_gauss_CUP(a, b, NULL, p);
1073 else {
1074 a = RgM_shallowcopy(a);
1075 b = RgM_shallowcopy(b);
1076 x = Flm_gauss_sp(a, b, NULL, p);
1077 }
1078 if (!x) return gc_NULL(av);
1079 return gerepileupto(av, x);
1080 }
1081
1082 static GEN
Flm_inv_i(GEN a,ulong * detp,ulong p,long inplace)1083 Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
1084 pari_sp av = avma;
1085 long n = lg(a) - 1;
1086 GEN b, x;
1087 if (n == 0) return cgetg(1, t_MAT);
1088 b = matid_Flm(nbrows(a));
1089 if (n >= Flm_CUP_LIMIT)
1090 x = Flm_gauss_CUP(a, b, detp, p);
1091 else {
1092 if (!inplace)
1093 a = RgM_shallowcopy(a);
1094 x = Flm_gauss_sp(a, b, detp, p);
1095 }
1096 if (!x) return gc_NULL(av);
1097 return gerepileupto(av, x);
1098 }
1099
1100 GEN
Flm_inv_sp(GEN a,ulong * detp,ulong p)1101 Flm_inv_sp(GEN a, ulong *detp, ulong p) {
1102 return Flm_inv_i(a, detp, p, 1);
1103 }
1104
1105 GEN
Flm_inv(GEN a,ulong p)1106 Flm_inv(GEN a, ulong p) {
1107 return Flm_inv_i(a, NULL, p, 0);
1108 }
1109
1110 GEN
Flm_Flc_gauss(GEN a,GEN b,ulong p)1111 Flm_Flc_gauss(GEN a, GEN b, ulong p) {
1112 pari_sp av = avma;
1113 GEN z = Flm_gauss(a, mkmat(b), p);
1114 if (!z) return gc_NULL(av);
1115 if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
1116 return gerepileuptoleaf(av, gel(z,1));
1117 }
1118
1119 GEN
Flm_adjoint(GEN A,ulong p)1120 Flm_adjoint(GEN A, ulong p)
1121 {
1122 pari_sp av = avma;
1123 GEN R, C, U, P, C1, U1, v, c, d;
1124 long r, i, q, n = lg(A)-1, m;
1125 ulong D;
1126 ulong pi = get_Fl_red(p);
1127 if (n == 0) return cgetg(1, t_MAT);
1128 r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
1129 m = nbrows(A);
1130 if (r == n)
1131 {
1132 GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
1133 return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
1134 }
1135 if (r < n-1) return zero_Flm(n, m);
1136 for (q = n, i = 1; i < n; i++)
1137 if (R[i] != i) { q = i; break; }
1138 C1 = matslice(C, 1, q-1, 1, q-1);
1139 c = vecslice(Flm_row(C, q), 1, q-1);
1140 c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
1141 d = cgetg(m+1, t_VECSMALL);
1142 for (i=1; i<q; i++) uel(d,i) = ucoeff(c,1,i);
1143 uel(d,q) = p-1;
1144 for (i=q+1; i<=m; i++) uel(d,i) = 0;
1145 U1 = vecslice(U, 1, n-1);
1146 v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
1147 v = vecsmall_append(v, p-1);
1148 D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
1149 for (i = 1; i <= n-1; i++)
1150 D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
1151 d = Flv_Fl_mul(d, D, p);
1152 return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
1153 }
1154
1155 static GEN
Flm_invimage_CUP(GEN A,GEN B,ulong p)1156 Flm_invimage_CUP(GEN A, GEN B, ulong p) {
1157 pari_sp av = avma;
1158 GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
1159 long r;
1160 ulong pi = get_Fl_red(p);
1161 r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
1162 Rc = indexcompl(R, nbrows(B));
1163 C1 = rowpermute(C, R);
1164 C2 = rowpermute(C, Rc);
1165 B1 = rowpermute(B, R);
1166 B2 = rowpermute(B, Rc);
1167 Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
1168 if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
1169 return NULL;
1170 Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
1171 zero_Flm(lg(A) - 1 - r, lg(B) - 1));
1172 X = rowpermute(Y, perm_inv(P));
1173 return gerepileupto(av, X);
1174 }
1175
1176 GEN
Flm_invimage_i(GEN A,GEN B,ulong p)1177 Flm_invimage_i(GEN A, GEN B, ulong p)
1178 {
1179 GEN d, x, X, Y;
1180 long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
1181
1182 if (!nB) return cgetg(1, t_MAT);
1183 if (maxss(nA, nB) >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
1184 return Flm_invimage_CUP(A, B, p);
1185
1186 x = Flm_ker(shallowconcat(Flm_neg(A,p), B), p);
1187 /* AX = BY, Y in strict upper echelon form with pivots = 1.
1188 * We must find T such that Y T = Id_nB then X T = Z. This exists iff
1189 * Y has at least nB columns and full rank */
1190 nY = lg(x)-1;
1191 if (nY < nB) return NULL;
1192 Y = rowslice(x, nA+1, nA+nB); /* nB rows */
1193 d = cgetg(nB+1, t_VECSMALL);
1194 for (i = nB, j = nY; i >= 1; i--, j--)
1195 {
1196 for (; j>=1; j--)
1197 if (coeff(Y,i,j)) { d[i] = j; break; }
1198 if (!j) return NULL;
1199 }
1200 /* reduce to the case Y square, upper triangular with 1s on diagonal */
1201 Y = vecpermute(Y, d);
1202 x = vecpermute(x, d);
1203 X = rowslice(x, 1, nA);
1204 return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
1205 }
1206 GEN
Flm_invimage(GEN A,GEN B,ulong p)1207 Flm_invimage(GEN A, GEN B, ulong p)
1208 {
1209 pari_sp av = avma;
1210 GEN X = Flm_invimage_i(A,B,p);
1211 if (!X) return gc_NULL(av);
1212 return gerepileupto(av, X);
1213 }
1214
1215 GEN
Flm_Flc_invimage(GEN A,GEN y,ulong p)1216 Flm_Flc_invimage(GEN A, GEN y, ulong p)
1217 {
1218 pari_sp av = avma;
1219 long i, l = lg(A);
1220 GEN M, x;
1221 ulong t;
1222
1223 if (l==1) return NULL;
1224 if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
1225 M = cgetg(l+1,t_MAT);
1226 for (i=1; i<l; i++) gel(M,i) = gel(A,i);
1227 gel(M,l) = y; M = Flm_ker(M,p);
1228 i = lg(M)-1; if (!i) return gc_NULL(av);
1229
1230 x = gel(M,i); t = x[l];
1231 if (!t) return gc_NULL(av);
1232
1233 setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
1234 if (t!=1) x = Flv_Fl_mul(x, t, p);
1235 return gerepileuptoleaf(av, x);
1236 }
1237