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