Integer matrix diagonalization

George Havas^{∗} and Bohdan S. Majewski^{†} Key Centre for Software Technology Department of Computer Science The University of Queensland Queensland 4072, Australia

Abstract

We consider algorithms for computing the Smith normal form of integer matrices. Various different strategies have been proposed, primarily trying to avoid the major obstacle that occurs in such computations — explosive growth in size of intermediate entries. We present a new algorithm with excellent performance. We investigate the complexity of such computations, indicating relationships with NP-complete problems. We also describe new heuristics which perform well in practice. We present experimental evidence which shows our algorithm outperforming the previous methods.

1 Introduction

Integer matrices A and B are equivalent if there exist unimodular matrices P and Q such that PAQ = B. Matrices P and Q correspond to elementary row and column operations: negating a row (or column); adding an integer multiple of one row (or column) to another; or interchanging two rows (or columns). It follows from a result of Smith [Smi61] that for any integer matrix there exists a unique diagonal matrix, with s_{ij} = 0 for i ≠ j , such that s_{i−1i−1} divides s_{ii}. This matrix is called the Smith normal form of the given matrix and has many important applications.

Smith gave a method for computing the Smith normal form, but his approach is only of theoretical merit. Over the years different strategies have been proposed, primarily trying to avoid the major obstacle that occurs in such computations — explosive growth in size of intermediate entries. A comprehensive bibliography and a number of methods are examined in [HHR93]. Here we present a new algorithm with excellent performance.

Modular methods can be used to solve these problem, but the complexity bounds on such computations have relatively high degree when compared to integer based methods. That is why we study integer methods in this paper.

∗E-mail: havas@cs.uq.oz.au; partially supported by the Australian Research Council

†E-mail: bohdan@cs.uq.oz.au

It is often the case that, for a matrix that has quite small entries in both the initial and the Smith normal form, the intermediate entries will be spectacularly large. Frumkin [Fru76] observed that for simple Gaussian elimination a worst case bound on the length of the entries is exponential. Examples are presented in [HHR93] which show that known polynomially bounded algorithms lead to quite large entries. Consequently, a careless algorithm, while executing a polynomial number of operations may have exponential complexity, simply because of the size of the numbers it operates on. In this paper we investigate the complexity of such computations, indicating relationships with NP-complete problems. As a consequence we suggest and justify a new heuristic strategy, with a polynomially bounded number of operations, which performs well in practice.

We use the following notation. For a m × n integer matrix A we denote its i-th row by a_{i∗} and its j -th column by a_{∗j} . The absolute value of x is denoted by |x|, while det(A) stands for the determinant of A. Matrix A may be alternatively written as

A = [a_{∗1}, …, a_{∗n}] = [ a_{1∗} ; ⋮ ; a_{m∗} ] .

2 Previous methods

A straightforward algorithm for computing the Smith normal form of an integer matrix is as follows. The first stage of the reduction is to compute an equivalent form

[ \begin{matrix}d_{1} & 0 & … & 0 \\ 0 \\ ⋮ & & A′ \\ 0 \end{matrix} ]

where d_{1} divides every entry in the submatrix A′.

If A is the zero matrix we are finished. If not, choose a nonzero entry (the pivot) and move it to a_{11} by suitable row and column interchanges, and make it positive.

While there is an entry a_{1j} in the first row not divisible by a_{11}: compute a_{1j} = a_{11}q + r; subtract q times the first column from the j th column; interchange the first and j th columns. Do the same with rows and columns interchanged. This may create new nondivisible entries in the row, so keep interchanging the role of rows and columns till no such nondivisible entries remain.

After this, a_{11} divides every entry in its row and column. Subtract suitable multiples of the first column (row) from the other columns (rows) to replace each entry in the first row (column), except a_{11}, by zero. Then we have the correct shape. If the divisibility condition is satisfied, we have finished. If not, there is an entry a_{ij} such that a_{11} does not divide a_{ij} : then add the ith row to the first row and return to the while statement. (In practice many algorithms simply compute a diagonal form first, sorting out the divisibility along the diagonal as a final step. In such algorithms this process is delayed to the end.) Finally reduce A′ recursively.

