Generate integer matrices with integer eigenvalues I want to generate $500$ random integer matrices with integer eigenvalues. 
Thanks to this post, I know how to generate a random matrix with whole eigenvalues:


*

*Generate a diagonal matrix $D$ with the desired (integer)
eigenvalues.

*Generate an invertible matrix $A$ of the same size as $D$.

*Record the matrix $A^{-1} D A$.
However, the problem is that this generates a lot of float matrices, and I'd actually like to have both integer matrices and integer eigenvalues. 
The float values are introduced by A.inverse(). According to this post, inverse matrices have whole integer values only when the determinant of the original matrix is 1 or -1 (and therefore an orthogonal matrix).
I tried using the C++ library AlgLib, which has a rmatrixrndorthogonal function that uses G.W. Stewart's 1980 algorithm to generate a random uniformly distributed (Haar) orthogonal matrix. I also tried using R's pracma library, which has a randorthofunction that also generates orthogonal matrices. However, both functions generate matrices with float values. 
Is there a way for me to generate orthogonal, integer matrices?
 A: Let $x$ be a column of our desired orthogonal integer matrix. 
To satisfy the condition that $$\sum_{i=1}^n x_i^2=1$$ where $x_i \in \mathbb{Z}$
we require exactly one $x_i$ to satisfy $|x_i|=1$ and $x_j = 0, j \neq i$.
Hence the columns of such matrices consists of $\{\pm e_1, \ldots, \pm e_n \}$ where $e_i$ is the standard unit vectors.
Remark: 
Suppose $P$ is one such matrix. and suppose $Pe_i = \pm e_j$,
$$e_i^T PDP^Te_i=e_j^TDe_j=d_j$$ 
Hence you are just permutating the diagonal entries. 
Schur decomposition might be of interest to you. 
A: If you only want orthogonal matrices, it's super easy. Here is a uniform sampling procedure for all $n\times n$ integer orthogonal matrices:


*

*Permute the index set $\{1,2,\ldots,n\}$ at random (using, e.g. Fisher-Yates shuffle). Let us call the permuted indices $\{\sigma(1),\sigma(2),\ldots,\sigma(n)\}$.

*Allocate an $n\times n$ zero matrix $A$. Then for $i=1,2,\ldots,n$, set the $(i,\sigma(i))$-th entry (or, the $(\sigma(i),i)$-th entry if the array uses column-major indexing; of course, you also need to make appropriate adjustments if the array uses zero-based indexing rather than one-based indexing) of $A$ to $+1$ or $-1$ at random.

A: Orthogonal matrices with integer values can have only entries $1,0 ,-1$. Every matrix which has only exactly one entry $= +1$ or $-1$ in a column and a row  where the rest entries are zeros is orthogonal.
Easy to check that dot products of columns (or rows) are in this case always equal $0$ and norms for columns (rows) are equal $1$.
The example
$\begin{bmatrix}
 0 &   1 &  0 &  0 \\ 
 -1 &  0 &  0 &  0 \\
0  & 0 &  -1 &  0\\
 0 &  0  & 0 &  1 
\end{bmatrix}$  $\dots$ 
A: Question 1. What do you mean by "I want to generate RANDOM integer matrices with integer eigenvalues" ? 
It seems to be a difficult question. Let $Z$ be the subset of $M_n(\mathbb{Z})$ constituted of matrices that have only integer eigenvalues.
Step 1. It is easy: Let G be an $N\times N$ real matrix whose entries are independent identically distributed
standard normal random variables $G_{i,j}\sim N(0,1)$ (cf. Ginibre 1965).
Step 2. Of course, we can choose a similar discrete probability measure over $\mathbb{Z}$ and randomly construct elements of $M_n(\mathbb{Z})$. It remains to obtain matrices with eigenvalues in $\mathbb{Z}$.
Unfortunately, if $U_N$ is the number of real eigenvalues of a matrix $G$ (step 1), then $E(U_N)=\sqrt{2N/\pi}+O(1)$ when $N\rightarrow +\infty$; here we want $N$ real eigenvalues! Moreover the real eigenvalues of a matrix $G$ cannot be randomly chosen: cf. for example  https://arxiv.org/abs/1512.01449
Then we have two problems: * what is the probability measure induced by $M_n(\mathbb{Z})$ over $Z$ ? 
On the other hand, we can write our matrix $A=P^{-1}DP$ where $D$ is an integer  diagonal matrix (if $A$ has multiple eigenvalues, then we may choose $D$ as a triangular matrix; cf. below); because of what is written above, I think that * the diagonal of $D$ cannot be randomly chosen. 
Question 2. Suppose that we solved the problem of the choice of $D$ in $P^{-1}DP$. Then how to choose $P$ so that $P^{-1}DP\in M_n(\mathbb{Z})$ ? The simplest way is to impose $P\in SL_n(\mathbb{Z})$ (the matrices with $\det=-1$ are useless). Note that the matrix $\begin{pmatrix}-4&8\\-3&6\end{pmatrix}$ cannot be written $P^{-1}diag(0,2)P$ where $P\in SL_2(\mathbb{Z})$ but it can be written $P^{-1}DP$ where $D$ is triangular. Thus, it is better to choose $D$ triangular; I don't know if every matrix in $Z$ can be  written $P^{-1}DP$ where $P\in SL_n(\mathbb{Z})$ and $D$ triangular in $M_n(\mathbb{Z})$. If by extraordinary, the answer is yes, then it remains to randomly obtain an element $P\in SL_n(\mathbb{Z})$.


*

*$SL_n(\mathbb{Z})$ is generated by the transvections. Then we can randomly choose transvections $I_N\pm E_{i,j}$ and consider their product; I think that if $N$ is large, then we must consider many transvections.


*We can also write $A=L_1U_1L_2U_2\cdots$ where the $L_i$ (resp. the $U_i$) are lower-triangular (resp. upper-triangular) with diagonal elements $\pm 1$; I don't think that this randomization is very good but it is easy to write!
