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
14 void
fmpz_mat_hnf_classical(fmpz_mat_t H,const fmpz_mat_t A)15 fmpz_mat_hnf_classical(fmpz_mat_t H, const fmpz_mat_t A)
16 {
17 slong i, i0, j, j2, k, l, m, n;
18 fmpz_t min, q;
19
20 m = fmpz_mat_nrows(A);
21 n = fmpz_mat_ncols(A);
22
23 fmpz_init(q);
24 fmpz_mat_set(H, A);
25
26 for (j = 0, k = 0, l = (n - m) * (n > m); n - j != l; j++, k++)
27 {
28 int col_finished = 1;
29 for (i = k + 1; (i < m) && col_finished; i++)
30 col_finished = fmpz_is_zero(fmpz_mat_entry(H, i, j));
31 if (col_finished)
32 {
33 if (fmpz_sgn(fmpz_mat_entry(H, k, j)) < 0)
34 {
35 for (j2 = j; j2 < n; j2++)
36 fmpz_neg(fmpz_mat_entry(H, k, j2),
37 fmpz_mat_entry(H, k, j2));
38 }
39 if (fmpz_is_zero(fmpz_mat_entry(H, k, j)))
40 {
41 k--;
42 if (l > 0)
43 l--;
44 }
45 else
46 {
47 /* reduce first entries of column j with row k */
48 for (i = 0; i < k; i++)
49 {
50 fmpz_fdiv_q(q, fmpz_mat_entry(H, i, j),
51 fmpz_mat_entry(H, k, j));
52 for (j2 = j; j2 < n; j2++)
53 {
54 fmpz_submul(fmpz_mat_entry(H, i, j2), q,
55 fmpz_mat_entry(H, k, j2));
56 }
57 }
58 }
59 }
60 else
61 {
62 i0 = 0;
63 fmpz_init(min);
64 /* locate non-zero entry in column j below k with lowest absolute
65 value */
66 for (i = k + 1; i < m; i++)
67 {
68 if (fmpz_is_zero(fmpz_mat_entry(H, i, j)))
69 continue;
70 if (fmpz_is_zero(min) ||
71 fmpz_cmpabs(fmpz_mat_entry(H, i, j), min) < 0)
72 {
73 i0 = i;
74 fmpz_abs(min, fmpz_mat_entry(H, i, j));
75 }
76 }
77 /* move the row found to row k */
78 if (i0 > k)
79 fmpz_mat_swap_rows(H, NULL, i0, k);
80
81 if (fmpz_sgn(fmpz_mat_entry(H, k, j)) < 0)
82 {
83 for (j2 = j; j2 < n; j2++)
84 {
85 fmpz_neg(fmpz_mat_entry(H, k, j2),
86 fmpz_mat_entry(H, k, j2));
87 }
88 }
89 /* reduce lower entries of column j with row k */
90 for (i = k + 1; i < m; i++)
91 {
92 fmpz_fdiv_q(q, fmpz_mat_entry(H, i, j),
93 fmpz_mat_entry(H, k, j));
94 for (j2 = j; j2 < n; j2++)
95 {
96 fmpz_submul(fmpz_mat_entry(H, i, j2), q,
97 fmpz_mat_entry(H, k, j2));
98 }
99 }
100 /* don't move to the next column yet */
101 j--;
102 k--;
103 fmpz_clear(min);
104 }
105 }
106
107 fmpz_clear(q);
108 }
109