The algorithm outlined above does not specify how pivots should be chosen, one feature that is exploited by the heuristics presented in [HHR93]. The general idea there is to consider various mathematical combinations of row and column metrics associated with potential pivots. Thus, pivots are chosen which minimize a combination of associated row and column metrics. Each combination can be interpreted as an estimate of some local property of the matrix.

In the following sections we investigate a different aspect of the basic algorithm. Namely we assume that we have selected the pivot, and now we want to force a number of entries of A to 0. As we will see this operation has certain flexibility, which, when properly exploited, can lead to optimal or close to optimal behavior of the algorithm.

3 Macro-step analysis

Our primary goal in Smith normal form computations is to minimize the size of the intermediate entry with maximum magnitude. Here we choose to perform an easier task, where we minimize the size of the intermediate entries during each macro-step. By a macro-step we mean all of the operations required to force all entries in a row and column to zero, except for the entry on the main diagonal. We need to execute min(m, n) macro-steps to compute the Smith normal form of a matrix. Minimizing the size of the entries during each macro-step does not guarantee global optimality, since it is conceivable that the calculation should not proceed macro-step by macrostep. However we will see that even this simplified problem resists good polynomial time solutions.

While investigating algorithms for computing the Smith normal form we have observed that a method with appalling performance is one that looks only at two vectors at the time. Given two columns a_{∗1} and a_{∗j} the method replaces the first column by x_{1}a_{∗1} + x_{2}a_{∗j} and the second column by a_{1j} / gcd(a_{11}, a_{1j} )a_{∗1} − a_{11}/ gcd(a_{11}, a_{1j} )a_{∗j} . The multipliers x_{1} and x_{2} are obtained from gcd(a_{11}, a_{1j} ) = x_{1}a_{11} + x_{2}a_{1j} . A quick analysis reveals the reason why the performance of this method tends to be bad. If the first q entries, a_{11}, …, a_{1q} are such that gcd(a_{11}, …, a_{1p−1}) > gcd(a_{11}, …, a_{1p}) ≠ a_{1p}, for all p ≤ q, then the first column after q steps will be

a_{∗1} ∏_{k=1}^{q−1} x_{2k−1} + ∑_{p=2}^{q} (a_{∗p}x_{2(p−1)} ∏_{k=p}^{q−1} x_{2k−1})

where the atomic multipliers are computed from the equality: x_{2k−1} gcd(a_{11}, …, a_{1p−1})+ x_{2k}a_{1p} = gcd(a_{11}, …, a_{1p}). The cumulative multipliers (given by the product of the atomic multipliers), as indicated in [Bra70], tend to be rather large. Consequently, the expected size of the entries in the first column after the first q elements are forced to zero is substantial. This naturally has a bad impact on the remaining reduction, leading possibly to an unwanted explosion in the size of the intermediate entries. Observe that if a_{11} and a_{21} are relatively prime (as occurs with probability more than 60% for random integers [Knu73]), the average size of the multipliers is O(a_{21}) and O(a_{11}). As a consequence some entries in the transformed first row can be as large as O(a_{11}a_{21}). Using the transformed row in subsequent transformations may cause a more than linear increase in the size of new entries, and consequently start the snowballing effect.

Notice that the above method suffers from the narrowness of its horizon. Obtaining optimal multipliers, that is multipliers such that |x_{1}| ≤ |a_{1j} |/2 and |x_{2}| ≤ |a_{11}|/2 is all the algorithm can do to minimize the intermediate entries in matrix A. This does not provide enough freedom and, consequently, the size of the intermediate entries grows very fast. The conclusion we draw from the above is that forcing any entry to 0 should be done in a wider context, possibly taking into account all vectors of A.

