Line data Source code
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
19 388512 : Flv_to_ZV(GEN x)
20 4096926 : { pari_APPLY_type(t_VEC, utoi(x[i])) }
21 :
22 : GEN
23 28810047 : Flc_to_ZC(GEN x)
24 240459356 : { pari_APPLY_type(t_COL, utoi(x[i])) }
25 :
26 : GEN
27 11411035 : Flm_to_ZM(GEN x)
28 39272723 : { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) }
29 :
30 : GEN
31 14 : Flc_to_ZC_inplace(GEN z)
32 : {
33 14 : long i, l = lg(z);
34 252 : for (i=1; i<l; i++) gel(z,i) = utoi(z[i]);
35 14 : settyp(z, t_COL);
36 14 : return z;
37 : }
38 :
39 : GEN
40 0 : Flm_to_ZM_inplace(GEN z)
41 : {
42 0 : long i, l = lg(z);
43 0 : for (i=1; i<l; i++) Flc_to_ZC_inplace(gel(z, i));
44 0 : return z;
45 : }
46 :
47 : static GEN
48 1797801 : Flm_solve_upper_1(GEN U, GEN B, ulong p, ulong pi)
49 1797801 : { return Flm_Fl_mul_pre(B, Fl_inv(ucoeff(U, 1, 1), p), p, pi); }
50 :
51 : static GEN
52 2790478 : Flm_rsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
53 : {
54 2790478 : ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
55 2790478 : ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
56 2790590 : ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
57 2790574 : ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
58 2790563 : GEN B1 = rowslice(B, 1, 1);
59 2790438 : GEN B2 = rowslice(B, 2, 2);
60 2790443 : GEN X2 = Flm_Fl_mul_pre(B2, dinv, p, pi);
61 2790508 : GEN X1 = Flm_Fl_mul_pre(Flm_sub(B1, Flm_Fl_mul_pre(X2, b, p, pi), p),
62 : ainv, p, pi);
63 2790519 : return vconcat(X1, X2);
64 : }
65 :
66 : /* solve U*X = B, U upper triangular and invertible */
67 : static GEN
68 6687325 : Flm_rsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
69 : {
70 6687325 : long n = lg(U) - 1, n1;
71 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
72 6687325 : pari_sp av = avma;
73 :
74 6687325 : if (n == 0) return B;
75 6687318 : if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
76 5590900 : if (n == 2) return Flm_rsolve_upper_2(U, B, p, pi);
77 2800426 : n1 = (n + 1)/2;
78 2800426 : U2 = vecslice(U, n1 + 1, n);
79 2801064 : U22 = rowslice(U2, n1 + 1, n);
80 2800985 : B2 = rowslice(B, n1 + 1, n);
81 2800994 : X2 = Flm_rsolve_upper_pre(U22, B2, p, pi);
82 2801162 : U12 = rowslice(U2, 1, n1);
83 2801033 : B1 = rowslice(B, 1, n1);
84 2801021 : B1 = Flm_sub(B1, Flm_mul_pre(U12, X2, p, pi), p);
85 2800962 : if (gc_needed(av, 1)) gerepileall(av, 2, &B1, &X2);
86 2800962 : U11 = matslice(U, 1,n1, 1,n1);
87 2800995 : X1 = Flm_rsolve_upper_pre(U11, B1, p, pi);
88 2801214 : X = vconcat(X1, X2);
89 2801239 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
90 2801239 : return X;
91 : }
92 :
93 : static GEN
94 2537336 : Flm_lsolve_upper_2(GEN U, GEN B, ulong p, ulong pi)
95 : {
96 2537336 : ulong a = ucoeff(U, 1, 1), b = ucoeff(U, 1, 2), d = ucoeff(U, 2, 2);
97 2537336 : ulong D = Fl_mul_pre(a, d, p, pi), Dinv = Fl_inv(D, p);
98 2537364 : ulong ainv = Fl_mul_pre(d, Dinv, p, pi);
99 2537352 : ulong dinv = Fl_mul_pre(a, Dinv, p, pi);
100 2537343 : GEN B1 = vecslice(B, 1, 1);
101 2537328 : GEN B2 = vecslice(B, 2, 2);
102 2537296 : GEN X1 = Flm_Fl_mul_pre(B1, ainv, p, pi);
103 2537335 : GEN X2 = Flm_Fl_mul_pre(Flm_sub(B2, Flm_Fl_mul_pre(X1, b, p, pi), p),
104 : dinv, p, pi);
105 2537330 : return shallowconcat(X1, X2);
106 : }
107 :
108 : /* solve X*U = B, U upper triangular and invertible */
109 : static GEN
110 5466377 : Flm_lsolve_upper_pre(GEN U, GEN B, ulong p, ulong pi)
111 : {
112 5466377 : long n = lg(U) - 1, n1;
113 : GEN U2, U11, U12, U22, B1, B2, X1, X2, X;
114 5466377 : pari_sp av = avma;
115 :
116 5466377 : if (n == 0) return B;
117 5466377 : if (n == 1) return Flm_solve_upper_1(U, B, p, pi);
118 4764948 : if (n == 2) return Flm_lsolve_upper_2(U, B, p, pi);
119 2227612 : n1 = (n + 1)/2;
120 2227612 : U2 = vecslice(U, n1 + 1, n);
121 2227865 : U11 = matslice(U, 1,n1, 1,n1);
122 2227837 : U12 = rowslice(U2, 1, n1);
123 2227839 : U22 = rowslice(U2, n1 + 1, n);
124 2227833 : B1 = vecslice(B, 1, n1);
125 2227809 : B2 = vecslice(B, n1 + 1, n);
126 2227819 : X1 = Flm_lsolve_upper_pre(U11, B1, p, pi);
127 2227842 : B2 = Flm_sub(B2, Flm_mul_pre(X1, U12, p, pi), p);
128 2227821 : if (gc_needed(av, 1)) gerepileall(av, 3, &B2, &U22, &X1);
129 2227821 : X2 = Flm_lsolve_upper_pre(U22, B2, p, pi);
130 2227871 : X = shallowconcat(X1, X2);
131 2227876 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
132 2227876 : return X;
133 : }
134 :
135 : static GEN
136 4795437 : Flm_rsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
137 : {
138 4795437 : GEN X1 = rowslice(A, 1, 1);
139 4795474 : GEN X2 = Flm_sub(rowslice(A, 2, 2),
140 4795394 : Flm_Fl_mul_pre(X1, ucoeff(L, 2, 1), p, pi), p);
141 4795258 : 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
147 10981454 : Flm_rsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
148 : {
149 10981454 : long m = lg(L) - 1, m1, n;
150 : GEN L1, L11, L21, L22, A1, A2, X1, X2, X;
151 10981454 : pari_sp av = avma;
152 :
153 10981454 : if (m == 0) return zero_Flm(0, lg(A) - 1);
154 10981454 : if (m == 1) return rowslice(A, 1, 1);
155 9385907 : if (m == 2) return Flm_rsolve_lower_unit_2(L, A, p, pi);
156 4590474 : m1 = (m + 1)/2;
157 4590474 : n = nbrows(L);
158 4591204 : L1 = vecslice(L, 1, m1);
159 4591223 : L11 = rowslice(L1, 1, m1);
160 4591125 : L21 = rowslice(L1, m1 + 1, n);
161 4591079 : A1 = rowslice(A, 1, m1);
162 4591118 : X1 = Flm_rsolve_lower_unit_pre(L11, A1, p, pi);
163 4591327 : A2 = rowslice(A, m1 + 1, n);
164 4591061 : A2 = Flm_sub(A2, Flm_mul_pre(L21, X1, p, pi), p);
165 4591000 : if (gc_needed(av, 1)) gerepileall(av, 2, &A2, &X1);
166 4591000 : L22 = matslice(L, m1+1,n, m1+1,m);
167 4591094 : X2 = Flm_rsolve_lower_unit_pre(L22, A2, p, pi);
168 4591200 : X = vconcat(X1, X2);
169 4591374 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
170 4591376 : return X;
171 : }
172 :
173 : static GEN
174 1938662 : Flm_lsolve_lower_unit_2(GEN L, GEN A, ulong p, ulong pi)
175 : {
176 1938662 : GEN X2 = vecslice(A, 2, 2);
177 1938665 : GEN X1 = Flm_sub(vecslice(A, 1, 1),
178 1938663 : Flm_Fl_mul_pre(X2, ucoeff(L, 2, 1), p, pi), p);
179 1938664 : return shallowconcat(X1, X2);
180 : }
181 :
182 : /* solve L*X = A, L square lower triangular with ones on the diagonal */
183 : static GEN
184 4503553 : Flm_lsolve_lower_unit_pre(GEN L, GEN A, ulong p, ulong pi)
185 : {
186 4503553 : long m = lg(L) - 1, m1;
187 : GEN L1, L2, L11, L21, L22, A1, A2, X1, X2, X;
188 4503553 : pari_sp av = avma;
189 :
190 4503553 : if (m <= 1) return A;
191 4015028 : if (m == 2) return Flm_lsolve_lower_unit_2(L, A, p, pi);
192 2076369 : m1 = (m + 1)/2;
193 2076369 : L2 = vecslice(L, m1 + 1, m);
194 2076384 : L22 = rowslice(L2, m1 + 1, m);
195 2076378 : A2 = vecslice(A, m1 + 1, m);
196 2076375 : X2 = Flm_lsolve_lower_unit_pre(L22, A2, p, pi);
197 2076383 : if (gc_needed(av, 1)) X2 = gerepilecopy(av, X2);
198 2076383 : L1 = vecslice(L, 1, m1);
199 2076380 : L21 = rowslice(L1, m1 + 1, m);
200 2076382 : A1 = vecslice(A, 1, m1);
201 2076381 : A1 = Flm_sub(A1, Flm_mul_pre(X2, L21, p, pi), p);
202 2076382 : L11 = rowslice(L1, 1, m1);
203 2076382 : if (gc_needed(av, 1)) gerepileall(av, 3, &A1, &L11, &X2);
204 2076382 : X1 = Flm_lsolve_lower_unit_pre(L11, A1, p, pi);
205 2076386 : X = shallowconcat(X1, X2);
206 2076386 : if (gc_needed(av, 1)) X = gerepilecopy(av, X);
207 2076386 : return X;
208 : }
209 :
210 : /* destroy A */
211 : static long
212 3584573 : Flm_CUP_basecase(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
213 : {
214 3584573 : long i, j, k, m = nbrows(A), n = lg(A) - 1, pr, pc;
215 : ulong u, v;
216 :
217 3584574 : if (P) *P = identity_perm(n);
218 3584568 : *R = cgetg(m + 1, t_VECSMALL);
219 17883590 : for (j = 1, pr = 0; j <= n; j++)
220 : {
221 20515614 : for (pr++, pc = 0; pr <= m; pr++)
222 : {
223 138788881 : for (k = j; k <= n; k++) { v = ucoeff(A, pr, k); if (!pc && v) pc = k; }
224 19124558 : if (pc) break;
225 : }
226 15690134 : if (!pc) break;
227 14299030 : (*R)[j] = pr;
228 14299030 : if (pc != j)
229 : {
230 2293193 : swap(gel(A, j), gel(A, pc));
231 2293193 : if (P) lswap((*P)[j], (*P)[pc]);
232 : }
233 14299030 : u = Fl_inv(ucoeff(A, pr, j), p);
234 97134919 : for (i = pr + 1; i <= m; i++)
235 : {
236 82835886 : v = Fl_mul_pre(ucoeff(A, i, j), u, p, pi);
237 82835111 : ucoeff(A, i, j) = v;
238 82835111 : v = Fl_neg(v, p);
239 383308984 : for (k = j + 1; k <= n; k++)
240 300477300 : ucoeff(A, i, k) = Fl_addmul_pre(ucoeff(A, i, k),
241 300473158 : ucoeff(A, pr, k), v, p, pi);
242 : }
243 : }
244 3584560 : setlg(*R, j);
245 3584581 : *C = vecslice(A, 1, j - 1);
246 3584563 : if (U) *U = rowpermute(A, *R);
247 3584562 : return j - 1;
248 : }
249 :
250 : static const long Flm_CUP_LIMIT = 8;
251 :
252 : static long
253 3473322 : Flm_CUP_pre(GEN A, GEN *R, GEN *C, GEN *U, GEN *P, ulong p, ulong pi)
254 : {
255 3473322 : 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 3473319 : pari_sp av = avma;
259 :
260 3473319 : if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
261 : /* destroy A; not called at the outermost recursion level */
262 2435222 : return Flm_CUP_basecase(A, R, C, U, P, p, pi);
263 1038097 : m1 = (minss(m, n) + 1)/2;
264 1038096 : A1 = rowslice(A, 1, m1);
265 1038085 : A2 = rowslice(A, m1 + 1, m);
266 1038084 : r1 = Flm_CUP_pre(A1, &R1, &C1, &U1, &P1, p, pi);
267 1038094 : if (r1 == 0)
268 : {
269 27197 : r2 = Flm_CUP_pre(A2, &R2, &C2, &U2, &P2, p, pi);
270 27197 : *R = cgetg(r2 + 1, t_VECSMALL);
271 36915 : for (i = 1; i <= r2; i++) (*R)[i] = R2[i] + m1;
272 27197 : *C = vconcat(zero_Flm(m1, r2), C2);
273 27197 : *U = U2;
274 27197 : *P = P2;
275 27197 : r = r2;
276 : }
277 : else
278 : {
279 1010897 : U11 = vecslice(U1, 1, r1);
280 1010893 : U12 = vecslice(U1, r1 + 1, n);
281 1010893 : T21 = vecslicepermute(A2, P1, 1, r1);
282 1010894 : T22 = vecslicepermute(A2, P1, r1 + 1, n);
283 1010894 : C21 = Flm_lsolve_upper_pre(U11, T21, p, pi);
284 1010897 : if (gc_needed(av, 1))
285 0 : gerepileall(av, 7, &R1, &C1, &P1, &U11, &U12, &T22, &C21);
286 1010897 : B2 = Flm_sub(T22, Flm_mul_pre(C21, U12, p, pi), p);
287 1010891 : r2 = Flm_CUP_pre(B2, &R2, &C2, &U2, &P2, p, pi);
288 1010887 : r = r1 + r2;
289 1010887 : *R = cgetg(r + 1, t_VECSMALL);
290 6787004 : for (i = 1; i <= r1; i++) (*R)[i] = R1[i];
291 6318358 : for (; i <= r; i++) (*R)[i] = R2[i - r1] + m1;
292 1010885 : *C = shallowconcat(vconcat(C1, C21),
293 : vconcat(zero_Flm(m1, r2), C2));
294 1010899 : *U = shallowconcat(vconcat(U11, zero_Flm(r2, r1)),
295 : vconcat(vecpermute(U12, P2), U2));
296 1010904 : *P = cgetg(n + 1, t_VECSMALL);
297 6787074 : for (i = 1; i <= r1; i++) (*P)[i] = P1[i];
298 12642830 : for (; i <= n; i++) (*P)[i] = P1[P2[i - r1] + r1];
299 : }
300 1038097 : if (gc_needed(av, 1)) gerepileall(av, 4, R, C, U, P);
301 1038098 : return r;
302 : }
303 :
304 : static long
305 1149327 : Flm_echelon_gauss(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
306 1149327 : { 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
310 1065040 : indexcompl(GEN v, long n)
311 : {
312 1065040 : long i, j, k, m = lg(v) - 1;
313 1065040 : GEN w = cgetg(n - m + 1, t_VECSMALL);
314 23071701 : for (i = j = k = 1; i <= n; i++)
315 22006661 : if (j <= m && v[j] == i) j++; else w[k++] = i;
316 1065040 : return w;
317 : }
318 :
319 : /* column echelon form */
320 : static long
321 1943290 : Flm_echelon_pre(GEN A, GEN *R, GEN *C, ulong p, ulong pi)
322 : {
323 1943290 : 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 1943290 : pari_sp av = avma;
327 :
328 1943290 : if (m < Flm_CUP_LIMIT || n < Flm_CUP_LIMIT)
329 1149326 : return Flm_echelon_gauss(Flm_copy(A), R, C, p, pi);
330 :
331 793964 : n1 = (n + 1)/2;
332 793964 : A1 = vecslice(A, 1, n1);
333 793964 : A2 = vecslice(A, n1 + 1, n);
334 793964 : r1 = Flm_echelon_pre(A1, &R1, &C1, p, pi);
335 793964 : if (!r1) return Flm_echelon_pre(A2, R, C, p, pi);
336 719085 : if (r1 == m) { *R = R1; *C = C1; return r1; }
337 :
338 713960 : R1c = indexcompl(R1, m);
339 713960 : C11 = rowpermute(C1, R1);
340 713956 : C21 = rowpermute(C1, R1c);
341 713959 : A12 = rowpermute(A2, R1);
342 713959 : A22 = rowpermute(A2, R1c);
343 713959 : M12 = Flm_rsolve_lower_unit_pre(C11, A12, p, pi);
344 713961 : B2 = Flm_sub(A22, Flm_mul_pre(C21, M12, p, pi), p);
345 713960 : r2 = Flm_echelon_pre(B2, &R2, &C2, p, pi);
346 713960 : if (!r2) { *R = R1; *C = C1; r = r1; }
347 : else
348 : {
349 668506 : R2 = perm_mul(R1c, R2);
350 668506 : C2 = rowpermute(vconcat(zero_Flm(r1, r2), C2),
351 : perm_inv(vecsmall_concat(R1, R1c)));
352 668507 : r = r1 + r2;
353 668507 : *R = cgetg(r + 1, t_VECSMALL);
354 668507 : *C = cgetg(r + 1, t_MAT);
355 8380000 : for (j = j1 = j2 = 1; j <= r; j++)
356 7711494 : if (j2 > r2 || (j1 <= r1 && R1[j1] < R2[j2]))
357 : {
358 4083177 : gel(*C, j) = gel(C1, j1);
359 4083177 : (*R)[j] = R1[j1++];
360 : }
361 : else
362 : {
363 3628317 : gel(*C, j) = gel(C2, j2);
364 3628317 : (*R)[j] = R2[j2++];
365 : }
366 : }
367 713960 : if (gc_needed(av, 1)) gerepileall(av, 2, R, C);
368 713960 : return r;
369 : }
370 :
371 : static void /* assume m < p */
372 11845659 : _Fl_addmul(GEN b, long k, long i, ulong m, ulong p, ulong pi)
373 : {
374 11845659 : uel(b,k) = Fl_addmul_pre(uel(b, k), m, uel(b, i), p, pi);
375 11845690 : }
376 : static void /* same m = 1 */
377 667137 : _Fl_add(GEN b, long k, long i, ulong p)
378 : {
379 667137 : uel(b,k) = Fl_add(uel(b,k), uel(b,i), p);
380 667137 : }
381 : static void /* assume m < p && SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
382 4427083 : _Fl_addmul_OK(GEN b, long k, long i, ulong m, ulong p)
383 : {
384 4427083 : uel(b,k) += m * uel(b,i);
385 4427083 : if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
386 4427083 : }
387 : static void /* assume SMALL_ULONG(p) && (! (b[i] & b[k] & HIGHMASK)) */
388 2196555 : _Fl_add_OK(GEN b, long k, long i, ulong p)
389 : {
390 2196555 : uel(b,k) += uel(b,i);
391 2196555 : if (uel(b,k) & HIGHMASK) uel(b,k) %= p;
392 2196555 : }
393 :
394 : /* assume 0 <= a[i,j] < p */
395 : static GEN
396 440343 : Fl_get_col_OK(GEN a, GEN b, long li, ulong p)
397 : {
398 440343 : GEN u = cgetg(li+1,t_VECSMALL);
399 440342 : ulong m = uel(b,li) % p;
400 : long i,j;
401 :
402 440342 : uel(u,li) = (m * ucoeff(a,li,li)) % p;
403 1700402 : for (i = li-1; i > 0; i--)
404 : {
405 1260060 : m = p - uel(b,i)%p;
406 4040111 : for (j = i+1; j <= li; j++) {
407 2780051 : if (m & HIGHBIT) m %= p;
408 2780051 : m += ucoeff(a,i,j) * uel(u,j); /* 0 <= u[j] < p */
409 : }
410 1260060 : m %= p;
411 1260060 : if (m) m = ((p-m) * ucoeff(a,i,i)) % p;
412 1260060 : uel(u,i) = m;
413 : }
414 440342 : return u;
415 : }
416 : static GEN
417 2163983 : Fl_get_col(GEN a, GEN b, long li, ulong p)
418 : {
419 2163983 : GEN u = cgetg(li+1,t_VECSMALL);
420 2163981 : ulong m = uel(b,li) % p;
421 : long i,j;
422 :
423 2163981 : uel(u,li) = Fl_mul(m, ucoeff(a,li,li), p);
424 5726424 : for (i=li-1; i>0; i--)
425 : {
426 3562441 : m = b[i]%p;
427 9118016 : for (j = i+1; j <= li; j++)
428 5555575 : m = Fl_sub(m, Fl_mul(ucoeff(a,i,j), uel(u,j), p), p);
429 3562441 : if (m) m = Fl_mul(m, ucoeff(a,i,i), p);
430 3562442 : uel(u,i) = m;
431 : }
432 2163983 : return u;
433 : }
434 :
435 : static GEN
436 667223 : 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 667223 : n = lg(x)-1;
443 667223 : m=nbrows(x); r=0;
444 :
445 667223 : c = zero_zv(m);
446 667222 : d = cgetg(n+1, t_VECSMALL);
447 667220 : a = 0; /* for gcc -Wall */
448 3277227 : for (k=1; k<=n; k++)
449 : {
450 9459651 : for (j=1; j<=m; j++)
451 8271785 : if (!c[j])
452 : {
453 4441806 : a = ucoeff(x,j,k) % p;
454 4441806 : if (a) break;
455 : }
456 2663899 : if (j > m)
457 : {
458 1187873 : if (deplin==1) {
459 53896 : c = cgetg(n+1, t_VECSMALL);
460 176185 : for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k) % p;
461 105221 : c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
462 53896 : return c;
463 : }
464 1133977 : r++; d[k] = 0;
465 : }
466 : else
467 : {
468 1476026 : ulong piv = p - Fl_inv(a, p); /* -1/a */
469 1476030 : c[j] = k; d[k] = j;
470 1476030 : ucoeff(x,j,k) = p-1;
471 1476030 : if (piv != 1)
472 3181515 : for (i=k+1; i<=n; i++) ucoeff(x,j,i) = (piv * ucoeff(x,j,i)) % p;
473 7591092 : for (t=1; t<=m; t++)
474 : {
475 6115062 : if (t == j) continue;
476 :
477 4639035 : piv = ( ucoeff(x,t,k) %= p );
478 4639035 : if (!piv) continue;
479 2169412 : if (piv == 1)
480 2548818 : for (i=k+1; i<=n; i++) _Fl_add_OK(gel(x,i),t,j, p);
481 : else
482 4888399 : for (i=k+1; i<=n; i++) _Fl_addmul_OK(gel(x,i),t,j,piv, p);
483 : }
484 : }
485 : }
486 613328 : if (deplin==1) return NULL;
487 :
488 613321 : y = cgetg(r+1, t_MAT);
489 1747296 : for (j=k=1; j<=r; j++,k++)
490 : {
491 1133977 : GEN C = cgetg(n+1, t_VECSMALL);
492 :
493 2202418 : gel(y,j) = C; while (d[k]) k++;
494 9357963 : for (i=1; i<k; i++)
495 8223988 : if (d[i])
496 2337276 : uel(C,i) = ucoeff(x,d[i],k) % p;
497 : else
498 5886712 : uel(C,i) = 0UL;
499 7745623 : uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
500 : }
501 613319 : if (deplin == 2) {
502 0 : GEN pc = cgetg(n - r + 1, t_VECSMALL); /* indices of pivot columns */
503 0 : for (i = j = 1; j <= n; j++)
504 0 : if (d[j]) pc[i++] = j;
505 0 : return mkvec2(y, pc);
506 : }
507 613319 : return y;
508 : }
509 :
510 : /* in place, destroy x */
511 : static GEN
512 814542 : 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 814542 : n = lg(x)-1;
518 814542 : if (!n) return cgetg(1,t_MAT);
519 814535 : if (SMALL_ULONG(p)) return Flm_ker_gauss_OK(x, p, deplin);
520 147312 : pi = get_Fl_red(p);
521 :
522 147312 : m=nbrows(x); r=0;
523 :
524 147312 : c = zero_zv(m);
525 147312 : d = cgetg(n+1, t_VECSMALL);
526 147312 : a = 0; /* for gcc -Wall */
527 474560 : for (k=1; k<=n; k++)
528 : {
529 854525 : for (j=1; j<=m; j++)
530 751176 : if (!c[j])
531 : {
532 567251 : a = ucoeff(x,j,k);
533 567251 : if (a) break;
534 : }
535 327255 : if (j > m)
536 : {
537 103349 : if (deplin==1) {
538 7 : c = cgetg(n+1, t_VECSMALL);
539 21 : for (i=1; i<k; i++) c[i] = ucoeff(x,d[i],k);
540 7 : c[k] = 1; for (i=k+1; i<=n; i++) c[i] = 0;
541 7 : return c;
542 : }
543 103342 : r++; d[k] = 0;
544 : }
545 : else
546 : {
547 223906 : ulong piv = p - Fl_inv(a, p); /* -1/a */
548 223906 : c[j] = k; d[k] = j;
549 223906 : ucoeff(x,j,k) = p-1;
550 223906 : if (piv != 1)
551 405220 : for (i=k+1; i<=n; i++)
552 187237 : ucoeff(x,j,i) = Fl_mul_pre(piv, ucoeff(x,j,i), p, pi);
553 971929 : for (t=1; t<=m; t++)
554 : {
555 748023 : if (t == j) continue;
556 :
557 524117 : piv = ucoeff(x,t,k);
558 524117 : if (!piv) continue;
559 308616 : if (piv == 1)
560 46124 : for (i=k+1; i<=n; i++) _Fl_add(gel(x,i),t,j,p);
561 : else
562 710347 : for (i=k+1; i<=n; i++) _Fl_addmul(gel(x,i),t,j,piv,p, pi);
563 : }
564 : }
565 : }
566 147305 : if (deplin==1) return NULL;
567 :
568 147298 : y = cgetg(r+1, t_MAT);
569 250640 : for (j=k=1; j<=r; j++,k++)
570 : {
571 103342 : GEN C = cgetg(n+1, t_VECSMALL);
572 :
573 166193 : gel(y,j) = C; while (d[k]) k++;
574 218087 : for (i=1; i<k; i++)
575 114745 : if (d[i])
576 79019 : uel(C,i) = ucoeff(x,d[i],k);
577 : else
578 35726 : uel(C,i) = 0UL;
579 166977 : uel(C,k) = 1UL; for (i=k+1; i<=n; i++) uel(C,i) = 0UL;
580 : }
581 147298 : if (deplin == 2) {
582 142444 : GEN pc = cgetg(n - r + 1, t_VECSMALL); /* indices of pivot columns */
583 442500 : for (i = j = 1; j <= n; j++)
584 300056 : if (d[j]) pc[i++] = j;
585 142444 : return mkvec2(y, pc);
586 : }
587 4854 : return y;
588 : }
589 :
590 : GEN
591 271610 : Flm_intersect_i(GEN x, GEN y, ulong p)
592 : {
593 271610 : long j, lx = lg(x);
594 : GEN z;
595 :
596 271610 : if (lx==1 || lg(y)==1) return cgetg(1,t_MAT);
597 271596 : z = Flm_ker(shallowconcat(x,y), p);
598 1470154 : for (j=lg(z)-1; j; j--) setlg(gel(z,j),lx);
599 271596 : return Flm_mul(x,z,p);
600 : }
601 : GEN
602 0 : Flm_intersect(GEN x, GEN y, ulong p)
603 : {
604 0 : pari_sp av = avma;
605 0 : return gerepileupto(av, Flm_image(Flm_intersect_i(x, y, p), p));
606 : }
607 :
608 : static GEN
609 319997 : Flm_ker_echelon(GEN x, ulong p, long pivots) {
610 319997 : pari_sp av = avma;
611 : GEN R, Rc, C, C1, C2, S, K;
612 319997 : long n = lg(x) - 1, r;
613 319997 : ulong pi = get_Fl_red(p);
614 319997 : r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
615 319997 : Rc = indexcompl(R, n);
616 319997 : C1 = rowpermute(C, R);
617 319997 : C2 = rowpermute(C, Rc);
618 319997 : S = Flm_lsolve_lower_unit_pre(C1, C2, p, pi);
619 319997 : K = vecpermute(shallowconcat(Flm_neg(S, p), matid_Flm(n - r)),
620 : perm_inv(vecsmall_concat(R, Rc)));
621 319997 : K = Flm_transpose(K);
622 319996 : if (pivots)
623 32218 : return gerepilecopy(av, mkvec2(K, R));
624 287778 : return gerepileupto(av, K);
625 : }
626 :
627 : static GEN
628 30803 : Flm_deplin_echelon(GEN x, ulong p) {
629 30803 : pari_sp av = avma;
630 : GEN R, Rc, C, C1, C2, s, v;
631 30803 : long i, n = lg(x) - 1, r;
632 30803 : ulong pi = get_Fl_red(p);
633 30803 : r = Flm_echelon_pre(Flm_transpose(x), &R, &C, p, pi);
634 30803 : if (r == n) return gc_NULL(av);
635 30796 : Rc = indexcompl(R, n);
636 30796 : i = Rc[1];
637 30796 : C1 = rowpermute(C, R);
638 30796 : C2 = rowslice(C, i, i);
639 30796 : s = Flm_row(Flm_lsolve_lower_unit_pre(C1, C2, p, pi), 1);
640 30796 : v = vecsmallpermute(vecsmall_concat(Flv_neg(s, p), vecsmall_ei(n - r, 1)),
641 : perm_inv(vecsmall_concat(R, Rc)));
642 30796 : return gerepileuptoleaf(av, v);
643 : }
644 :
645 : static GEN
646 1165345 : Flm_ker_i(GEN x, ulong p, long deplin, long inplace) {
647 1165345 : if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
648 350800 : switch(deplin) {
649 287778 : case 0: return Flm_ker_echelon(x, p, 0);
650 30803 : case 1: return Flm_deplin_echelon(x, p);
651 32219 : case 2: return Flm_ker_echelon(x, p, 1);
652 : }
653 814545 : return Flm_ker_gauss(inplace? x: Flm_copy(x), p, deplin);
654 : }
655 :
656 : GEN
657 615289 : Flm_ker_sp(GEN x, ulong p, long deplin) { return Flm_ker_i(x, p, deplin, 1); }
658 : GEN
659 550056 : Flm_ker(GEN x, ulong p) { return Flm_ker_i(x, p, 0, 0); }
660 : GEN
661 0 : Flm_deplin(GEN x, ulong p) { return Flm_ker_i(x, p, 1, 0); }
662 :
663 : /* in place, destroy a, SMALL_ULONG(p) is TRUE */
664 : static ulong
665 2310 : Flm_det_gauss_OK(GEN a, long nbco, ulong p)
666 : {
667 2310 : long i,j,k, s = 1;
668 2310 : ulong q, x = 1;
669 :
670 9534 : for (i=1; i<nbco; i++)
671 : {
672 8841 : for(k=i; k<=nbco; k++)
673 : {
674 8659 : ulong c = ucoeff(a,k,i) % p;
675 8659 : ucoeff(a,k,i) = c;
676 8659 : if (c) break;
677 : }
678 23191 : for(j=k+1; j<=nbco; j++) ucoeff(a,j,i) %= p;
679 7406 : if (k > nbco) return ucoeff(a,i,i);
680 7224 : if (k != i)
681 : { /* exchange the lines s.t. k = i */
682 3122 : for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
683 784 : s = -s;
684 : }
685 7224 : q = ucoeff(a,i,i);
686 :
687 7224 : if (x & HIGHMASK) x %= p;
688 7224 : x *= q;
689 7224 : q = Fl_inv(q,p);
690 24080 : for (k=i+1; k<=nbco; k++)
691 : {
692 16856 : ulong m = ucoeff(a,i,k) % p;
693 16856 : if (!m) continue;
694 :
695 9618 : m = p - ((m*q)%p);
696 37576 : for (j=i+1; j<=nbco; j++)
697 : {
698 27958 : ulong c = ucoeff(a,j,k);
699 27958 : if (c & HIGHMASK) c %= p;
700 27958 : ucoeff(a,j,k) = c + m*ucoeff(a,j,i);
701 : }
702 : }
703 : }
704 2128 : if (x & HIGHMASK) x %= p;
705 2128 : q = ucoeff(a,nbco,nbco);
706 2128 : if (q & HIGHMASK) q %= p;
707 2128 : x = (x*q) % p;
708 2128 : if (s < 0 && x) x = p - x;
709 2128 : return x;
710 : }
711 :
712 : /* in place, destroy a */
713 : static ulong
714 53345 : Flm_det_gauss(GEN a, ulong p)
715 : {
716 53345 : long i,j,k, s = 1, nbco = lg(a)-1;
717 53345 : ulong pi, q, x = 1;
718 :
719 53345 : if (SMALL_ULONG(p)) return Flm_det_gauss_OK(a, nbco, p);
720 51035 : pi = get_Fl_red(p);
721 325543 : for (i=1; i<nbco; i++)
722 : {
723 274910 : for(k=i; k<=nbco; k++)
724 274910 : if (ucoeff(a,k,i)) break;
725 274507 : if (k > nbco) return ucoeff(a,i,i);
726 274507 : if (k != i)
727 : { /* exchange the lines s.t. k = i */
728 1090 : for (j=i; j<=nbco; j++) lswap(ucoeff(a,i,j), ucoeff(a,k,j));
729 212 : s = -s;
730 : }
731 274507 : q = ucoeff(a,i,i);
732 :
733 274507 : x = Fl_mul_pre(x, q, p, pi);
734 274507 : q = Fl_inv(q,p);
735 1176234 : for (k=i+1; k<=nbco; k++)
736 : {
737 901726 : ulong m = ucoeff(a,i,k);
738 901726 : if (!m) continue;
739 :
740 842424 : m = Fl_mul_pre(m, q, p, pi);
741 4329374 : for (j=i+1; j<=nbco; j++)
742 3486950 : ucoeff(a,j,k) = Fl_sub(ucoeff(a,j,k), Fl_mul_pre(m,ucoeff(a,j,i), p, pi), p);
743 : }
744 : }
745 51036 : if (s < 0) x = Fl_neg(x, p);
746 51036 : return Fl_mul(x, ucoeff(a,nbco,nbco), p);
747 : }
748 :
749 : static ulong
750 37377 : Flm_det_CUP(GEN a, ulong p) {
751 : GEN R, C, U, P;
752 37377 : long i, n = lg(a) - 1, r;
753 : ulong d;
754 37377 : ulong pi = get_Fl_red(p);
755 37377 : r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi);
756 37378 : if (r < n)
757 42 : d = 0;
758 : else {
759 37336 : d = perm_sign(P) == 1? 1: p-1;
760 429901 : for (i = 1; i <= n; i++)
761 392565 : d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
762 : }
763 37378 : return d;
764 : }
765 :
766 : static ulong
767 90722 : Flm_det_i(GEN x, ulong p, long inplace) {
768 90722 : pari_sp av = avma;
769 : ulong d;
770 90722 : if (lg(x) - 1 >= Flm_CUP_LIMIT)
771 37377 : d = Flm_det_CUP(x, p);
772 : else
773 53345 : d = Flm_det_gauss(inplace? x: Flm_copy(x), p);
774 90722 : return gc_ulong(av, d);
775 : }
776 :
777 : ulong
778 90722 : Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); }
779 : ulong
780 0 : Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); }
781 :
782 : /* Destroy x */
783 : static GEN
784 4397343 : Flm_gauss_pivot(GEN x, ulong p, long *rr)
785 : {
786 : GEN c,d;
787 : long i,j,k,r,t,n,m;
788 :
789 4397343 : n=lg(x)-1; if (!n) { *rr=0; return NULL; }
790 :
791 4396034 : m=nbrows(x); r=0;
792 4396029 : d=cgetg(n+1,t_VECSMALL);
793 4396018 : c = zero_zv(m);
794 18907731 : for (k=1; k<=n; k++)
795 : {
796 36109802 : for (j=1; j<=m; j++)
797 34309058 : if (!c[j])
798 : {
799 18418994 : ucoeff(x,j,k) %= p;
800 18418994 : if (ucoeff(x,j,k)) break;
801 : }
802 14511674 : if (j>m) { r++; d[k]=0; }
803 : else
804 : {
805 12710885 : ulong piv = p - Fl_inv(ucoeff(x,j,k), p);
806 12710968 : c[j]=k; d[k]=j;
807 29214585 : for (i=k+1; i<=n; i++)
808 16503620 : ucoeff(x,j,i) = Fl_mul(piv, ucoeff(x,j,i), p);
809 60432003 : for (t=1; t<=m; t++)
810 47721201 : if (!c[t]) /* no pivot on that line yet */
811 : {
812 20430924 : piv = ucoeff(x,t,k);
813 20430924 : if (piv)
814 : {
815 12029285 : ucoeff(x,t,k) = 0;
816 37274887 : for (i=k+1; i<=n; i++)
817 25245741 : ucoeff(x,t,i) = Fl_add(ucoeff(x,t,i),
818 25245765 : Fl_mul(piv,ucoeff(x,j,i),p),p);
819 : }
820 : }
821 41925311 : for (i=k; i<=n; i++) ucoeff(x,j,i) = 0; /* dummy */
822 : }
823 : }
824 4396057 : *rr = r; return gc_const((pari_sp)d, d);
825 : }
826 :
827 : static GEN
828 273930 : Flm_pivots_CUP(GEN x, ulong p, long *rr)
829 : {
830 273930 : long i, n = lg(x) - 1, r;
831 273930 : GEN R, C, U, P, d = zero_zv(n);
832 273930 : ulong pi = get_Fl_red(p);
833 273930 : r = Flm_CUP_pre(x, &R, &C, &U, &P, p, pi);
834 3017425 : for(i = 1; i <= r; i++)
835 2743495 : d[P[i]] = R[i];
836 273930 : *rr = n - r; return gc_const((pari_sp)d, d);
837 : }
838 :
839 : GEN
840 4671272 : Flm_pivots(GEN x, ulong p, long *rr, long inplace)
841 : {
842 4671272 : if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT)
843 273930 : return Flm_pivots_CUP(x, p, rr);
844 4397342 : return Flm_gauss_pivot(inplace? x: Flm_copy(x), p, rr);
845 : }
846 :
847 : long
848 39844 : Flm_rank(GEN x, ulong p)
849 : {
850 39844 : pari_sp av = avma;
851 : long r;
852 39844 : if (lg(x) - 1 >= Flm_CUP_LIMIT && nbrows(x) >= Flm_CUP_LIMIT) {
853 : GEN R, C;
854 9688 : ulong pi = get_Fl_red(p);
855 9688 : return gc_long(av, Flm_echelon_pre(x, &R, &C, p, pi));
856 : }
857 30156 : (void) Flm_pivots(x, p, &r, 0);
858 30156 : return gc_long(av, lg(x)-1 - r);
859 : }
860 :
861 : /* assume dim A >= 1, A invertible + upper triangular, 1s on diagonal,
862 : * reduced mod p */
863 : static GEN
864 1274 : Flm_inv_upper_1_ind(GEN A, long index, ulong p)
865 : {
866 1274 : long n = lg(A)-1, i = index, j;
867 1274 : GEN u = const_vecsmall(n, 0);
868 1274 : u[i] = 1;
869 1274 : if (SMALL_ULONG(p))
870 3269 : for (i--; i>0; i--)
871 : {
872 2009 : ulong m = ucoeff(A,i,i+1) * uel(u,i+1); /* j = i+1 */
873 7035 : for (j=i+2; j<=n; j++)
874 : {
875 5026 : if (m & HIGHMASK) m %= p;
876 5026 : m += ucoeff(A,i,j) * uel(u,j);
877 : }
878 2009 : u[i] = Fl_neg(m % p, p);
879 : }
880 : else
881 21 : for (i--; i>0; i--)
882 : {
883 7 : ulong m = Fl_mul(ucoeff(A,i,i+1),uel(u,i+1), p); /* j = i+1 */
884 7 : for (j=i+2; j<=n; j++) m = Fl_add(m, Fl_mul(ucoeff(A,i,j),uel(u,j),p), p);
885 7 : u[i] = Fl_neg(m, p);
886 : }
887 1274 : return u;
888 : }
889 : static GEN
890 490 : Flm_inv_upper_1(GEN A, ulong p)
891 : {
892 : long i, l;
893 490 : GEN B = cgetg_copy(A, &l);
894 1764 : for (i = 1; i < l; i++) gel(B,i) = Flm_inv_upper_1_ind(A, i, p);
895 490 : return B;
896 : }
897 : /* destroy a, b */
898 : static GEN
899 118584 : Flm_gauss_sp_OK(GEN a, GEN b, ulong *detp, ulong p)
900 : {
901 118584 : long i, j, k, li, bco, aco = lg(a)-1, s = 1;
902 118584 : ulong det = 1;
903 : GEN u;
904 :
905 118584 : li = nbrows(a);
906 118584 : bco = lg(b)-1;
907 416987 : for (i=1; i<=aco; i++)
908 : {
909 : ulong invpiv;
910 : /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
911 1028328 : for (k = 1; k < i; k++) ucoeff(a,k,i) %= p;
912 473688 : for (k = i; k <= li; k++)
913 : {
914 473675 : ulong piv = ( ucoeff(a,k,i) %= p );
915 473675 : if (piv)
916 : {
917 416969 : ucoeff(a,k,i) = Fl_inv(piv, p);
918 416975 : if (detp)
919 : {
920 0 : if (det & HIGHMASK) det %= p;
921 0 : det *= piv;
922 : }
923 416975 : break;
924 : }
925 : }
926 : /* found a pivot on line k */
927 416988 : if (k > li) return NULL;
928 416981 : if (k != i)
929 : { /* swap lines so that k = i */
930 38045 : s = -s;
931 157275 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
932 205989 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
933 : }
934 416981 : if (i == aco) break;
935 :
936 298404 : invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
937 947176 : for (k=i+1; k<=li; k++)
938 : {
939 648773 : ulong m = ( ucoeff(a,k,i) %= p );
940 648773 : if (!m) continue;
941 :
942 183362 : m = Fl_mul(m, invpiv, p);
943 183361 : if (m == 1) {
944 142454 : for (j=i+1; j<=aco; j++) _Fl_add_OK(gel(a,j),k,i, p);
945 225301 : for (j=1; j<=bco; j++) _Fl_add_OK(gel(b,j),k,i, p);
946 : } else {
947 509167 : for (j=i+1; j<=aco; j++) _Fl_addmul_OK(gel(a,j),k,i,m, p);
948 845634 : for (j=1; j<=bco; j++) _Fl_addmul_OK(gel(b,j),k,i,m, p);
949 : }
950 : }
951 : }
952 118582 : if (detp)
953 : {
954 0 : det %= p;
955 0 : if (s < 0 && det) det = p - det;
956 0 : *detp = det;
957 : }
958 118582 : u = cgetg(bco+1,t_MAT);
959 558919 : for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col_OK(a,gel(b,j), aco,p);
960 118576 : return u;
961 : }
962 :
963 : /* destroy a, b */
964 : static GEN
965 2157962 : Flm_gauss_sp_i(GEN a, GEN b, ulong *detp, ulong p)
966 : {
967 2157962 : long i, j, k, li, bco, aco = lg(a)-1, s = 1;
968 2157962 : ulong det = 1;
969 : GEN u;
970 : ulong pi;
971 2157962 : if (!aco) { if (detp) *detp = 1; return cgetg(1,t_MAT); }
972 2157962 : if (SMALL_ULONG(p)) return Flm_gauss_sp_OK(a, b, detp, p);
973 2039379 : pi = get_Fl_red(p);
974 2039381 : li = nbrows(a);
975 2039380 : bco = lg(b)-1;
976 5308001 : for (i=1; i<=aco; i++)
977 : {
978 : ulong invpiv;
979 : /* Fl_get_col wants 0 <= a[i,j] < p for all i,j */
980 5547168 : for (k = i; k <= li; k++)
981 : {
982 5547168 : ulong piv = ucoeff(a,k,i);
983 5547168 : if (piv)
984 : {
985 5308005 : ucoeff(a,k,i) = Fl_inv(piv, p);
986 5307999 : if (detp) det = Fl_mul_pre(det, piv, p, pi);
987 5307999 : break;
988 : }
989 : }
990 : /* found a pivot on line k */
991 5307999 : if (k > li) return NULL;
992 5307999 : if (k != i)
993 : { /* swap lines so that k = i */
994 182103 : s = -s;
995 677958 : for (j=i; j<=aco; j++) swap(gcoeff(a,i,j), gcoeff(a,k,j));
996 369941 : for (j=1; j<=bco; j++) swap(gcoeff(b,i,j), gcoeff(b,k,j));
997 : }
998 5307999 : if (i == aco) break;
999 :
1000 3268615 : invpiv = p - ucoeff(a,i,i); /* -1/piv mod p */
1001 8247570 : for (k=i+1; k<=li; k++)
1002 : {
1003 4978949 : ulong m = ucoeff(a,k,i);
1004 4978949 : if (!m) continue;
1005 :
1006 4141225 : m = Fl_mul_pre(m, invpiv, p, pi);
1007 4141228 : if (m == 1) {
1008 637886 : for (j=i+1; j<=aco; j++) _Fl_add(gel(a,j),k,i, p);
1009 443501 : for (j=1; j<=bco; j++) _Fl_add(gel(b,j),k,i, p);
1010 : } else {
1011 11370117 : for (j=i+1; j<=aco; j++) _Fl_addmul(gel(a,j),k,i,m, p, pi);
1012 7895936 : for (j=1; j<=bco; j++) _Fl_addmul(gel(b,j),k,i,m, p, pi);
1013 : }
1014 : }
1015 : }
1016 2039380 : if (detp)
1017 : {
1018 0 : if (s < 0 && det) det = p - det;
1019 0 : *detp = det;
1020 : }
1021 2039380 : u = cgetg(bco+1,t_MAT);
1022 4203363 : for (j=1; j<=bco; j++) gel(u,j) = Fl_get_col(a,gel(b,j), aco,p);
1023 2039382 : return u;
1024 : }
1025 :
1026 : static GEN
1027 1085542 : Flm_gauss_from_CUP(GEN b, GEN R, GEN C, GEN U, GEN P, ulong p, ulong pi, ulong *detp)
1028 : {
1029 1085542 : GEN Y = Flm_rsolve_lower_unit_pre(rowpermute(C, R), rowpermute(b, R), p, pi);
1030 1085550 : GEN X = rowpermute(Flm_rsolve_upper_pre(U, Y, p, pi), perm_inv(P));
1031 1085550 : if (detp)
1032 : {
1033 918316 : ulong d = perm_sign(P) == 1? 1: p-1;
1034 918316 : long i, r = lg(R);
1035 5670982 : for (i = 1; i < r; i++)
1036 4752665 : d = Fl_mul_pre(d, ucoeff(U, i, i), p, pi);
1037 918317 : *detp = d;
1038 : }
1039 1085551 : return X;
1040 : }
1041 :
1042 : static GEN
1043 167245 : Flm_gauss_CUP(GEN a, GEN b, ulong *detp, ulong p) {
1044 : GEN R, C, U, P;
1045 167245 : long n = lg(a) - 1, r;
1046 167245 : ulong pi = get_Fl_red(p);
1047 167245 : if (nbrows(a) < n || (r = Flm_CUP_pre(a, &R, &C, &U, &P, p, pi)) < n)
1048 14 : return NULL;
1049 167231 : return Flm_gauss_from_CUP(b, R, C, U, P, p, pi, detp);
1050 : }
1051 :
1052 : GEN
1053 2205148 : Flm_gauss_sp(GEN a, GEN b, ulong *detp, ulong p) {
1054 2205148 : pari_sp av = avma;
1055 : GEN x;
1056 2205148 : if (lg(a) - 1 >= Flm_CUP_LIMIT)
1057 47186 : x = Flm_gauss_CUP(a, b, detp, p);
1058 : else
1059 2157962 : x = Flm_gauss_sp_i(a, b, detp, p);
1060 2205153 : if (!x) return gc_NULL(av);
1061 2205146 : return gerepileupto(av, x);
1062 : }
1063 :
1064 : GEN
1065 2143698 : Flm_gauss(GEN a, GEN b, ulong p) {
1066 2143698 : pari_sp av = avma;
1067 : GEN x;
1068 2143698 : if (lg(a) - 1 >= Flm_CUP_LIMIT)
1069 104327 : x = Flm_gauss_CUP(a, b, NULL, p);
1070 : else {
1071 2039371 : a = RgM_shallowcopy(a);
1072 2039368 : b = RgM_shallowcopy(b);
1073 2039369 : x = Flm_gauss_sp(a, b, NULL, p);
1074 : }
1075 2143700 : if (!x) return gc_NULL(av);
1076 2143686 : return gerepileupto(av, x);
1077 : }
1078 :
1079 : static GEN
1080 63628 : Flm_inv_i(GEN a, ulong *detp, ulong p, long inplace) {
1081 63628 : pari_sp av = avma;
1082 63628 : long n = lg(a) - 1;
1083 : GEN b, x;
1084 63628 : if (n == 0) return cgetg(1, t_MAT);
1085 63628 : b = matid_Flm(nbrows(a));
1086 63628 : if (n >= Flm_CUP_LIMIT)
1087 15732 : x = Flm_gauss_CUP(a, b, detp, p);
1088 : else {
1089 47896 : if (!inplace)
1090 46300 : a = RgM_shallowcopy(a);
1091 47896 : x = Flm_gauss_sp(a, b, detp, p);
1092 : }
1093 63628 : if (!x) return gc_NULL(av);
1094 63621 : return gerepileupto(av, x);
1095 : }
1096 :
1097 : GEN
1098 1827 : Flm_inv_sp(GEN a, ulong *detp, ulong p) {
1099 1827 : return Flm_inv_i(a, detp, p, 1);
1100 : }
1101 :
1102 : GEN
1103 61801 : Flm_inv(GEN a, ulong p) {
1104 61801 : return Flm_inv_i(a, NULL, p, 0);
1105 : }
1106 :
1107 : GEN
1108 21 : Flm_Flc_gauss(GEN a, GEN b, ulong p) {
1109 21 : pari_sp av = avma;
1110 21 : GEN z = Flm_gauss(a, mkmat(b), p);
1111 21 : if (!z) return gc_NULL(av);
1112 14 : if (lg(z) == 1) { set_avma(av); return cgetg(1,t_VECSMALL); }
1113 14 : return gerepileuptoleaf(av, gel(z,1));
1114 : }
1115 :
1116 : GEN
1117 918328 : Flm_adjoint(GEN A, ulong p)
1118 : {
1119 918328 : pari_sp av = avma;
1120 : GEN R, C, U, P, C1, U1, v, c, d;
1121 918328 : long r, i, q, n = lg(A)-1, m;
1122 : ulong D;
1123 918328 : ulong pi = get_Fl_red(p);
1124 918328 : if (n == 0) return cgetg(1, t_MAT);
1125 918328 : r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
1126 918326 : m = nbrows(A);
1127 918326 : if (r == n)
1128 : {
1129 918312 : GEN X = Flm_gauss_from_CUP(matid_Flm(m), R, C, U, P, p, pi, &D);
1130 918320 : return gerepileupto(av, Flm_Fl_mul_pre(X, D, p, pi));
1131 : }
1132 14 : if (r < n-1) return zero_Flm(n, m);
1133 28 : for (q = n, i = 1; i < n; i++)
1134 14 : if (R[i] != i) { q = i; break; }
1135 14 : C1 = matslice(C, 1, q-1, 1, q-1);
1136 14 : c = vecslice(Flm_row(C, q), 1, q-1);
1137 14 : c = Flm_lsolve_lower_unit_pre(C1, Flm_transpose(mkmat(c)), p, pi);
1138 14 : d = cgetg(m+1, t_VECSMALL);
1139 28 : for (i=1; i<q; i++) uel(d,i) = ucoeff(c,1,i);
1140 14 : uel(d,q) = p-1;
1141 21 : for (i=q+1; i<=m; i++) uel(d,i) = 0;
1142 14 : U1 = vecslice(U, 1, n-1);
1143 14 : v = gel(Flm_rsolve_upper_pre(U1, mkmat(gel(U,n)), p, pi),1);
1144 14 : v = vecsmall_append(v, p-1);
1145 14 : D = perm_sign(P) != (odd(q+n)?-1:1) ? p-1 : 1;
1146 28 : for (i = 1; i <= n-1; i++)
1147 14 : D = Fl_mul_pre(D, ucoeff(U1, i, i), p, pi);
1148 14 : d = Flv_Fl_mul(d, D, p);
1149 14 : return rowpermute(Flc_Flv_mul(v, d, p),perm_inv(P));
1150 : }
1151 :
1152 : static GEN
1153 287 : Flm_invimage_CUP(GEN A, GEN B, ulong p) {
1154 287 : pari_sp av = avma;
1155 : GEN R, Rc, C, U, P, B1, B2, C1, C2, X, Y, Z;
1156 : long r;
1157 287 : ulong pi = get_Fl_red(p);
1158 287 : r = Flm_CUP_pre(A, &R, &C, &U, &P, p, pi);
1159 287 : Rc = indexcompl(R, nbrows(B));
1160 287 : C1 = rowpermute(C, R);
1161 287 : C2 = rowpermute(C, Rc);
1162 287 : B1 = rowpermute(B, R);
1163 287 : B2 = rowpermute(B, Rc);
1164 287 : Z = Flm_rsolve_lower_unit_pre(C1, B1, p, pi);
1165 287 : if (!gequal(Flm_mul_pre(C2, Z, p, pi), B2))
1166 14 : return NULL;
1167 273 : Y = vconcat(Flm_rsolve_upper_pre(vecslice(U, 1, r), Z, p, pi),
1168 273 : zero_Flm(lg(A) - 1 - r, lg(B) - 1));
1169 273 : X = rowpermute(Y, perm_inv(P));
1170 273 : return gerepileupto(av, X);
1171 : }
1172 :
1173 : GEN
1174 931 : Flm_invimage_i(GEN A, GEN B, ulong p)
1175 : {
1176 : GEN d, x, X, Y;
1177 931 : long i, j, nY, nA = lg(A)-1, nB = lg(B)-1;
1178 :
1179 931 : if (!nB) return cgetg(1, t_MAT);
1180 784 : if (maxss(nA, nB) >= Flm_CUP_LIMIT && nbrows(B) >= Flm_CUP_LIMIT)
1181 287 : return Flm_invimage_CUP(A, B, p);
1182 :
1183 497 : x = Flm_ker(shallowconcat(Flm_neg(A,p), B), p);
1184 : /* AX = BY, Y in strict upper echelon form with pivots = 1.
1185 : * We must find T such that Y T = Id_nB then X T = Z. This exists iff
1186 : * Y has at least nB columns and full rank */
1187 497 : nY = lg(x)-1;
1188 497 : if (nY < nB) return NULL;
1189 497 : Y = rowslice(x, nA+1, nA+nB); /* nB rows */
1190 497 : d = cgetg(nB+1, t_VECSMALL);
1191 1778 : for (i = nB, j = nY; i >= 1; i--, j--)
1192 : {
1193 1295 : for (; j>=1; j--)
1194 1288 : if (coeff(Y,i,j)) { d[i] = j; break; }
1195 1288 : if (!j) return NULL;
1196 : }
1197 : /* reduce to the case Y square, upper triangular with 1s on diagonal */
1198 490 : Y = vecpermute(Y, d);
1199 490 : x = vecpermute(x, d);
1200 490 : X = rowslice(x, 1, nA);
1201 490 : return Flm_mul(X, Flm_inv_upper_1(Y,p), p);
1202 : }
1203 : GEN
1204 889 : Flm_invimage(GEN A, GEN B, ulong p)
1205 : {
1206 889 : pari_sp av = avma;
1207 889 : GEN X = Flm_invimage_i(A,B,p);
1208 889 : if (!X) return gc_NULL(av);
1209 889 : return gerepileupto(av, X);
1210 : }
1211 :
1212 : GEN
1213 131616 : Flm_Flc_invimage(GEN A, GEN y, ulong p)
1214 : {
1215 131616 : pari_sp av = avma;
1216 131616 : long i, l = lg(A);
1217 : GEN M, x;
1218 : ulong t;
1219 :
1220 131616 : if (l==1) return NULL;
1221 131476 : if (lg(y) != lgcols(A)) pari_err_DIM("Flm_Flc_invimage");
1222 131476 : M = cgetg(l+1,t_MAT);
1223 1080350 : for (i=1; i<l; i++) gel(M,i) = gel(A,i);
1224 131476 : gel(M,l) = y; M = Flm_ker(M,p);
1225 131476 : i = lg(M)-1; if (!i) return gc_NULL(av);
1226 :
1227 128243 : x = gel(M,i); t = x[l];
1228 128243 : if (!t) return gc_NULL(av);
1229 :
1230 128236 : setlg(x,l); t = Fl_inv(Fl_neg(t,p),p);
1231 128236 : if (t!=1) x = Flv_Fl_mul(x, t, p);
1232 128236 : return gerepileuptoleaf(av, x);
1233 : }
|
