I think there are some errors, either book errata, or transcription errors, in the stated algorithm. Firstly, $k$ is the pivot column, so it needs to be decremented (whenever nonzero) at the end of the first $i$ loop. Secondly, for each iteration in this loop, $j$ needs to be looped over from $k-1$ down to $1$, to clear out the lower triangular part, so it must be initialized inside the $i$ loop.
I am also skeptical of whether special handling of xgcd is necessary or even helpful, since it always returns a "positive definite" $d$ (an admitted abuse of terminology), in the sense that it is (or should be) always positive given $(a,b)\ne(0,0)$. Because of this, we also do not need the $b < 0$ test (and correction) in the final (partial upper triangular) reductions. In fact, it seems rather artificial to require $v \le 0$ and $u > 0$ of xgcd. What about $(a,b)=(3,1)$ and $(-1,3)$? Also, your pseudocode for the xgcd wrapper only alters its output if $u$ or $v$ is zero, and I think the loop could be avoided by solving some linear equations (taking quotients and a minimum).
However, I don't have Cohen's book, and don't know how or why he requires $|u|$ and $|v|$ to be minimal. Does he mean that $\max\{|u|,|v|\}$ must be minimal? I admit, I haven't thought deeply about it. But I think it is a fact that xgcd already does a good job of producing a minimal solution.
I have corrected the algorithm, as I think it should be, and implemented it in sage. After testing described below, I also made a correction.
m = Matrix(ZZ,4,4,\
[[10,-2,-2, 4]\
,[ 0,-2,-1,-6]\
,[ 1, 2, 3,-2]\
,[ 0, 1,-1, 0]]).transpose(); m
def hnf(m,show=False):
nr = m.nrows()
nc = m.ncols()
if show:
print 'hnf: '+str(nr)+'x'+str(nc)
assert(nr == nc)
c = m.columns()
M = Matrix(ZZ,nc,nr,c).transpose()
if show:
print M
def r(n):
s = range(n)
s.reverse()
return s
k = nc - 1
pivots = []
for i in r(nr):
for j in r(k):
ck = c[k]; a = ck[i]
cj = c[j]; b = cj[i]
if b != 0:
d,u,v = xgcd(a,b)
if show:
print 'row '+str(i+1)+', cols '+str(j+1)+','+str(k+1)\
, 'xgcd:',u,'(',a,') + ',v,'(',b, ') = ',d
ad = a/d
bd = b/d
c[k] = ck*u + cj*v # if a != d:
c[j] = cj*ad - ck*bd
M = Matrix(ZZ,nc,nr,c).transpose()
if show:
print M
if c[k][i] != 0:
pivots.append((i,k))
if c[k][i] < 0: # this is necessary for the last pivot column,
c[k] = -c[k] # since it cannot utilize xgcd to "normalize"
k = k - 1
for (i,k) in pivots:
b = c[k][i]; assert(b > 0)
if show:
print 'pivot('+str(i+1)+','+str(k+1)+')='+str(b)
for j in [k+1..nc-1]:
q = floor(c[j][i] / b)
c[j] -= c[k]*q
M = Matrix(ZZ,nc,nr,c).transpose()
if show:
print M
pivots.reverse()
#if show:
# print 'pivots:', pivots
return M, pivots
h,p = hnf(m,True)
Here is the output (for your test case):
hnf: 4x4
[10 0 1 0]
[-2 -2 2 1]
[-2 -1 3 -1]
[ 4 -6 -2 0]
row 4, cols 3,4 xgcd: 0 ( 0 ) + -1 ( -2 ) = 2
[10 0 0 -1]
[-2 -2 1 -2]
[-2 -1 -1 -3]
[ 4 -6 0 2]
row 4, cols 2,4 xgcd: 1 ( 2 ) + 0 ( -6 ) = 2
[ 10 -3 0 -1]
[ -2 -8 1 -2]
[ -2 -10 -1 -3]
[ 4 0 0 2]
row 4, cols 1,4 xgcd: 1 ( 2 ) + 0 ( 4 ) = 2
[ 12 -3 0 -1]
[ 2 -8 1 -2]
[ 4 -10 -1 -3]
[ 0 0 0 2]
row 3, cols 2,3 xgcd: -1 ( -1 ) + 0 ( -10 ) = 1
[12 3 0 -1]
[ 2 18 -1 -2]
[ 4 0 1 -3]
[ 0 0 0 2]
row 3, cols 1,3 xgcd: 1 ( 1 ) + 0 ( 4 ) = 1
[12 3 0 -1]
[ 6 18 -1 -2]
[ 0 0 1 -3]
[ 0 0 0 2]
row 2, cols 1,2 xgcd: 0 ( 18 ) + 1 ( 6 ) = 6
[33 12 0 -1]
[ 0 6 -1 -2]
[ 0 0 1 -3]
[ 0 0 0 2]
pivot(4,4)=2
[33 12 0 -1]
[ 0 6 -1 -2]
[ 0 0 1 -3]
[ 0 0 0 2]
pivot(3,3)=1
[33 12 0 -1]
[ 0 6 -1 -5]
[ 0 0 1 0]
[ 0 0 0 2]
pivot(2,2)=6
[33 12 12 11]
[ 0 6 5 1]
[ 0 0 1 0]
[ 0 0 0 2]
pivot(1,1)=33
[33 12 12 11]
[ 0 6 5 1]
[ 0 0 1 0]
[ 0 0 0 2]
Sage also has a hermite_form function for matrices,
var('s')
for s in ['pari', 'pari0', 'pari4', 'padic', 'default', 'ntl']:
print str(s)
print m.hermite_form(algorithm=s)
but each implementation produces the matrix below
$$
\left(\matrix{
2 & 0 & 2 & 27 \\
0 & 1 & 1 & 64 \\
0 & 0 & 3 & 45 \\
0 & 0 & 0 & 66
}\right)
$$
which, as you point out in a comment, is arrived at by iterating over $i$ increasing, i.e. by pivoting from the top (left) downwards,
rather than decreasing from the bottom (right) upwards.
We can verify this directly using the transformation
$\phi:m\mapsto rm^tr=(rmr)^t$, for $r$ the (orthogonal) matrix
$$
r=
\left(\begin{array}{rrrr}
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0
\end{array}\right),
$$
which transposes $m$ by both diagonals. Our (Cohen's) HNF (which pivots from the bottom right upward) and sage's (which forms pivots from the top left downwards) can be obtained from each other by pre- and post- transforming by $\phi:$
r = Matrix(ZZ,4,4,[\
floor(1/(1+n%3))*\
(1-floor(1/(1+n%15)))\
for n in range(16)])
print r
phi = lambda x: (r*x*r).transpose()
print phi(m).hermite_form() == phi(h)
phi(phi(m).hermite_form())
By the way, if you don't mind a little excursion, the group of transformations of square matrices generated under composition by the (self-inverse) maps $\color{red}{t}:m\mapsto m^t$ (transpose) and $\color{blue}{r}:m\mapsto rm$ (row reversal) or $\color{green}{c}:m\mapsto mr$ (column reversal) is in fact isomorphic to the dihedral group $D_4$ of the square. We can visualize this group by its action on a fixed matrix, e.g. the $2\times2$ matrix $e=\left(\matrix{1&2\\3&4}\right)$.
Each transformation is represented by, or identified with, (a path from $e$ to) a vertex. The order two elements consist of the surface edges and interior diagonals, while the rotations consist of the face diagonals, and their cosets partition the vertices into two intersecting tetrahedra. This is also a commutative diagram, and thus can be used to verify such formulas as $rc=cr$, $rt=tc$, and $ct=tr$ as well as $\phi=rtc=trc$ from above.
Anyway after some testing
ranks = [0,0,0,0,0]
passed = 0
failed = 0
for n in range(10^2):
m = random_matrix(ZZ,4)
h, p = hnf(phi(m))
if phi(h) == m.hermite_form():
ranks[len(p)] += 1
passed += 1
else:
print 'Test FAIL!'
print 'm:'; print m
print 'm"=phi(m):'; print phi(m)
print 'hnf(m"):'; print h
print 'hnf(m")":'; print phi(h)
print 'm.hnf():'; print m.hermite_form()
failed += 1
break
print n+1,'tests performed:',passed,'passed,',failed,'failed.'
print 'rank counts:',ranks
I found that I needed to sometimes correct the sign of the first pivot column, and incorporated this into the hnf function above (along with making the print output optional but off by default). In any case, in sage, you can just define $\phi$ and use phi(phi(m).hermite_form()); after all, why should we add a seventh implementation (except to learn!)?
Best Answer
Your equations: $$a(m) = a_1 + m d_1$$ $$b(n) = b_1 + n d_2 $$
You want $a(m) = b(n)$ or $a(m)-b(n)=0$, so it may be written as $$(a_1-b_1) + m(d_1) +n(-d_2) = 0$$ or $$ m(d_1) +n(-d_2) = (b_1-a_1) $$
You want $n$ and $m$ minimal and that solves that. This is of course very similar to the EGCD, but that the value of $b_1 - a_1$ is the desired value intead of the value $1$. EGCD sovles it for the value of $1$ (or the gcd of $d_1$ and $d_2$).
This is actually a lattice reduction problem since you are interested in the minimum integer values for $m$ and $n$. That is an involved problem, but this is the lowest dimension and thus is relatively easy.
The method I have written in the past used a matrix form to solve it. I started with the matrix $$\pmatrix{1 & 0 & d_1 \\ 0 & 1 & -d_2}$$
which represents the equations \begin{align} (1)d_1 + (0)(-d_2) = & \phantom{-} d_1 \\ (0)d_1 + (1)(-d_2) = & -d_2 \\ \end{align}
Each row of the matrix gives valid numbers for your equation, the first element in the row is the number for $m$, the second is for $n$ and the third is the value of your equation for those $m$ and $n$. Now if you combine the rows (such as row one equals row one minus row two) then you still get valid numbers. The goal then is to find a combination resulting in the desired value of $b_1 - a_1$ in the final column.
If you use EGCD you can start with this matrix: $$\pmatrix{d_2 \over g & d_1 \over g & 0 \\ u & -v & g}$$
where $u$, $v$, and $g$ are the outputs of the EGCD (with $g$ the gcd) since EGCD gives you $ud_1 + vd_2 = g$
In your example you would have: $$\pmatrix{16 & 3 & 0 \\ -5 & -1 & 5}$$ From there you can scale the last row to get $kg = (b_1 - a_1)$ for some integer $k$, then to find the minimal values use the first row to reduce, since the zero in the first row will not affect the result.
Again, for your example, $k=13$ which gives $$\pmatrix{16 & 3 & 0 \\ -65 & -13 & 65}$$
Adding $5$ times the first row gives $$\pmatrix{15 & 2 & 65}$$
Which represents the $16$th and $3$rd terms (count starts at $0$) respectively.