In this analysis we use some well-performing heuristics for selecting pivot elements. The result of Rose and Tarjan [RT78] shows that finding the best method of pivot selection in Gaussian elimination to minimize fill-in is already an NP-complete problem. It follows that pivot selection to minimize entry size in Smith normal form computation, which is achieved by a form of Gaussian elimination, will be difficult.

In each macro-step, one column and row are cleared, i.e., all entries in the column and row (except for the entry on the main diagonal) are forced to zero. The remaining nonzero entry must, once this step is complete, be the greatest common divisor of all other entries in the row and column at the beginning of the step. Hence we consider the following optimization problem:

We are given an integer matrix

B = [ \begin{matrix}b_{11} & … & b_{1n} \\ ⋮ & ... & ⋮ \\ b_{m1} & … & b_{mn} \end{matrix} ] = [b_{∗1}, b_{∗2}, …, b_{∗n}]

where b_{∗i} is the i-th column of B. Find an integer row vector x = [x_{1}, …, x_{m}], with x_{1} ≠ 0, such that

x · b_{∗1} = gcd(b_{11}, …, b_{n1})

minimizing

max_{i,j} |b_{ij} − \frac{x · b_{∗j}}{x · b_{∗1}} b_{i1}| .

By casting the (m − i) × (n − i) submatrix of A into B and solving the above problem we can perform a macro-step in the computation of the Smith normal form of A, achieving optimal or close to optimal performance in terms of the size of the intermediate entries.

An equivalent form of expression for the new value of b_{ij} is

b_{ij} ← ∑_{k=1}^{m} x_{k} \frac{det ([ \begin{matrix}b_{k1} & b_{kj} \\ b_{i1} & b_{ij} \end{matrix} ])}{gcd(b_{11}, …, b_{m1})}.  (1)

Consider now a single new row of B. Denote by d_{ij}^{T} = [d_{ij 1}, …, d_{ijn}] the vector of coefficients associated with the x_{k}’s in (1), d_{ijk} = (b_{k1}b_{ij} − b_{kj} b_{i1})/ gcd(b_{11}, …, b_{m1}). Each b_{ij} can be expressed now as d_{ij} · x. Hence the i-th row of B, b_{i∗}, can be computed as

b_{i∗} = T_{i}x = ∑_{i=1}^{m} t_{∗i}x_{i}  (2)

where T_{i} = [t_{pq}], with t_{pq} = d_{ipq}. The last equation presents a rather unpleasant surprise. The problem of minimizing the maximum entry in b_{i∗} where b_{i∗} is expressed as a linear combination of n vectors has been proved to be NP-complete [vEB81]. In our case the task is even more complicated as we do not have a total freedom in selecting values for the x_{i}’s. Although it is possible to express each x_{i} as ∑_{j} c_{ij} q_{j} , for some constants c_{ij} and integer values q_{j} (cf. [NZM91, §5.2] or [Bla63]), and therefore we can generate any legal set of multipliers, we can not “tune” vector x, by adjusting any single coordinate, without changing the remaining. There is a close relationship between the above optimization problem and the NP-complete CLOSEST VECTOR problem of [vEB81].

An alternative way to approach the problem of keeping the intermediate entries small during Smith normal form computation might be to choose only a small fixed number of vectors k > 2 and use them to perform a single macro-step. Intuitively we must gain something by looking at k > 2 vectors at each step instead of just 2. Unfortunately, we do not know how large k should be. This follows from the fact that the gcd of the leading entries of these k vectors must divide all leading entries of the remaining vectors, which is equivalent to the requirement that the gcd of the leading entries of the selected k vectors must be equal to the gcd of the leading entries of all n vectors. However, even choosing the smallest k with such property is not an easy task. More rigorously we have the following result.

Theorem 1 Given a set of integers U and a positive constant K ≤ |U |, the task of finding a subset R ⊆ U , such that |R| ≤ K and gcd(r ∈ R) = gcd(u ∈ U ) is NP-complete.

