# Computing Hermite Normal Form using Extended Euclidean Algorithm (modulo $D$)

I am attempting to compute the Hermite Normal Form of a matrix $A$, again, only this time using the Extended Euclidean Algorithm modulo the determinant of the matrix. This was a method contributed by Hafner and McCurley [HM91], a variation of which is provided below (ref. [Coh93]).

$b \leftarrow \det(A)\\ \text{for } i \leftarrow n \text{ to } 1 \text{ do}\\ \quad \text{compute } (d, u, v) \text{ such that } ua_{i,i} + vb = d = \gcd(a_{i,i},b)\\ \quad H_{i} \leftarrow (uA_{i} \mod |b|)\\ \quad \text{if } d = |b| \text{ then}\\ \quad\quad h_{i,i} \leftarrow d\\ \quad \text{end if}\\ \quad b \leftarrow b / d\\ \text{end for}\\$

$\text{for } i \leftarrow (n - 1) \text{ to } 1 \text{ do}\\ \quad\text{for } j \leftarrow (i + 1) \text{ to } n \text{ do}\\ \quad\quad q \leftarrow \lfloor h_{i,j} / h_{i,i} \rfloor\\ \quad\quad H_{j} \leftarrow H_{j} - q \times H_{i}\\ \quad\text{end for}\\ \text{end for}$

Trying this for the following simple $2 \times 2$ square matrix,

$\left( \begin{array}{rr} 4 & 2 \\ 3 & 1 \\ \end{array} \right)$

the output that I expect is,

$\left( \begin{array}{rr} 2 & 0 \\ 0 & 1 \\ \end{array} \right).$

On the first iteration of the algorithm, we have $b = \det(A) = -2$, and

$\begin{array}{rcl} d, u, v & = & \text{XGCD}(a_{i,i}, b) \\ & = & \text{XGCD}(1,-2) \\ & = & 1, 1, 0 \\ \end{array}$

The last column of the HNF matrix, $H$, is therefore calculated thus,

$\begin{array}{rcl} H_{i} & = & uA_{i} \mod |b| \\ & = & 1 \times \left( \begin{array}{c} 2 \\ 1 \\ \end{array} \right) \mod 2 \\ & = & \left( \begin{array}{c} 0 \\ 1 \\ \end{array} \right) \\ \end{array}$

On the second iteration of the algorithm, we have $b = b / d = -2 / 1 = -2$, and

$\begin{array}{rcl} d, u, v & = & \text{XGCD}(a_{i,i}, b) \\ & = & \text{XGCD}(4,-2) \\ & = & 2, 0, -1 \\ \end{array}$

The first column of the HNF matrix, $H$, is therefore calculated thus,

$\begin{array}{rcl} H_{i} & = & uA_{i} \mod |b| \\ & = & 0 \times \left( \begin{array}{c} 4 \\ 3 \\ \end{array} \right) \mod 2 \\ & = & \left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) \\ \end{array}$

But, because $d = |b|$, $h_{1,1}$ is assigned the value of $2$, therefore producing the correct output,

$\begin{array}{rcl} H_{1} & = & \left( \begin{array}{c} 2 \\ 0 \\ \end{array} \right) \\ \end{array}$

When I apply the algorithm to larger dimensions, though, I do not get the result that I expect. For example, providing the following matrices as input,

$A = \left( \begin{array}{rrrr} 10 & 0 & 1 & 0 \\ -2 & -2 & 2 & 1 \\ -2 & -1 & 3 & -1 \\ 4 & -6 & -2 & 0 \\ \end{array} \right), \text{ } B = \left( \begin{array}{rrrr} -8 & 5 & -3 & 1 \\ 0 & 0 & -3 & -4 \\ 7 & -3 & -5 & 2 \\ -6 & -6 & 2 & 1 \\ \end{array} \right)$

results in the following matrices as output, respectively,

$\text{HNF}(A) = \left( \begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 396 \\ \end{array} \right), \text{ } \text{HNF}(B) = \left( \begin{array}{rrrr} 1 & 0 & 2299 & 1 \\ 0 & 2873 & 2299 & 2869 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 2298 & 1 \\ \end{array} \right)$

The expected result are

