I am trying to calculate the Hermite Normal Form of a square $n \times n$ matrix using the Extended Euclidean Algorithm to compute the columns of the HNF matrix, rather than the standard (column) matrix operations. According to Cohen [Coh93], this is possible using the following algorithm, where $a_{i,j}$ denotes the coefficients of a matrix $A$, $A_{j}$ denotes column $j$ of the matrix $A$, and $B$ is an auxiliary column vector.
$ \def\indent{\qquad} k \leftarrow n\\ j \leftarrow n\\ \text{for } i \leftarrow n \text{ to } 1 \text{ do}\\ \text{ compute } (u, v, d) \text{ such that } ua_{i,k} + va_{i,j} = d = \gcd(a_{i,k},a_{i,j}) \text{ with } |u| \text{ and } |v| \text{ minimal (see below)}\\ \text{ while } j \geq 1 \text{ do}\\ \indent\text{ if } a_{i,j} \neq 0 \text{ then}\\ \indent\indent B \leftarrow uA_{k} + vA_{j}\\ \indent\indent A_{j} \leftarrow (a_{i,k}/d)A_{j} - (a_{i,j}/d)A_{k}\\ \indent\indent A_{k} \leftarrow B\\ \indent\text{ end if}\\ \indent\text{ } j \leftarrow j - 1\\ \text{ end while}\\ \text{end for}\\ \vspace{5mm}\\ \text{// Final Reductions}\\ k \leftarrow n\\ \text{for } i \leftarrow n \text{ to } 1 \text{ do}\\ \text{ } b \leftarrow a_{i,k}\\ \text{ if } b < 0 \text{ then}\\ \indent\text{ } A_{k} \leftarrow -A_{k}\\ \indent\text{ } b \leftarrow -b\\ \text{ end if}\\ \text{ if } b = 0 \text{ then}\\ \indent\text{ continue}\\ \text{ end if}\\ \text{ for } j \leftarrow k+1 \text{ to } n \text{ do}\\ \indent\text{ } q \leftarrow \lfloor a_{i,j} / b \rfloor\\ \indent\text{ } A_{j} \leftarrow A_{j} - qA_{k}\\ \text{ end for}\\ \text{ } k \leftarrow k - 1\\ \text{end for}\\ $
Note that the Euclidean step requires that $|u|$ and $|v|$ be minimal. This means that of all the possible unique pairs, $(u,v)$ be chosen such that
$ -\frac{|a|}{d} < v \text{sign}(b) \leq 0 \text{ and } 1 \leq u \text{sign}(a) \leq \frac{|b|}{d} $
My algorithm for calculating the minimum $|u|$ and $|v|$ therefore consists of the following:
$ d, u, v = \text{XGCD}(a_{i,k},a_{i,j}); \text{// Using XGCD defined by Magma} \\ x_{0} \leftarrow u \\ y_{0} \leftarrow v \\ \text{if } (u = 0) \text{ or } (v = 0) then \\ \text{ for } i \leftarrow 0 \text{ to } d \text{ do} \\ \indent\text{ } u \leftarrow x_{0} - i \times b / d \\ \indent\text{ } v \leftarrow y_{0} + i \times a / d \\ \indent\text{ if } -\frac{|a|}{d} < v \leq 0 \text{ and } 1 \leq u \leq \frac{|b|}{d} \text{ then} \\ \indent\indent\text{ break} \\ \indent\text{ end if} \\ \text{ end for} \\ \text{end if} \\ \text{return } d, u, v \\ $
Here, I have assumed that Magma's XGCD returns the correct sign for $u$ and $v$.
To test this algorithm, I am using the following sample matrix as input to the procedure:
$ \left( \begin{array}{rrrr} 10 & 0 & 1 & 0 \\ -2 & -2 & 2 & 1 \\ -2 & -1 & 3 & -1 \\ 4 & -6 & -2 & 0 \\ \end{array} \right) $
I expect the following output from the procedure:
$ \left( \begin{array}{rrrr} 33 & 12 & 12 & 11 \\ 0 & 6 & 5 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 2 \\ \end{array} \right) $
Instead, I am getting the following output:
$ \left( \begin{array}{rrrr} 3168 & 48 & 48 & 47 \\ 0 & 24 & 22 & 18 \\ 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 2 \\ \end{array} \right) $
It is apparent that there is no problem with the steps that reduce the lower triangular part of the matrix, nor with the steps that reduce the upper triangular part of the matrix (that is, the Final Reduction steps). The problem therefore lies with the Euclidean step.
My question is whether or not I have correctly implemented the procedure that returns the minimal values for $|u|$ and $|v|$? Note that I have used the following link relating to Linear Diophantine Equations as a guideline to implementing this procedure: Euclidean Algorithm. Is this correct, or should I be tackling this problem a different way? Either way, where could I find a better reference that better describes the problem at hand? Any help/advice would be appreciated.
References:
[Coh93] H. Cohen, A Course in Computational Number Theory, Vol. 138 of Graduate Texts in Mathematics. Springer, Heidelberg, 1993.