Proof. (Sketch) There exists a polynomial transformation from MINIMUM COVER [GJ79, [SP5]] to the above problem. The transformation relies on constructing a set of numbers, such that the i-th number corresponds to the i-th collection in the MINIMUM COVER. The numbers are constructed from the first n prime numbers, in such a way that the j -th prime is a factor of the i-th number if and only if the i-th collection does not contain the j -th member of the set for which the minimum cover is sought. A simple argument proves that finding a set R gives a solution to MINIMUM COVER and vice versa. (For the details of the proof see [HM93].)

4 A greedy approach

Since problem of minimizing the intermediate entries during Smith normal form computations is very likely to be intractable, we redirect our attention to methods, that although they cannot guarantee optimal performance all the time, at least perform reasonably well in practice. As a starting point we select the approach described in [HHR93].

Firstly we notice that a single macro-step can be slightly modified, generally improving the size of the intermediate entries. While from theoretical point it is quite convenient to consider the process as consisting of two parts: (i) bringing the first entry in the first row to the gcd of the all nonzero entries in the first column and then (ii) using the first row to force the leading entries to zero, in practice this approach requires a gcd method that computes very small multipliers. Observe that in equation (1) the size of a new entry in B depends directly on the size of the multipliers. In a separate paper we prove that minimizing the size of the multipliers is intractable [HM93]. Consequently this macro-step structure leads directly to hard problems. A simple alteration however brings quite substantial benefits. Consider Blankinship’s method [Bla63]. At each step the operator row is used to reduce only one row. If, instead, the appropriate multiples of the operator row are subtracted from all other rows the resulting matrix has a much smaller entries.

Example. Consider the following set {f_{n}, f_{n+1}, f_{n+2} − 1}, where f_{i} is the i-th Fibonacci number. The matrix given by Blankinship’s algorithm for such a set is

1  f_{n−1}  f_{n−2}  0
0  f_{n+1}  f_{n}  0
0 (f_{n+2} − 1)f_{n−1} (f_{n+2} − 1)f_{n−2} 1

Observe that entries in the third row are quite large. Now suppose we perform the same basic sequence of the operations, however we subtract the operator row from both other rows. (During Blankinship’s odd numbered steps we subtract the first row from the second and third, and during the unchanged even numbered steps we subtract the second row from the first.) As the result we obtain

1 f_{n−1}  f_{n−2} 0
0 f_{n+1} f_{n} 0
0 f_{n+1} + 1 f_{n−1} − 1 1

Notice the big reduction in the size of the entries in the last row. Even better performance is obtained when we use the best remainder strategy ([HHR93]). Then we obtain

−1  −1  −1  1
0 f_{n−4} − 5 f_{n−4} + 3 −f_{n−4}
0 f_{n−2} − 2 f_{n−2} + 1 −f_{n−2}

The advantage of using the operator row on all rows instead of just one is quite pronounced. ⊓ ⊔

The essence of the above example is more formally expressed by the following lemma:

Lemma 2 There is a modification of Blankinship’s algorithm such that, for a given set of n integers {a_{1}, … a_{n}}, the largest entry in the final matrix does not exceed max{a_{i} + a_{j} }.

Proof. (Sketch) For simplicity we assume that a_{1} and a_{2} are relatively prime, and we observe the changes which occur in only one additional row of matrix C , say the 3rd row. Consider the following modification. Initially we set an auxiliary variable t to ⌊a_{3}/a_{2}⌋ and set a_{3} to a_{3} mod a_{2}. We select row 1 as the operator row, and apply it to both other rows. As a consequence of the initial reduction, a_{3} ≤ a_{2} − 1 and |c_{31}| ≤ |c_{21}| (and also 0 = |c_{32}| ≤ |c_{22}| = 1). In the next step, row 2 is selected as the operator row and is applied to row 1 and 3. Again we have a_{3} ≤ a_{1} − 1 and therefore after this step |c_{31}| ≤ |c_{11}| and |c_{32}| ≤ |c_{12}|. The process continues, with row 1 being the operator row in every odd step and row 2 being the operator row in every even step. Regardless which of the two rows is the operator row, the value of a_{3} is always strictly less than the other number. From that and the properties of the algorithm, it follows (by induction on the number of steps) that the above inequalities hold in odd and even steps, respectively. Once the final matrix is computed we need to correct c_{32}, as initially we reduced a_{3} with a_{2}. Thus c_{32} := c_{32} − t. As we assumed that a_{1} and a_{2} are coprime, we infer that the absolute values of the entries in either row 1 or 2 are equal to a_{2} and a_{1}. As the final entries in the third row are not greater than the entries in the row with the larger entries we conclude that max{c_{ij} } ≤ max{a_{1}, a_{2}} + a_{3}. In the general setting we have the result given in the lemma. ⊓ ⊔