$\text{HNF}(A) = \left( \begin{array}{rrrr} 33 & 12 & 12 & 11 \\ 0 & 6 & 5 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 2 \\ \end{array} \right), \text{ } \text{HNF}(B) = \left( \begin{array}{rrrr} 2873 & 1663 & 286 & 335 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)$

respectively.

I am certain it is an implementation issue, however, the question I would have to ask is whether a coefficient with the value of zero on the diagonal will make any difference to the algorithm, especially since no condition exists to cater for it like in the previous (related) question I asked concerning the computation of hermite normal forms using the extended euclidean algorithm.

Furthermore, how can the lower triangular portion of the matrix be reduced if there is no step like $A_{j} \leftarrow (a_{i,k}/d)A_{j} - (a_{i,j}/d)A_{k}$ to do this?

References:

[Coh93] A Course in Computational Algebraic Number Theory, Vol. 138 of Graduate Texts in Mathematics. Springer, Heidelberg, 1993.

[HM91] Asymptotically Fast Triangularization of Matrices over Rings. SIAM Journal of Computing, 20:1068-1083, December 1991

-
This definitely ain't geometry, dude... –  Guess who it is. Mar 29 '12 at 5:53
Couldn't find a label that was more suitable to Euclidean Algorithm, and I don't have enough points to create a new label. Thanks for the correction. –  Bill Mar 29 '12 at 6:01

I have done some further research and discovered that the above algorithm is to be used as a post-process to a modified version of the algorithm mentioned in the previous question I raised concerning the computation of Hermite normal forms using the extended Euclidean algorithm. The modification entails reducing the columns of the matrix modulo the determinant of the input matrix $A$ as is illustrated below.

$D \leftarrow \det(A)\\ k \leftarrow n\\ \text{for } i \leftarrow n \text{ to } 1 \text{ do}\\ \quad\text{compute } (u, v, d) \text{ such that } ua_{i,k} + va_{i,j} = d = \gcd(a_{i,k},a_{i,j})\\ \quad\text{for } j \leftarrow k \text{ to } 1 \text{ do}\\ \qquad\text{if } a_{i,j} \neq 0 \text{ then}\\ \qquad\quad B \leftarrow (uA_{k} + vA_{j}) \mod |D|\\ \qquad\quad A_{j} \leftarrow ((a_{i,k}/d)A_{j} - (a_{i,j}/d)A_{k}) \mod |D|\\ \qquad\quad A_{k} \leftarrow B \mod |D|\\ \qquad\text{end if}\\ \quad\text{end for}\\ \quad k \leftarrow k - 1\\ \text{end for}\\$

$\text{// Final Reductions}\\ k \leftarrow n\\ \text{for } i \leftarrow n \text{ to } 1 \text{ do}\\ \quad b \leftarrow a_{i,k}\\ \quad\text{if } b < 0 \text{ then}\\ \qquad A_{k} \leftarrow -A_{k}\\ \qquad b \leftarrow -b\\ \quad{end if}\\ \quad\text{if } b = 0 \text{ then}\\ \qquad\text{continue}\\ \quad\text{end if}\\ \quad\text{for } j \leftarrow k+1 \text{ to } n \text{ do}\\ \qquad q \leftarrow \lfloor a_{i,j} / b \rfloor\\ \qquad A_{j} \leftarrow (A_{j} - qA_{k}) \mod |D|\\ \quad\text{end for}\\ \quad k \leftarrow k - 1\\ \text{end for}\\$

Inputting $A$ into this modified version of the algorithm produces the correct HNF matrix.

$\text{HNF}(A) = \left( \begin{array}{rrrr} 33 & 12 & 12 & 11 \\ 0 & 6 & 5 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 2 \\ \end{array} \right)$

Inputting $B$, however, results in the following matrix, which is not quite the HNF of B.

$\text{HNF}(B) = \left( \begin{array}{rrrr} 0 & 2116 & 739 & 335 \\ 0 & 3 & 2 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)$

This is where the the above algorithm documented in this question is used. It transforms the above (intermediate) matrix into the proper HNF matrix for B.

$\text{HNF}(B) = \left( \begin{array}{rrrr} 2873 & 1663 & 286 & 335 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)$

-