1 /*
2 Copyright (C) 2014 Alex J. Best
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 "fmpz_mat.h"
13
_eliminate_col(fmpz_mat_t S,slong i,const fmpz_t mod)14 static void _eliminate_col(fmpz_mat_t S, slong i, const fmpz_t mod)
15 {
16 slong j, k, m, n;
17 fmpz * t;
18 fmpz_t b, g, u, v, r1g, r2g;
19
20 m = S->r;
21 n = S->c;
22
23 if (i == m - 1)
24 {
25 fmpz_gcd(fmpz_mat_entry(S, i, i), fmpz_mat_entry(S, i, i), mod);
26 return;
27 }
28
29 fmpz_init(g);
30 fmpz_init(u);
31 fmpz_init(b);
32 fmpz_init(r1g);
33 fmpz_init(r2g);
34
35 if (!fmpz_is_zero(fmpz_mat_entry(S, i, i)))
36 {
37 fmpz_init(v);
38
39 fmpz_xgcd(g, u, v, fmpz_mat_entry(S, i + 1, i),
40 fmpz_mat_entry(S, i, i));
41 fmpz_divexact(r1g, fmpz_mat_entry(S, i + 1, i), g);
42 fmpz_divexact(r2g, fmpz_mat_entry(S, i, i), g);
43 for (j = i; j < n; j++)
44 {
45 fmpz_mul(b, u, fmpz_mat_entry(S, i + 1, j));
46 fmpz_addmul(b, v, fmpz_mat_entry(S, i, j));
47 fmpz_mul(fmpz_mat_entry(S, i, j), r1g,
48 fmpz_mat_entry(S, i, j));
49 fmpz_submul(fmpz_mat_entry(S, i, j), r2g,
50 fmpz_mat_entry(S, i + 1, j));
51 fmpz_set(fmpz_mat_entry(S, i + 1, j), b);
52 }
53
54 fmpz_clear(v);
55 }
56
57 /* compute extended gcd of entries in column i */
58 t = _fmpz_vec_init(m - i - 1);
59
60 fmpz_set(g, fmpz_mat_entry(S, i + 1, i));
61 fmpz_one(t);
62 for (j = 2; j < m - i; j++)
63 {
64 fmpz_xgcd(g, u, t + j - 1, g, fmpz_mat_entry(S, i + j, i));
65 for (k = 0; k < j - 1; k++)
66 fmpz_mul(t + k, t + k, u);
67 }
68
69 /* set row i to have gcd in col i */
70 for (k = i + 1; k < m; k++)
71 {
72 fmpz_mod(t + k - i - 1, t + k - i - 1, mod);
73 for (j = i; j < n; j++)
74 fmpz_addmul(fmpz_mat_entry(S, i, j), t + k - i - 1,
75 fmpz_mat_entry(S, k, j));
76 }
77
78 _fmpz_vec_clear(t, m - i - 1);
79
80 /* reduce each row k with row i */
81 if (!fmpz_is_zero(g)) /* if g = 0 then don't need to reduce */
82 {
83 for (k = i + 1; k < m; k++)
84 {
85 fmpz_divexact(r1g, fmpz_mat_entry(S, k, i), g);
86 fmpz_neg(r1g, r1g);
87 for (j = i; j < n; j++)
88 fmpz_addmul(fmpz_mat_entry(S, k, j), r1g,
89 fmpz_mat_entry(S, i, j));
90 }
91 for (k = i + 1; k < m; k++)
92 fmpz_mod(fmpz_mat_entry(S, k, i), fmpz_mat_entry(S, k, i), mod);
93 }
94 for (j = i; j < m; j++)
95 for (k = i + 1; k < n; k++)
96 fmpz_fdiv_r(fmpz_mat_entry(S, j, k), fmpz_mat_entry(S, j, k), mod);
97 fmpz_gcd(fmpz_mat_entry(S, i, i), fmpz_mat_entry(S, i, i), mod);
98
99 fmpz_clear(b);
100 fmpz_clear(g);
101 fmpz_clear(u);
102 fmpz_clear(r1g);
103 fmpz_clear(r2g);
104 }
105
_eliminate_row(fmpz_mat_t S,slong i,const fmpz_t mod)106 static void _eliminate_row(fmpz_mat_t S, slong i, const fmpz_t mod)
107 {
108 slong j, k, m, n;
109 fmpz * t;
110 fmpz_t b, g, u, v, r1g, r2g, halfmod;
111
112 m = S->r;
113 n = S->c;
114
115 if (i == n - 1)
116 {
117 fmpz_gcd(fmpz_mat_entry(S, i, i), fmpz_mat_entry(S, i, i), mod);
118 return;
119 }
120
121 fmpz_init(g);
122 fmpz_init(u);
123 fmpz_init(b);
124 fmpz_init(r1g);
125 fmpz_init(r2g);
126 fmpz_init(halfmod);
127 fmpz_fdiv_q_2exp(halfmod, mod, 1);
128
129 if (!fmpz_is_zero(fmpz_mat_entry(S, i, i)))
130 {
131 fmpz_init(v);
132
133 fmpz_xgcd(g, u, v, fmpz_mat_entry(S, i, i + 1),
134 fmpz_mat_entry(S, i, i));
135 fmpz_divexact(r1g, fmpz_mat_entry(S, i, i + 1), g);
136 fmpz_divexact(r2g, fmpz_mat_entry(S, i, i), g);
137 for (j = i; j < m; j++)
138 {
139 fmpz_mul(b, u, fmpz_mat_entry(S, j, i + 1));
140 fmpz_addmul(b, v, fmpz_mat_entry(S, j, i));
141 fmpz_mul(fmpz_mat_entry(S, j, i), r1g,
142 fmpz_mat_entry(S, j, i));
143 fmpz_submul(fmpz_mat_entry(S, j, i), r2g,
144 fmpz_mat_entry(S, j, i + 1));
145 fmpz_set(fmpz_mat_entry(S, j, i + 1), b);
146 }
147
148 fmpz_clear(v);
149 }
150
151 /* compute extended gcd of entries in row i */
152 t = _fmpz_vec_init(n - i - 1);
153
154 fmpz_set(g, fmpz_mat_entry(S, i, i + 1));
155 fmpz_one(t);
156 for (j = 2; j < n - i; j++)
157 {
158 fmpz_xgcd(g, u, t + j - 1, g, fmpz_mat_entry(S, i, i + j));
159 for (k = 0; k < j - 1; k++)
160 fmpz_mul(t + k, t + k, u);
161 }
162
163 /* reduce col i to have gcd in row i */
164 for (k = i + 1; k < n; k++)
165 {
166 fmpz_mod(t + k - i - 1, t + k - i - 1, mod);
167 for (j = i; j < m; j++)
168 fmpz_addmul(fmpz_mat_entry(S, j, i), t + k - i - 1,
169 fmpz_mat_entry(S, j, k));
170 }
171
172 _fmpz_vec_clear(t, n - i - 1);
173
174 /* reduce each col k with col i */
175 if (!fmpz_is_zero(g)) /* if g = 0 then don't need to reduce */
176 {
177 for (k = i + 1; k < n; k++)
178 {
179 fmpz_divexact(r1g, fmpz_mat_entry(S, i, k), g);
180 fmpz_neg(r1g, r1g);
181 for (j = i; j < m; j++)
182 fmpz_addmul(fmpz_mat_entry(S, j, k), r1g,
183 fmpz_mat_entry(S, j, i));
184 }
185 }
186 for (j = i + 1; j < m; j++)
187 for (k = i; k < n; k++)
188 fmpz_fdiv_r(fmpz_mat_entry(S, j, k), fmpz_mat_entry(S, j, k), mod);
189 fmpz_gcd(fmpz_mat_entry(S, i, i), fmpz_mat_entry(S, i, i), mod);
190
191 fmpz_clear(b);
192 fmpz_clear(g);
193 fmpz_clear(u);
194 fmpz_clear(r1g);
195 fmpz_clear(r2g);
196 fmpz_clear(halfmod);
197 }
198
fmpz_mat_snf_iliopoulos(fmpz_mat_t S,const fmpz_mat_t A,const fmpz_t mod)199 void fmpz_mat_snf_iliopoulos(fmpz_mat_t S, const fmpz_mat_t A, const fmpz_t mod)
200 {
201 slong i, k, n;
202 int done;
203
204 n = FLINT_MIN(A->c, A->r);
205
206 fmpz_mat_set(S, A);
207
208 for (i = 0; i < A->r; i++)
209 for (k = 0; k < A->c; k++)
210 fmpz_mod(fmpz_mat_entry(S, i, k), fmpz_mat_entry(S, i, k), mod);
211
212 for (k = 0; k != n; k++)
213 {
214 do
215 {
216 _eliminate_row(S, k, mod);
217 _eliminate_col(S, k, mod);
218 done = 1;
219 if (fmpz_is_zero(fmpz_mat_entry(S, k, k)))
220 {
221 for (i = k + 1; i < A->c && done; i++)
222 done = fmpz_is_zero(fmpz_mat_entry(S, k, i));
223 }
224 else
225 {
226 for (i = k + 1; i < A->c && done; i++)
227 done = fmpz_divisible(fmpz_mat_entry(S, k, i),
228 fmpz_mat_entry(S, k, k));
229 }
230 }
231 while (!done);
232 for (i = k + 1; i < A->c; i++)
233 fmpz_zero(fmpz_mat_entry(S, k, i));
234 }
235
236 fmpz_mat_snf_diagonal(S, S);
237 }
238