If we apply Blankinship’s algorithm to the first column of matrix B (i.e., we set a_{i} = b_{i1}) then the final matrix given by Blankinship’s algorithm, without the first column is a unimodular matrix P with the property that PB has the same shape as B after a single macro-step. Combining all the above together with the results presented in [HHR93] leads to the following heuristic.

Input: an m × n matrix (B). Output: a matrix equivalent to the input matrix with b_{11} ≠ 0 and b_{i1} = b_{1j} = 0, for i = 2, …, m and j = 2, …, n (or the zero matrix).

Method: If all elements of B are equal to zero then stop. Otherwise choose such a nonzero element of the matrix B for which the product of the Euclidean column and row norms is minimal. If there is more than one such element, choose one with minimal absolute value. Move the selected element to position (1, 1) by row and column interchanges. This is a pivot selection method which differs from earlier methods in that the primary emphasis is on the row and column metrics rather than on the pivot size. The aim remains to restrain entry growth, and this emphasis gives better performance. Unit pivots are still implicitly preferred since they contribute less to both the row and column norms than nonunit pivots. Further, this heuristic remains fully functional in the absence of unit pivots.

For each nondiagonal entry in the first column (b_{i1}, i > 1) compute a “leastremainder” r_{i} = b_{i1} − q_{i}b_{11} of this element divided by b_{11}, with integer quotient q_{i} and such that |r_{i}| ≤ b_{11}/2. Repeat this process with the entries in the first row, computing r′_{j} = b_{1j} −q′_{j}b_{11}. Notice that here we implement the suggestion in [HHR93] to choose a least-remainder rather than more simply taking an arbitrary remainder in the larger range [−(b_{11} − 1), b_{11} − 1]. The beneficial effect of this decision in normal gcd calculations was proved by Kronecker [Kro01, p. 118] and is quantified in [Ros92, p. 67]. It extends to our computations in a natural way.

For i = 2, …, m subtract q_{i} times row one from row i, and likewise subtract q′_{j} times column one from column j . If all elements in the first row and column are zero then stop. (Note that unit pivots lead to fast stopping.) Otherwise, search through the first column and row for the best new pivot, in terms of the product of the row and column norms. As before, if there is more than one potential pivot with the same norm product, choose one of minimal magnitude. Move the new pivot to b_{11} and keep repeating the process till it stops.

Notice that this algorithm is symmetric in a sense. Once the new pivot has been selected both row and column reductions are done. However the order in which row and column reductions are executed does not matter, as the results are the same. Thus, for any submatrix consisting of four elements

[ \begin{matrix}b_{11} & b_{1i} \\ b_{j 1} & b_{ji} \end{matrix} ]

for j , i ≥ 2, element b_{ji} is set to b_{ji} − q_{1}b_{1i} − q_{2}r_{1}, if row reduction precedes column reduction, and is set to b_{ji} − q_{2}b_{j 1} − q_{1}r_{2} otherwise, where b_{j 1} = q_{1}b_{11} + r_{1} and b_{1i} = q_{2}b_{11} + r_{2}. This proves that the new value of b_{ji} is the same, regardless of the relative order of the reductions. The correctness of the algorithm follows from the fact that the pivot is selected to be nonzero and the algorithm has to stop after at most log b_{11} reductions. (The new pivot is always selected from elements that have no more than half the magnitude of the old pivot.)

