1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/*
Copyright (C) 2014 Alex J. Best
Copyright (C) 2026 Edgar Costa
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include "fmpz.h"
#include "fmpz_mat.h"
/*
Compute SNF via iterative Hermite normal form.
Phase 1: alternately compute row-HNF and column-HNF (via transpose)
until the matrix is diagonal.
Phase 2: fix divisibility chain on the diagonal using gcd/lcm.
Phase 3: make diagonal entries non-negative.
This is the same algorithm as fmpz_mat_snf_transform but without
tracking the transformation matrices U and V.
*/
static void
_fmpz_mat_snf_iterative_hermite(fmpz_mat_t S, const fmpz_mat_t A)
{
slong m = fmpz_mat_nrows(A);
slong n = fmpz_mat_ncols(A);
slong d = FLINT_MIN(m, n);
slong j, k;
slong max_iters, iter;
fmpz_mat_t X, Xt;
fmpz_t dd, pp, qq;
if (d == 0)
{
fmpz_mat_zero(S);
return;
}
fmpz_mat_init(X, m, n);
fmpz_mat_init(Xt, n, m);
fmpz_init(dd);
fmpz_init(pp);
fmpz_init(qq);
fmpz_mat_set(X, A);
max_iters = _fmpz_mat_snf_iter_bound(A);
for (iter = 0; !fmpz_mat_is_diagonal(X); iter++)
{
if (iter >= max_iters)
flint_throw(FLINT_ERROR,
"(fmpz_mat_snf): Phase 1 exceeded iteration bound "
"(%wd). Likely bug in fmpz_mat_hnf or unexpected "
"input; please report.\n", max_iters);
/* Row HNF */
fmpz_mat_hnf(X, X);
if (fmpz_mat_is_diagonal(X))
break;
/* Column HNF via transpose */
fmpz_mat_transpose(Xt, X);
fmpz_mat_hnf(Xt, Xt);
fmpz_mat_transpose(X, Xt);
}
/*
Phase 2: fix divisibility chain on diagonal.
For each pair (j, k) with j < k, if X[j,j] does not divide X[k,k],
replace with gcd and lcm.
*/
for (j = 0; j < d; j++)
{
if (fmpz_is_one(fmpz_mat_entry(X, j, j)))
continue;
for (k = j + 1; k < d; k++)
{
if (fmpz_is_zero(fmpz_mat_entry(X, k, k)))
continue;
if (!fmpz_is_zero(fmpz_mat_entry(X, j, j))
&& fmpz_divisible(fmpz_mat_entry(X, k, k),
fmpz_mat_entry(X, j, j)))
continue;
fmpz_gcd(dd, fmpz_mat_entry(X, j, j),
fmpz_mat_entry(X, k, k));
fmpz_divexact(pp, fmpz_mat_entry(X, j, j), dd);
fmpz_divexact(qq, fmpz_mat_entry(X, k, k), dd);
/* X[j,j] = gcd, X[k,k] = lcm = p*q*d */
fmpz_set(fmpz_mat_entry(X, j, j), dd);
fmpz_mul(fmpz_mat_entry(X, k, k), pp, qq);
fmpz_mul(fmpz_mat_entry(X, k, k),
fmpz_mat_entry(X, k, k), dd);
}
}
/*
Phase 3: make diagonal entries non-negative.
*/
for (j = 0; j < d; j++)
{
if (fmpz_sgn(fmpz_mat_entry(X, j, j)) < 0)
fmpz_neg(fmpz_mat_entry(X, j, j),
fmpz_mat_entry(X, j, j));
}
fmpz_mat_set(S, X);
fmpz_mat_clear(Xt);
fmpz_mat_clear(X);
fmpz_clear(dd);
fmpz_clear(pp);
fmpz_clear(qq);
}
void
fmpz_mat_snf(fmpz_mat_t S, const fmpz_mat_t A)
{
if (fmpz_mat_is_diagonal(A))
{
fmpz_mat_snf_diagonal(S, A);
return;
}
_fmpz_mat_snf_iterative_hermite(S, A);
}