Generation of unimodular matrices with bounded elements Does anybody know what is the algorithm for generating random unimodular matrices (integer matrices with determinant $\pm 1$) whose elements do not exceed a given bound? Such an algorithm is mentioned here, and the following reference is provided: 

Jürgen Hausen, Generating Problems in Linear Algebra, MapleTech, Vol. 1, No. 2, 1994.

However, this paper seems to be no longer accessible online. If the algorithm is based on the Hermite normal form, how do they ensure that the elements of the generated matrix are bounded by a given positive integer? Many thanks in advance for any insights :-)
 A: I tried the following approach


*

*generate a unimodular matrix $A$ (possibly with large entries)

*fix it to reduce the abs size of the entries


In 1, I generate $n\times n$ matrices $L$ and $U$ with 1 on the diagonal; multiply them and randomly permute rows and columns.
In 2, repeat until "convergence" the following step. Pick the largest element in absolute value, assume it's $a_{ij}$. Pick a random non zero entry $a_{kj}$ on the $j$-th column, $k\ne i$, to reduce the size of $a_{ij}$ with $a'_{i,.} = a_{i,.} + a_{k,.}$.
genunimod <- function(n,k=5){
r <- trunc(sqrt(k))+1
A <- matrix(sample(-r:r,n**2,rep=T),n,n)
L <- matrix(0,n,n)
U <- matrix(0,n,n)
l <- outer(1:n,1:n,">")
u <- outer(1:n,1:n,"<")
L[l] <- A[l]
U[u] <- A[u]
diag(U) <- 1
diag(L) <- 1

p1 <- sample(1:n) #permutation shuffle rows
p2 <- sample(1:n) #permutation shuffle cols
P <- diag(rep(1,n)) # identity

A <- P[p1,] %*% L %*% U %*% P[p2,]
# this A is unimodular but may have large entries

m <- k+1
useless <- 0
steps <- 0
while(m>k & useless<=n**2){
    steps <- steps+1
    mm <- max(abs(A))
    v <- which(abs(A)==mm,arr.ind=T)[1,]        
    # use a random element to reduce element v by row
    k <- sample(which(A[,v[2]]!=0)[-v[1]],1)
    s <- sign(A[v[1],v[2]]*A[k,v[2]])
    A[v[1],] <- A[v[1],]-s*A[k,]
    m <- max(abs(A)) # new max term
    if(m>=mm) useless <- useless+1
    }
list(A=A,LU=P[p1,] %*% L %*% U %*% P[p2,],steps=steps,useless=useless)
}

The code can fail in several ways: the largest element can, in principle, be located in an otherwise null column; a single step can be unsuccessful in reducing the largest entry and I stop computations if more than $n^2$ such useless steps occur.
Efficiency can be vastly improved using less primitive methods (say, using "pivoting" ideas instead of randomly picking $k$, using columns and rows and so on).
It takes 1.2 seconds (on a MacBook Air) to generate 1000 $5\times5$ unimodular matrices with entries $\leq 10$. I never had an "otherwise zero" column (good!) and experienced 69 "failures" (7\%, less good!) due to too many useless steps. Hence, 1.2 seconds were indeed needed to create 931 matrices. The digits have the following bell-shaped distribution. It should be symmetric and more weight is given to small numbers.

A: This isn't a complete answer, but it's a start.  The citation is wrong as given, in that the author's name is misspelled.  See here (go to vol.1 no.2), which unfortunately has the desired link broken, but here is a fixed version of that link. Alas, the paper is absent, but its abstract is: 

This article presents a method for creating exercises in linear algebra by making use of Maple's random number generator. The main tool for doing that is a Maple procedure which generates unimodular integral matrices. The presentation also considers standard exercises for basic topics such as linear independence, bases and the rank of a matrix.

Most importantly, the author's name is spelled Hausen (not Hansen).  He is currently a professor at the University of Tübingen.  His website is here, and his email address can be found there.  Although his research interests appear to have departed from Maple programming, I would suggest contacting him directly with your inquiry.
A: First, the unimodular matrix $U$ that you get from doing a Hermite decomposition $A=U R$ ($U$ is unimodular, $R$ is lower-triangular) might have some very large entries. Numerically, I found that for a matrix $A$ of dimension $n\times n$, with its entries chosen uniformly at random from $\{0,\pm1\}$, the maximum entry of $U$ behaves roughly as
$$ \log \max(U) \approx 2(n-10) \log 2, \qquad n\geq 20, $$
so for large dimensions, the straightforward approach doesn't work so well.
Second, if you don't care all that much about getting a specific posterior distribution over the unimodular matrices, there is a very simple algorithm. Let $U$ be an upper-triangular matrix with ones on the diagonal and whose above-diagonal entries are chosen uniformly at random from an integer interval $I$, with maximum magnitude $a = \max\limits_{x\in I}|x|$.
Clearly, $\det U=1$. So, if you form a matrix $A = U_1 U_2^\top$, it will have determinant $1$, and its entries will be at most $(n-1)^2 a^2$. However, its entries will not be uniformly distributed; for instance, the top-left entry will always be $1$. It is also possible to take matrices such as $U_1 U_2^\top U_3 U_4^\top$, which will still have an uneven distribution of entries, but a little less so.
To satisfy an upper bound on matrix elements, it should be enough to pick a sufficiently small value of $a$, and possibly generate several matrices until one of them has sufficiently small elements.
Third, another way is to let $J_k$ be the identity matrix with the nonzero entries in its $k$-th column replaced with random integers. Then $\det J_k=1$ and the matrix $A=J_1J_2\cdots J_n$ will have $\det A=1$, although different columns will have different distributions.
If you wish to have determinant $-1$, choose randomly an odd permutation with a matrix $P$, $\det P=-1$, and use $AP$.