5 Performance

We have undertaken extensive analyses of the behaviour of the new algorithm. To exemplify this we start by using test data described in detail in [HHR93]. This allows us to compare readily the performance of the new heuristics with previous methods.

Our first example is matrix R_{1}, a 26 × 27 matrix. In the studies included in [HHR93] various maximum entries are reported, depending on strategy used. Thus after 12 macro-steps using the pivot on maximum strategy the largest entry had 1626 decimal digits. The computation completed with a 265 digit maximum entry when choosing to pivot on the first nonzero entry. This improved to 110 digits when pivoting on the minimum entry. The polynomial time Kannan-Bachem [KB79] algorithm produced an 11 digit maximum entry. The metric-based heuristics introduced in [HHR93] achieved maxima between 8501 and 800722. Our new heuristics outperform the best of these, producing a maximum entry 7269.

Our second example is matrix A_{1}, a 304 × 153 matrix, sparse and with small entries. In this case the Kannan-Bachem algorithm produced an 35 digit maximum entry and the metric-based strategies only sometimes succeeded in 32 bit precision. The best result without special reduction techniques had a maximum entry 1895. Here the new method achieved a maximum of 190. In the past this kind of result had been only ever achieved with very substantial use of lattice basis reduction.

Our third example comes from the index 120 subgroup in the Heineken group. Even with metric-based strategies the previous methods all failed without extra lattice basis reduction steps. The new method succeeds with maximal magnitude entry 10606.

In all of these cases the new heuristics outperformed the best heuristic method presented in [HHR93]. While the comparison with the best strategy of Havas, Holt and Rees is not always very impressive, this greedy approach definitely improves on their ×_{2} method, which uses the same norm to control its decisions. Naturally, we can expect, as for any heuristic, that for certain inputs the choice of a different metric could lead to an improvement. However, in general, deciding the best metric a priori is difficult. This is the reason we use only one measure, regardless of the input.

In addition we tested our new method on a number of artificial matrices generated in the following way. First a diagonal matrix with desired Smith normal form is created. Next we use a scrambling program which with probability 0.6 negates a randomly selected row and with probability 0.4 adds to a randomly selected row a constant multiple of another randomly selected row. The constant is randomly selected in the range [−3, +3]. The number of times these operations are performed on input matrix is randomly selected to be between 97 and 154. (The choice of parameters here was dictated by the desire to avoid overflow of 32 bit integers during the scrambling operations.) Matrices created this way were then used as input to the new method. Observe that, as we used only row operations in creating these matrices, it should be (and often is) easier to obtain the Smith normal form by performing row operations before any column operations. However the new method, as it is designed for general case use, does not take advantage of this.

One indicative example was provided by a 10 × 10 matrix, called A_{21}, whose Smith normal form has the following diagonal entries: 1, 1, 1, 1, 6, 6, 60, 60, 60,

180. Only four of the ten diagonal entries are equal to 1, which means that any Smith normal form computation will run out of unit pivots after the first 4 macrosteps. The scrambled matrix had an initial maximum magnitude entry −10620. The new method applied to A_{21} gives a maximum magnitude entry of 13320 after the first macro-step. After the third macro-step the maximum increases to 19140. The remaining steps do not introduce any larger entries and 19140 is the largest intermediate entry generated by the algorithm. Havas, Holt and Reese [HHR93] noted that a potential candidate for proving exponential entry growth is the heuristic that always selects the maximum element as the next pivot. This heuristic, when run on A_{21}, generates a stunning maximum entry that was 11310 decimal digits long. Pivoting on the first nonzero entry generates a largest entry which was 30 digits long. The Kannan-Bachem algorithm generates a nine decimal digit number. The ×_{2} strategy of [HHR93] gives a six digit maximal magnitude entry, 101880.

We also discovered that, for a sample of randomly generated matrices (which frequently have a Smith normal form consisting of m − 1 units followed by a large number), the new heuristics never led to an entry larger than the largest number in the Smith normal form. Thus, in many cases we have observed the computation proceeding close to optimally.

