number theory - 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?

Any help/advice/feedback would be appreciated.

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