pari-2.13.1.tar.gz: pari-2.13.1/src/basemath/Flv.c

"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/Flv.c" (4 Jan 2021, 32119 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Flv.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report:

2.13.0_vs_2.13.1.

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 Flv_to_ZV(GEN x) 20 { pari_APPLY_type(t_VEC, utoi(x[i])) } 21 22 GEN 23 Flc_to_ZC(GEN x) 24 { pari_APPLY_type(t_COL, utoi(x[i])) } 25 26 GEN 27 Flm_to_ZM(GEN x) 28 { pari_APPLY_type(t_MAT, Flc_to_ZC(gel(x,i))) } 29 30 GEN 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 */ 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 */ 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)) */ 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)) */ 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 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 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 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 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 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 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 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 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 652 Flm_ker_sp(GEN x, ulong p, long deplin) { 653 return Flm_ker_i(x, p, deplin, 1); 654 } 655 656 GEN 657 Flm_ker(GEN x, ulong p) { 658 return Flm_ker_i(x, p, 0, 0); 659 } 660 661 GEN 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 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 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 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 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 781 Flm_det_sp(GEN x, ulong p) { return Flm_det_i(x, p, 1); } 782 ulong 783 Flm_det(GEN x, ulong p) { return Flm_det_i(x, p, 0); } 784 785 /* Destroy x */ 786 static GEN 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 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 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 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 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 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 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 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 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 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 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 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 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 1101 Flm_inv_sp(GEN a, ulong *detp, ulong p) { 1102 return Flm_inv_i(a, detp, p, 1); 1103 } 1104 1105 GEN 1106 Flm_inv(GEN a, ulong p) { 1107 return Flm_inv_i(a, NULL, p, 0); 1108 } 1109 1110 GEN 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 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 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 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 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 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 }