Finally we point out that this method still has its limitations. We have so far not successfully handled a matrix arising from an index 960 subgroup of the Heineken group. This computation can be completed using modular methods. It remains to be seen if we can do so with integer methods by adding appropriate lattice basis reduction steps. While reduction with guaranteed performance (cf. [LLL82]) might be too expensive to use in practice, it is reasonable to hope that some simpler but faster methods may help us reach the goal with acceptable performance.

It is worth noting that the Kannan-Bachem [KB79] algorithm has guaranteed polynomial time complexity as a result of performing similar reductions after forcing a number of entries to zero. However, their reductions are performed so often that they slow down their method significantly, both in practice and in formal complexity terms. An alternative approach is to use some triggering mechanism, which starts reductions every time the intermediate entries become too large, along lines discussed in [HHR93].

6 Concluding remarks

We have presented an new algorithm for computing the Smith normal form of an integer matrix. The algorithm is based on heuristics which reduce intermediate entry growth in Gaussian elimination over the integers. Practical tests indicate very good performance, outperforming previous integer methods on a wide range of examples.

We have also tested these ideas in other contexts. They have shown general applicability elsewhere. Similar methods reduce the difficulty of Hermite normal form computation for integer matrices and exact rational Gaussian elimination. We expect that the principles will extend to many exact matrix computation problems.

References

[Bla63] W.A. Blankinship. A new version of the Euclidean algorithm. Amer. Math. Mon., 70:742–745, 1963.

[Bra70] G.H. Bradley. Algorithm and bound for the greatest common divisor of n integers. Communications of the ACM, 13:433–436, 1970.

[Fru76] M.A. Frumkin. An application of modular arithmetic to the construction of algorithms for solving systems of linear equations. Soviet Math. Dokl., 17:1165–1168, 1976.

[GJ79] M.R. Garey and D.S. Johnson. Computers and Intractibility: A Guide to the Theory of NP-completness. W.H. Freeman, San Francisco, 1979.

[HHR93] G. Havas, D.F. Holt, and S. Rees. Recognizing badly presented Z -modules. Linear Algebra and its Applications, 192:137–163, 1993.

[HM93] G. Havas and B.S. Majewski. On the greatest common divisor of sets of integers. Technical Report TR0288, The University of Queensland, Brisbane, 1993.

[KB79] R. Kannan and R. Bachem. Polynomial time algorithms for computing Smith and Hermite normal forms of an integer matrix. SIAM J. Comput., 8:499–507, 1979.

[Knu73] D.E. Knuth. The Art of Computer Programming, Vol. 2: Seminumerical Algorithms. Addison-Wesley, Reading, Mass., 2nd edition, 1973.

[Kro01] L. Kronecker. Vorlesung über Zahlentheorie I. Teubner, Liepzig, 1901.

[LLL82] A.K. Lenstra, H.W. Lenstra Jr., and L. Lovász. Factoring polynomials with rational coefficients. Math. Ann., 261:515–534, 1982.

[NZM91] I. Niven, H.S. Zuckerman, and H.L. Montgomery. An Introduction to the Theory of Numbers. John Wiley & Sons, Inc., 5th edition, 1991.

[Ros92] K.H. Rosen. Elementary number theory and its applications. Addison-Wesley, Reading, Mass., 3rd edition, 1992.

[RT78] D.J. Rose and R.E. Tarjan. Algorithmic aspects of vertex elimination on directed graphs. SIAM J. Appl. Math., 34:176–197, 1978.

[Smi61] H.J.S. Smith. On systems of linear indeterminate equations and congruences. Philos. Trans. Royal Soc. London, cli:293–326, 1861. Reprinted in The Collected Mathematical Papers of Henry John Stephen Smith, Volume 1, New York: Chelsea (1965).

[vEB81] P. van Emde Boas. Another NP-complete partition problem and the complexity of computing short vectors in a lattice. Technical Report MI/UVA 81–04, The University of Amsterdam, Amsterdam, 1981.