1 /*
2 Copyright (C) 2011 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "flint.h"
13 #include "nmod_poly.h"
14 #include "nmod_poly_mat.h"
15
16 #define E(j,k) nmod_poly_mat_entry(B,j,k)
17
18 static __inline__ void
nmod_poly_mat_swap_rows(nmod_poly_mat_t mat,slong * perm,slong r,slong s)19 nmod_poly_mat_swap_rows(nmod_poly_mat_t mat, slong * perm, slong r, slong s)
20 {
21 if (r != s)
22 {
23 nmod_poly_struct * u;
24 slong t;
25
26 if (perm)
27 {
28 t = perm[s];
29 perm[s] = perm[r];
30 perm[r] = t;
31 }
32
33 u = mat->rows[s];
34 mat->rows[s] = mat->rows[r];
35 mat->rows[r] = u;
36 }
37 }
38
39 slong
nmod_poly_mat_fflu(nmod_poly_mat_t B,nmod_poly_t den,slong * perm,const nmod_poly_mat_t A,int rank_check)40 nmod_poly_mat_fflu(nmod_poly_mat_t B, nmod_poly_t den, slong * perm,
41 const nmod_poly_mat_t A, int rank_check)
42 {
43 nmod_poly_t t;
44 slong m, n, j, k, rank, r, pivot_row, pivot_col;
45
46 if (nmod_poly_mat_is_empty(A))
47 {
48 nmod_poly_one(den);
49 return 0;
50 }
51
52 nmod_poly_mat_set(B, A);
53 m = B->r;
54 n = B->c;
55 rank = pivot_row = pivot_col = 0;
56
57 nmod_poly_init(t, nmod_poly_mat_modulus(A));
58
59 while (pivot_row < m && pivot_col < n)
60 {
61 r = nmod_poly_mat_find_pivot_partial(B, pivot_row, m, pivot_col);
62
63 if (r == -1)
64 {
65 if (rank_check)
66 {
67 nmod_poly_zero(den);
68 rank = 0;
69 break;
70 }
71 pivot_col++;
72 continue;
73 }
74 else if (r != pivot_row)
75 nmod_poly_mat_swap_rows(B, perm, pivot_row, r);
76
77 rank++;
78
79 for (j = pivot_row + 1; j < m; j++)
80 {
81 for (k = pivot_col + 1; k < n; k++)
82 {
83 nmod_poly_mul(E(j, k), E(j, k), E(pivot_row, pivot_col));
84 nmod_poly_mul(t, E(j, pivot_col), E(pivot_row, k));
85 nmod_poly_sub(E(j, k), E(j, k), t);
86 if (pivot_row > 0)
87 nmod_poly_div(E(j, k), E(j, k), den);
88 }
89 }
90
91 nmod_poly_set(den, E(pivot_row, pivot_col));
92 pivot_row++;
93 pivot_col++;
94 }
95
96 nmod_poly_clear(t);
97 return rank;
98 }
99